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