Deleting unrepeatable threshold example, adding 3 repeatables: 3, 5 & 8 variables
[openmx:openmx.git] / models / nightly / thresholdModel1Factor5Variate.R
1 #
2 #  OpenMx Ordinal Data Example
3 #  Revision history:
4 #               Michael Neale 14 Aug 2010
5 #
6
7 # Step 1: load libraries
8 require(OpenMx)
9 require(MASS)
10 #
11 # Step 2: set up simulation parameters 
12 # Note: nVariables>=3, nThresholds>=1, nSubjects>=nVariables*nThresholds (maybe more)
13 # and model should be identified
14 #
15 nVariables<-5
16 nFactors<-1
17 nThresholds<-3
18 nSubjects<-500
19 isIdentified<-function(nVariables,nFactors) as.logical(1+sign((nVariables*(nVariables-1)/2) -  nVariables*nFactors + nFactors*(nFactors-1)/2))
20 # if this function returns FALSE then model is not identified, otherwise it is.
21 isIdentified(nVariables,nFactors)
22
23 loadings <- matrix(.7,nrow=nVariables,ncol=nFactors)
24 residuals <- 1 - (loadings * loadings)
25 sigma <- loadings %*% t(loadings) + vec2diag(residuals)
26 mu <- matrix(0,nrow=nVariables,ncol=1)
27 # Step 3: simulate multivariate normal data
28 set.seed(1234)
29 continuousData <- mvrnorm(n=nSubjects,mu,sigma)
30
31 # Step 4: chop continuous variables into ordinal data 
32 # with nThresholds+1 approximately equal categories, based on 1st variable
33 quants<-quantile(continuousData[,1],  probs = c((1:nThresholds)/(nThresholds+1)))
34 ordinalData<-matrix(0,nrow=nSubjects,ncol=nVariables)
35 for(i in 1:nVariables)
36 {
37 ordinalData[,i] <- cut(as.vector(continuousData[,i]),c(-Inf,quants,Inf))
38 }
39
40 # Step 5: make the ordinal variables into R factors
41 ordinalData <- mxFactor(as.data.frame(ordinalData),levels=c(1:(nThresholds+1)))
42
43 # Step 6: name the variables
44 fruitynames<-paste("banana",1:nVariables,sep="")
45 names(ordinalData)<-fruitynames
46
47
48 thresholdModel <- mxModel("thresholdModel",
49         mxMatrix("Full", nVariables, nFactors, values=0.2, free=TRUE, lbound=-.99, ubound=.99, name="L"),
50         mxMatrix("Unit", nVariables, 1, name="vectorofOnes"),
51         mxMatrix("Zero", 1, nVariables, name="M"),
52         mxAlgebra(vectorofOnes - (diag2vec(L %*% t(L))) , name="E"),
53         mxAlgebra(L %*% t(L) + vec2diag(E), name="impliedCovs"),
54         mxMatrix("Full", 
55             name="thresholdDeviations", nrow=nThresholds, ncol=nVariables,
56             values=.2,
57             free = TRUE, 
58             lbound = rep( c(-Inf,rep(.01,(nThresholds-1))) , nVariables),
59             dimnames = list(c(), fruitynames)),
60     mxMatrix("Lower",nThresholds,nThresholds,values=1,free=F,name="unitLower"),
61     mxAlgebra(unitLower %*% thresholdDeviations, name="thresholdMatrix"),
62             mxFIMLObjective(covariance="impliedCovs", means="M", dimnames = fruitynames, thresholds="thresholdMatrix"),
63             mxData(observed=ordinalData, type='raw')
64 )
65
66 summary(thresholdModelrun <- mxRun(thresholdModel))