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)} }