Compute condition number of information matrix
[openmx:openmx.git] / models / nightly / 3LatentMultiRegWithContinuousModerator-e.R
1 # ---------------------------------------------------------------------
2 # Program: 3LatentMultiRegWithModerator100521.R
3 #  Author: Steven M. Boker
4 #    Date: Sat May 22 11:13:51 EDT 2010
5 #
6 # This program tests variations on a latent variable multiple regression
7 #    using a standard RAM.
8 #
9 # ---------------------------------------------------------------------
10 # Revision History
11 #    -- Sat May 22 11:13:51 EDT 2010
12 #      Created 3LatentMultiRegWithModerator100521.R.
13 #
14 # ---------------------------------------------------------------------
15
16 # ----------------------------------
17 # Read libraries and set options.
18
19 library(OpenMx)
20
21 options(width=100)
22
23 # ---------------------------------------------------------------------
24 # Data for multiple regression of F3 on F1 and F2 with moderator variable Z.
25
26 numberSubjects <- 100000
27 numberIndicators <- 12
28 numberFactors <- 3
29
30 set.seed(10)
31
32 fixedBMatrixF <- matrix(c(.4, .2), 2, 1, byrow=TRUE)
33 randomBMatrixF <- matrix(c(.3, .5), 2, 1, byrow=TRUE)
34 XMatrixF <- matrix(rnorm(numberSubjects*2, mean=0, sd=1), numberSubjects, 2)
35 UMatrixF <- matrix(rnorm(numberSubjects*1, mean=0, sd=1), numberSubjects, 1)
36 Z <- matrix(rnorm(numberSubjects, mean=0, sd=1), nrow=numberSubjects, ncol=2)
37
38 XMatrix <- cbind(XMatrixF, XMatrixF %*% fixedBMatrixF + (XMatrixF*Z) %*% randomBMatrixF + UMatrixF)
39
40 BMatrix <- matrix(c( 1, .6, .7, .8,  0,  0,  0,  0,  0,  0,  0,  0,
41                      0,  0,  0,  0,  1, .5, .6, .7,  0,  0,  0,  0,
42                      0,  0,  0,  0,  0,  0,  0,  0,  1, .7, .6, .5), numberFactors, numberIndicators, byrow=TRUE)
43 UMatrix <- matrix(rnorm(numberSubjects*numberIndicators, mean=0, sd=1), numberSubjects, numberIndicators)
44 YMatrix <- XMatrix %*% BMatrix + UMatrix
45
46 cor(cbind(XMatrix,Z[,1]))
47
48 dimnames(YMatrix) <- list(NULL, paste("X", 1:numberIndicators, sep=""))
49
50 latentMultiRegModerated1 <- cbind(YMatrix,Z=Z[,1])
51
52 round(cor(latentMultiRegModerated1), 3)
53 round(cov(latentMultiRegModerated1), 3)
54
55 latentMultiRegModerated1[,'Z'] <- latentMultiRegModerated1[,'Z'] - mean(latentMultiRegModerated1[,'Z'])
56
57 numberFactors <- 3
58 numberIndicators <- 12
59 numberModerators <- 1
60 indicators <- paste("X", 1:numberIndicators, sep="")
61 moderators <- c("Z")
62 totalVars <- numberIndicators + numberFactors + numberModerators
63
64 # ----------------------------------
65 # Build an orthogonal simple structure factor model
66
67 latents <- paste("F", 1:numberFactors, sep="")
68
69 uniqueLabels <- paste("U_", indicators, sep="")
70 meanLabels <- paste("M_", latents, sep="")
71 factorVarLabels <- paste("Var_", latents, sep="")
72
73 latents1 <- latents[1]
74 indicators1 <- indicators[1:4]
75 loadingLabels1 <- paste("b_F1", indicators[1:4], sep="") 
76 latents2 <- latents[2]
77 indicators2 <- indicators[5:8]
78 loadingLabels2 <- paste("b_F2", indicators[5:8], sep="") 
79 latents3 <- latents[3]
80 indicators3 <- indicators[9:12]
81 loadingLabels3 <- paste("b_F3", indicators[9:12], sep="") 
82
83 threeLatentOrthogonal <- mxModel("threeLatentOrthogonal",
84     type="RAM",
85     manifestVars=c(indicators),
86     latentVars=c(latents,"dummy1"),
87     mxPath(from=latents1, to=indicators1, 
88            arrows=1, connect="all.pairs", 
89            free=TRUE, values=.2, 
90            labels=loadingLabels1),
91     mxPath(from=latents2, to=indicators2, 
92            arrows=1, connect="all.pairs", 
93            free=TRUE, values=.2, 
94            labels=loadingLabels2),
95     mxPath(from=latents3, to=indicators3, 
96            arrows=1, connect="all.pairs", 
97            free=TRUE, values=.2, 
98            labels=loadingLabels3),
99     mxPath(from=latents1, to=indicators1[1], 
100            arrows=1, 
101            free=FALSE, values=1),
102     mxPath(from=latents2, to=indicators2[1], 
103            arrows=1, 
104            free=FALSE, values=1),
105     mxPath(from=latents3, to=indicators3[1], 
106            arrows=1, 
107            free=FALSE, values=1),
108     mxPath(from=indicators, 
109            arrows=2, 
110            free=TRUE, values=.2, 
111            labels=uniqueLabels),
112     mxPath(from=latents,
113            arrows=2, 
114            free=TRUE, values=.8, 
115            labels=factorVarLabels),
116     mxPath(from="one", to=indicators, 
117            arrows=1, free=FALSE, values=0),
118     mxPath(from="one", to=c(latents), 
119            arrows=1, free=TRUE, values=.1, 
120            labels=meanLabels),
121     mxData(observed=latentMultiRegModerated1, type="raw")
122     )
123
124 threeLatentOrthogonalOut <- mxRun(threeLatentOrthogonal)