copied scripts from passing
[openmx:openmx.git] / demo / DefinitionMeans_PathRaw.R
1 require(OpenMx)\r
2 library(MASS)\r
3 #Definition Variable Test 3\r
4 #Author: Mike Neale\r
5 #Date: July 29 2009\r
6 \r
7 #This script is used to test the definition variable functionality in OpenMx\r
8 #The definition variable in this example is dichotomous, and describes two different groups\r
9 #These two groups are measured on two variables, x and y\r
10 #The group with a definition value of 1 has means of 1 and 2 for x and y\r
11 #The group with a definition value of 0 has means af zero for x and y \r
12 #The definition variable is then used to define a mean deviation of the group with definition value 1\r
13 \r
14 #make some data!\r
15 set.seed(200)\r
16 n = 500\r
17 \r
18 Sigma <- matrix(c(1,.5,.5,1),2,2)\r
19 group1<-mvrnorm(n=n, c(1,2), Sigma)\r
20 group2<-mvrnorm(n=n, c(0,0), Sigma)\r
21 \r
22 #put them both together, add a definition variable, and make an selection variables object\r
23 y<-rbind(group1,group2)\r
24 dimnames(y)[2]<-list(c("x","y"))\r
25 def<-rep(c(1,0),each=n)\r
26 selvars<-c("x","y")\r
27 #write data to a file for the mx script to read (not necessary for running in R)\r
28 write.table(cbind(y,def),file="temp-files/xydefmeans.rec",col.names=F,row.names=F)\r
29 \r
30 # Three covariance model matrices: \r
31 #  "cov" for the zero relationship group\r
32 #  "def" for the definition variable, \r
33 #   and "beta" for estimating difference between groups' covariances\r
34 # One common mean vector, "M"\r
35 \r
36 #define the model, including a FIML objective function, which will optimize the matrix S\r
37 defMeansModel<-mxModel("Definition Means via Paths", \r
38     type="RAM",\r
39     mxData(\r
40         observed=data.frame(y,def), \r
41         type="raw"),\r
42     manifestVars=c("x","y"),\r
43     latentVars="DefDummy",\r
44     mxPath(from=c("x","y"), \r
45         arrows=2,\r
46         free=TRUE,\r
47         values=c(1,.1,1),\r
48         labels=c("Varx","Vary")\r
49     ), # variances\r
50     mxPath(from="x", to="y",\r
51         arrows=2,\r
52         free=TRUE,\r
53         values=c(.1),\r
54         labels=c("Covxy")\r
55     ), # covariances\r
56     mxPath(from="one",\r
57         to=c("x","y","DefDummy"),\r
58         arrows=1,\r
59         free=c(TRUE,TRUE,FALSE),\r
60         values=c(1,1,1),\r
61         labels =c("meanx","meany","data.def")\r
62     ),\r
63     mxPath(from="DefDummy",\r
64         to=c("x","y"),\r
65         arrows=1,\r
66         free=c(TRUE,TRUE),\r
67         values=c(1,1),\r
68         labels =c("beta_1","beta_2")\r
69     )\r
70 )\r
71 \r
72 # Run the model\r
73 defMeansFit<-mxRun(defMeansModel)\r
74 defMeansFit@matrices\r
75 defMeansFit@algebras\r
76 \r
77 # Compare OpenMx estimates to summary statistics from raw data, remembering to knock off 1 and 2 from group 1's\r
78 # data, so as to estimate variance of combined sample without the mean correction.\r
79 \r
80 # First we compute some summary statistics from the data\r
81 ObsCovs <- cov(rbind(group1 - rep(c(1,2), each=n), group2))\r
82 ObsMeansGroup1 <- c(mean(group1[,1]), mean(group1[,2]))\r
83 ObsMeansGroup2 <- c(mean(group2[,1]), mean(group2[,2]))\r
84 \r
85 # Second we extract the parameter estimates and matrix algebra results from the model\r
86 Sigma <- mxEvaluate(S[1:2,1:2], defMeansFit)\r
87 Mu <- mxEvaluate(M[1:2], defMeansFit)\r
88 beta <- mxEvaluate(A[1:2,3], defMeansFit)\r
89 \r
90 # Third, we check to see if things are more or less equal\r
91 omxCheckCloseEnough(ObsCovs,Sigma,.01)\r
92 omxCheckCloseEnough(ObsMeansGroup1,as.vector(Mu+beta),.001)\r
93 omxCheckCloseEnough(ObsMeansGroup2,as.vector(Mu),.001)\r
94 \r