1"sim.structure" <- "sim.structural" <-
2function (fx=NULL,Phi=NULL,fy=NULL,f=NULL,n=0,uniq=NULL,raw=TRUE, items = FALSE, low=-2,high=2,d=NULL,cat=5,mu=0) {
3 cl <- match.call()
4require(MASS)
5
6
7 if(is.null(f)) { if(is.null(fy)) {f <- fx} else {
8    f <- super.matrix(fx,fy)} }
9  f <- as.matrix(f)
10  if(!is.null(Phi)) {if(length(Phi)==1) Phi <- matrix(c(1,Phi,Phi,1),2,2)}
11   #these are parameters for simulating items
12   nf <- ncol(f)
13   nvar <- nrow(f)
14if(is.null(d)) {d <- seq(low,high,(high-low)/(nvar/nf-1))
15            d <- rep(d,nf)} else {if(length(d)==1) d <- rep(d,nvar)}
16	a <- rep(1,nvar)
17
18  if(is.vector(f)) {f <- as.matrix(f)  #this is the case if doing a congeneric model
19                    Phi <- 1}
20  if(!is.null(Phi)) {
21  model <-  f %*% Phi %*%  t(f) #the model correlation matrix for oblique factors
22} else { model <- f%*% t(f)}
23if(is.null(uniq)) {diag(model) <- 1 } else { diag(model) <- uniq  + diag(model)}                      # put ones along the diagonal unless uniq is specified
24  nvar <- dim(f)[1]
25  if(is.null(rownames(model))) {colnames(model) <- rownames(model) <- paste("V",1:nvar,sep="")} #else {colnames(model) <- rownames(model) <- rownames(fx)}
26
27  if(n>0) {
28    mu <- rep(mu,nvar)
29  	observed <- mvrnorm(n = n, mu, Sigma=model, tol = 1e-6, empirical = FALSE)
30
31  	if(items) {observedp <- matrix(t(pnorm(a*t(observed)- d)),n,nvar)
32  	         observed[] <- rbinom(n*nvar, cat, observedp)}
33  r <- cor(observed) }
34
35  	reliability <- diag(f %*% t(f))
36  if(n<1) {results <- list(model=model,reliability=reliability) } else {
37  if (!raw) {results <- list( model=model,reliability=reliability,r=r,N=n )} else {
38             results <- list( model=model,reliability=reliability,r=r,observed= observed,N=n) } }
39  results$Call <- cl
40  class(results) <- c("psych", "sim")
41 return(results)}
42
43
44"sim" <-
45function (fx=NULL,Phi=NULL,fy=NULL,alpha=.8,lambda = 0,n=0,mu=NULL,raw=TRUE) {
46 cl <- match.call()
47require(MASS)
48#set up some default values
49if(is.null(fx)) {fx <- matrix(c(rep(c(.8,.7,.6,rep(0,12)),3),.8,.7,.6),ncol=4)
50   if(is.null(Phi)) {Phi <- diag(1,4,4)
51                     Phi <- alpha^abs(row(Phi) -col(Phi)) + lambda^2
52                     diag(Phi) <- max((alpha + lambda),1)
53	                Phi <- cov2cor(Phi)}
54   if(is.null(mu)) {mu <- c(0,.5,1,2)}                   }
55    if(is.null(fy)) {f <- fx} else {
56    f <- super.matrix(fx,fy)}
57
58   if(is.null(mu)) {mu <- rep(0,ncol(fx))}
59
60     means <- fx %*% mu
61
62  if(is.vector(f)) {f <- as.matrix(f)  #this is the case if doing a congeneric model
63                    Phi <- 1}
64  if(!is.null(Phi)) {
65  model <-  f %*% Phi %*%  t(f) #the model correlation matrix for oblique factors
66} else { model <- f%*% t(f)}
67diag(model)<- 1                       # put ones along the diagonal
68  nvar <- dim(f)[1]
69  colnames(model) <- rownames(model) <- paste("V",1:nvar,sep="")
70  if(n>0) {
71
72  	observed <- mvrnorm(n = n, means, Sigma=model, tol = 1e-6, empirical = FALSE)
73  r <- cor(observed) }
74  	reliability <- diag(f %*% t(f))
75  if(n<1) {results <- list(model=model,reliability=reliability) } else {
76  if (!raw) {results <- list( model=model,reliability=reliability,r=r,N=n )} else {
77             results <- list( model=model,reliability=reliability,r=r,observed= observed,N=n) } }
78  results$Call <- cl
79  class(results) <- c("psych", "sim")
80 return(results)}
81
82
83"sim.simplex" <-
84	function(nvar =12, alpha=.8,lambda=0,beta=1,mu=NULL, n=0) {
85	 cl <- match.call()
86	R <- matrix(0,nvar,nvar)
87	R[] <- alpha^abs(col(R)-row(R))*beta + lambda^2
88	diag(R) <- max((alpha * beta) + lambda,1)
89	R <- cov2cor(R)
90	colnames(R) <- rownames(R) <- paste("V",1:nvar,sep="")
91	require(MASS)
92	if(is.null(mu)) {mu <- rep(0,nvar)}
93	if(n>0) {
94	observed.scores <- mvrnorm(n = n, mu, Sigma=R, tol = 1e-6, empirical = FALSE)
95	observed <- cor(observed.scores)
96	results <- list(model=R,r=observed,observed=observed.scores)
97	results$Call <- cl
98	class(results) <- c("psych", "sim")} else {results <- R}
99
100	results
101	}
102
103
104#simulate major and minor factors
105"sim.minor" <-
106function(nvar=12,nfact=3,n=0,fbig=NULL,fsmall = c(-.2,.2),bipolar=TRUE) {
107if(is.null(fbig)) {loads <- c(.8,.6) } else {loads <- fbig}
108loads <- sample(loads,nvar/nfact,replace=TRUE)
109if(nfact == 1) {fx <- matrix(loads,ncol=1)} else {fx <- matrix(c(rep(c(loads,rep(0,nvar)),(nfact-1)),loads),ncol=nfact)}
110if(bipolar) fx <- 2*((sample(2,nvar,replace=TRUE) %%2)-.5) * fx
111fsmall  <- c(fsmall,rep(0,nvar/4))
112fs <- matrix(sample(fsmall,nvar*floor(nvar/2),replace=TRUE),ncol=floor(nvar/2))
113fload <- cbind(fx,fs)
114colnames(fload) <- c(paste("F",1:nfact,sep=""),paste("m",1:(nvar/2),sep=""))
115rownames(fload) <- paste("V",1:nvar,sep="")
116results <- sim(fload,n=n)
117         results$fload <- fload
118 class(results) <- c("psych", "sim")
119 return(results)
120}
121
122
123#similuate various structures and summarize them
124"sim.omega" <-
125function(nvar=12,nfact=3,n=0,fbig=NULL,fsmall = c(-.2,.2),bipolar=TRUE,om.fact=3,flip=TRUE,option="equal",ntrials=10) {
126results <- matrix(NA,nrow=ntrials,ncol=7)
127for (i in 1:ntrials) {
128x <- sim.minor(nvar=nvar,nfact=nfact,n=n,fbig=fbig,fsmall=fsmall,bipolar=bipolar)
129if(n > 0) {
130om <- try(omega(x$observed,om.fact,flip=flip,plot=FALSE,option=option))
131          ic <- ICLUST(x$observed,1,plot=FALSE)} else {om <- try(omega(x$model,om.fact,flip=flip,plot=FALSE,option=option))
132          ic <- ICLUST(x$model,1,plot=FALSE)}
133
134results[i,1] <- om$omega_h
135     loads <- om$schmid$sl
136   p2 <- loads[,ncol(loads)]
137         	mp2 <- mean(p2)
138         	vp2 <- var(p2)
139
140
141results[i,2] <-  mp2
142results[i,3] <-  sqrt(vp2)
143results[i,4] <- sqrt(vp2)/mp2
144results[i,5] <-ic$beta
145if(!is.null(om$schmid$RMSEA)) {results[i,6] <-om$schmid$RMSEA[1]} else {results[i,6] <- NA}
146if(!is.null(om$schmid$rms))results[i,7] <- om$schmid$rms
147
148
149}
150colnames(results) <- c("omega","mean p2","sd p2","coeff v","Beta","RMSEA","rms")
151return(results)
152}
153
154
155
156"sim.general" <-
157function(nvar=9,nfact=3, g=.3,r=.3,n=0) {
158require(MASS)
159  r1 <- matrix(r,nvar/nfact,nvar/nfact)
160  R <- matrix(g,nvar,nvar)
161  rf <- super.matrix(r1,r1)
162  if(nfact>2) {for (f in 1:(nfact-2)){
163  rf <- super.matrix(r1,rf)}}
164  R <- R + rf
165  diag(R) <- 1
166  colnames(R) <- rownames(R) <- paste((paste("V",1:(nvar/nfact),sep="")),rep(1:nfact,each=(nvar/nfact)),sep="gr")
167  if(n > 0) {x <-  mvrnorm(n = n, mu=rep(0,nvar), Sigma = R, tol = 1e-06,empirical = FALSE)
168   return(x)} else {
169  return(R)} }