copied scripts from passing
[openmx:openmx.git] / demo / BivariateSaturated_MatrixCov.R
1 require(OpenMx)
2
3 #Simulate Data
4 require(MASS)
5 set.seed(200)
6 rs=.5
7 xy <- mvrnorm (1000, c(0,0), matrix(c(1,rs,rs,1),2,2))
8 testData <- xy
9 selVars <- c('X','Y')
10 dimnames(testData) <- list(NULL, selVars)
11 summary(testData)
12 cov(testData)
13
14 #example 3: Saturated Model with Cov Matrices and Matrix-Style Input
15 bivSatModel3 <- mxModel("bivSat3",
16     mxMatrix(
17         type="Symm", 
18         nrow=2, 
19         ncol=2, 
20         free=T, 
21         values=c(1,.5,1), 
22         dimnames=list(selVars,selVars), 
23         name="expCov"
24     ),
25     mxData(
26         observed=cov(testData), 
27         type="cov", 
28         numObs=1000 
29     ),
30     mxMLObjective(
31         covariance="expCov"
32     )
33     )
34 bivSatFit3 <- mxRun(bivSatModel3)
35 EC3 <- mxEvaluate(expCov, bivSatFit3)
36 LL3 <- mxEvaluate(objective,bivSatFit3)
37 SL3 <- summary(bivSatFit3)$SaturatedLikelihood
38 Chi3 <- LL3-SL3
39
40 #example 3m: Saturated Model with Cov Matrices & Means and Matrix-Style Input
41 bivSatModel3m <- mxModel("bivSat3m",
42     mxMatrix(
43         type="Symm", 
44         nrow=2, 
45         ncol=2, 
46         free=T, 
47         values=c(1,.5,1), 
48         dimnames=list(selVars,selVars), 
49         name="expCov"
50     ),
51     mxMatrix(
52         type="Full", 
53         nrow=1, 
54         ncol=2, 
55         free=T, 
56         values=c(0,0), 
57         dimnames=list(NULL,selVars), 
58         name="expMean"
59     ),
60     mxData(
61         observed=cov(testData), 
62         type="cov", 
63         numObs=1000, 
64         means=colMeans(testData) 
65     ),
66     mxMLObjective(
67         covariance="expCov",
68         means="expMean"
69     )
70     )
71 bivSatFit3m <- mxRun(bivSatModel3m)
72 EM3m <- mxEvaluate(expMean, bivSatFit3m)
73 EC3m <- mxEvaluate(expCov, bivSatFit3m)
74 LL3m <- mxEvaluate(objective,bivSatFit3m)
75 SL3m <- summary(bivSatFit3m)$SaturatedLikelihood
76 Chi3m <- LL3m-SL3m
77
78
79 #Mx answers hard-coded
80 #example Mx..1: Saturated Model with Cov Matrices
81 Mx.EC1 <- matrix(c(1.0102951, 0.4818317, 0.4818317, 0.9945329),2,2)
82 Mx.LL1 <- -2.258885e-13
83
84 #example Mx..1m: Saturated Model with Cov Matrices & Means
85 Mx.EM1m <- matrix(c(0.03211648, -0.004883811),1,2)
86 Mx.EC1m <- matrix(c(1.0102951, 0.4818317, 0.4818317, 0.9945329),2,2)
87 Mx.LL1m <- -5.828112e-14
88
89
90 #Compare OpenMx results to Mx results (LL: likelihood; EC: expected covariance, EM: expected means)
91 #3:CovMat
92 omxCheckCloseEnough(Chi3,Mx.LL1,.001)
93 omxCheckCloseEnough(EC3,Mx.EC1,.001)
94 #3m:CovMPat 
95 omxCheckCloseEnough(Chi3m,Mx.LL1m,.001)
96 omxCheckCloseEnough(EC3m,Mx.EC1m,.001)
97 omxCheckCloseEnough(EM3m,Mx.EM1m,.001)