10001036

#####################
##Population Parameters
#####################
c.ignore= .1; ## c.ignore=0;
set.seed(0);
n.factor=3;
n.load = 4;
beta.ignore = c.ignore;
beta.max=.5;
beta.min=.4;
lambda.max =1.1;
lambda.min =.9;
lambda.fix =1;
lambda.ignore = c.ignore;
theta.ignore = c.ignore/2;
theta.max=0.45
theta.min=0.35
#####################
require(Rpad);require(R2HTML);HTMLon(); ## online output
#####################
#######################
##Population
#######################
Beta = matrix(0,nrow=n.factor, ncol=n.factor);
Beta[row(Beta)-col(Beta)>=1] =
runif(length(Beta[row(Beta)-col(Beta) >= 1]),max=beta.ignore,min=-beta.ignore);
Beta[row(Beta)-col(Beta)==1] =
runif(length(Beta[row(Beta)-col(Beta)==1]),max=beta.max,min=beta.min);
Beta
HTML(Beta);
Psi = matrix(0,nrow=n.factor, ncol=n.factor);
Psi[] =
runif(Psi[],max=beta.ignore,min=-beta.ignore);
Psi[row(Psi)-col(Psi)==0] =
runif(Psi[row(Psi)-col(Psi) == 0],max=1-beta.min^2,min=1-beta.max^2);
Psi[1,1]=Psi[1,1]+mean(c(beta.min,beta.max) ^2);
Psi = .5 * (Psi +t(Psi));
eigen(Psi)$values
Psi
HTML(Psi);
##eta = Beta*eta+zeta
## (I-Beta)*eta=zeta
## eta =(I-Beta)^-1 * zeta
## var(eta)= (I-Beta)^-1 *Psi* t((I-Beta)^-1)
Sigma.eta = solve(diag(n.factor)-Beta) %*%
Psi %*%
t(solve(diag(n.factor)-Beta));
Sigma.eta
HTML(Sigma.eta);
Lambda.y=diag(nrow(Beta)) %x%
matrix(c(lambda.fix,rep(-1,n.load-1)),nrow=n.load,ncol=1);
Lambda.y[0>Lambda.y] =
runif(Lambda.y[0>Lambda.y],min=lambda.min,max=lambda.max);
Lambda.y[Lambda.y==0] =
runif(Lambda.y[Lambda.y==0],min=-lambda.ignore,max=lambda.ignore);
Lambda.y;
HTML(Lambda.y)
Theta.y=diag(nrow(Lambda.y));
Theta.y[Theta.y>0]=
runif(Theta.y[Theta.y>0],min=theta.min,max=theta.max);
Theta.y[Theta.y==0]=
runif(Theta.y[Theta.y==0],min=-theta.ignore,max=theta.ignore);
Theta.y = .5 * (Theta.y +t(Theta.y));
Theta.y
HTML(Theta.y);
Sigma = Lambda.y %*% Sigma.eta %*% t(Lambda.y) + Theta.y;
Sigma = .5 * (Sigma +t(Sigma));
Sigma
HTML(Sigma);
#####################
#####################
######################
##Model
######################
rownames(Sigma) = colnames(Sigma) = paste('Y',1:ncol(Sigma),sep='');
model.r = matrix('',ncol=3,nrow=0);
class(model.r)='mod';
ij = (row(Beta)-col(Beta)==1);
model.Beta = t(rbind(
paste('eta',row(Beta)[ij], ' <- ',
'eta',col(Beta)[ij],sep=''),
paste('beta',row(Beta)[ij],'_',col(Beta)[ij],sep=''),
rep(NA,sum(ij))));
ij<-(row(Psi)==col(Psi));
model.Psi <- t(rbind(
paste('eta',row(Psi)[ij], ' <-> ',
'eta',col(Psi)[ij],sep=''),
paste('psi',row(Psi)[ij],'_',col(Psi)[ij],sep=''),
rep(NA,sum(ij))));
ij<-(Lambda.y==lambda.fix);
model.Lambda.y.fix <- t(rbind(
paste('Y',row(Lambda.y)[ij], ' <- ',
'eta',col(Lambda.y)[ij],sep=''),
rep(NA,sum(ij)),
rep(1,sum(ij))));
ij<-((Lambda.y!=lambda.fix)&(Lambda.y>lambda.ignore));
model.Lambda.y <- t(rbind(
paste('Y',row(Lambda.y)[ij], ' <- ',
'eta',col(Lambda.y)[ij],sep=''),
paste('lambda_y',row(Lambda.y)[ij],'_',col(Lambda.y)[ij],sep=''),
rep(NA,sum(ij))));
ij<-((Theta.y>theta.ignore)&(row(Theta.y)>=col(Theta.y)));
model.Theta.y <- t(rbind(
paste('Y',row(Theta.y)[ij], ' <-> ',
'Y',col(Theta.y)[ij],sep=''),
paste('theta_y',row(Theta.y)[ij],'_',col(Theta.y)[ij],sep=''),
rep(NA,sum(ij))));
model.r <-rbind(model.Beta, model.Psi,
model.Lambda.y.fix, model.Lambda.y,
model.Theta.y);
ij <- (col(Sigma)>=row(Sigma));
model.u <- t(rbind(
paste('Y',col(Sigma)[ij], ' <-> ',
'Y',row(Sigma)[ij],sep=''),
paste('sigma',col(Sigma)[ij],'_',row(Sigma)[ij],sep=''),
rep(NA,sum(ij))));
class(model.u)='mod';
#####################
HTML(model.r);
##HTML(model.u);
#####################
####################
##Modeled Population
####################
require(sem);
sem.Sigma<-sem(model.r,Sigma,101);
Sigma.tilde<-sem.Sigma$C;
theta.tilde<-sem.Sigma$coeff;
##Population modeled parameters
HTML(t(t(theta.tilde)));
##population misfit:
F0.tilde<-summary(sem.Sigma)$chisq /100;
RMSEA.tilde<-(F0.tilde/summary(sem.Sigma)$df)^.5;
HTMLh2(round(F0.tilde,4));
HTMLh2(round(RMSEA.tilde,4));
(res.tilde<-Sigma-Sigma.tilde);
HTML(round(res.tilde,4));
#####################
boxplot.matrix = function(M,ylim=c(-1,1)) {
M = as.matrix(M);
boxplot(c(M[row(M)>col(M)]),at=1,xlab='',ylab='',ylim=ylim);
points(rep(1,length(c(M[row(M)>col(M)]))),c(M[row(M)>col(M)]),pch='-',col='red');
boxplot.stats(c(M[row(M)>col(M)]));
}
boxplot(res.tilde);
####################
#####################
HTMLoff(); ## online output
#####################