############################################################################## # Program: Vergleich PLS und FIML # # copyright: Florian Schuberth # University of Wuerzburg # florian.schuberth@uni-wuerzburg.de ############################################################################### rm(list=ls(all=TRUE)) library(MASS); library(matrixcalc); library(lattice) library(lavaan) library(mvtnorm) library(polycor) library(OpenMx) # number of iterations it<-1 #Matrices where the results are collected #rows of the matrix aa<-19 result<-matrix(nrow=aa, ncol=it) #standardized estimates result_std<-matrix(nrow=aa, ncol=it) rownames(result)<-c("lambda_y_21","lambda_y_31","lambda_x_21","lambda_x_31","beta","Var(Zeta)","Var(Xi)","Thresh_x1_1","Thresh_x1_2","Thresh_x2_1","Thresh_x2_2","Thresh_x3_1","Thresh_x3_2","Thresh_y1_1","Thresh_y1_2","Thresh_y2_1","Thresh_y2_2","Thresh_y3_1","Thresh_y3_2") rownames(result_std)=tau_x1_2) {x1[iii]<-3} } ########### categorical indicator x2 ########## x2_stern<-lambda_x_21*xi1+delta2 #Initialize: x2<-12 for (iii in 1:ii) { if(x2_stern[iii] < tau_x2_1) {x2[iii]<-1} if (tau_x2_1 <= x2_stern[iii] && x2_stern[iii] < tau_x2_2) {x2[iii]<-2} if (x2_stern[iii]>=tau_x2_2) {x2[iii]<-3} } ############# categorical indicator x3 ######### x3_stern<-lambda_x_31*xi1+delta3 x3<-12 for (iii in 1:ii) { if(x3_stern[iii] < tau_x3_1) {x3[iii]<-1} if (tau_x3_1 <= x3_stern[iii] && x3_stern[iii] < tau_x3_2) {x3[iii]<-2} if (x3_stern[iii]>=tau_x3_2) {x3[iii]<-3} } ############# categorical indicator y1 ############### y1_stern<-lambda_y_11*eta1+eps1 y1<-12 for (iii in 1:ii) { if(y1_stern[iii] < tau_y1_1) {y1[iii]<-1} if (tau_y1_1 <= y1_stern[iii] && y1_stern[iii] < tau_y1_2) {y1[iii]<-2} if (y1_stern[iii]>=tau_y1_2) {y1[iii]<-3} } ############# categorical indicator y2 ################ y2_stern<-lambda_y_21*eta1+eps2 y2<-12 for (iii in 1:ii) { if(y2_stern[iii] < tau_y2_1) {y2[iii]<-1} if (tau_y2_1 <= y2_stern[iii] && y2_stern[iii] < tau_y2_2) {y2[iii]<-2} if (y2_stern[iii]>=tau_y2_2) {y2[iii]<-3} } ############# categorical indicator y3 ################ y3_stern<-lambda_y_31*eta1+eps3 y3<-12 for (iii in 1:ii) { if(y3_stern[iii] < tau_y3_1) {y3[iii]<-1} if (tau_y3_1 <= y3_stern[iii] && y3_stern[iii] < tau_y3_2) {y3[iii]<-2} if (y3_stern[iii]>=tau_y3_2) {y3[iii]<-3} } datmat<-data.frame(x1,x2,x3, y1,y2,y3); colnames(datmat)<-c("x1","x2","x3","y1","y2","y3") ###################################################################### #################### openMx Estimation ############# ##### categorical indicator #### y1<-mxFactor(y1, levels=c(1,2,3), ordered=TRUE ) y2<-mxFactor(y2, levels=c(1,2,3), ordered=TRUE ) y3<-mxFactor(y3, levels=c(1,2,3), ordered=TRUE ) x1<-mxFactor(x1, levels=c(1,2,3), ordered=TRUE ) x2<-mxFactor(x2, levels=c(1,2,3), ordered=TRUE ) x3<-mxFactor(x3, levels=c(1,2,3), ordered=TRUE ) datmatmx<-data.frame(x1,x2,x3, y1,y2,y3) #Indicators observed<-c("x1","x2","x3","y1","y2","y3") blockx<-c("x1","x2","x3") blocky<-c("y1","y2","y3") #latent variables latent<-c("eta1","xi1") ###### OpenMx in RAM notation ### model_mx<-mxModel("Title", mxData(observed=datmatmx, type="raw") , type="RAM", manifestVars=observed, latentVars=latent, #Block X mxPath(from=c("xi1"),to=c(blockx),arrows=1,free=c(F,T,T),connect='single', values=c(1, 2.3 ,1.6),labels=c("lambda_x_11","lambda_x_21","lambda_x_31")), #Block Y mxPath(from=c("eta1"),to=c(blocky),arrows=1,free=c(F,T,T), values=c(1,1.5,-4),labels=c("lambda_y_11","lambda_y_21","lambda_y_31")), #pathmodel mxPath(from=c("xi1"),to=c("eta1"),arrows=1,free=T,values=0.6,labels="beta"), #error variances of block X and block Y mxPath(from=observed,arrows=2 , free=c(F,F,F,F,F,F),connect='single' , labels=c("Var(delta1)","Var(delta2)","Var(delta3)","Var(eps1)","Var(eps2)","Var(eps3)"),lbound=0.0001,values=c(1 , 1 , 1 , 1 , 1 , 1)), #error terms of the latent varibale mxPath(from=c("xi1","eta1"), arrows=2 , free=c(T,T), values=c(1.3 , 0.9),labels=c("Var(Xi)","Var(Zeta)"),connect='single' ) , # means mxPath(from="one", to=c(observed,latent),arrows=1,free=c(F,F,F,F,F,F,F,F),values=c(0,0,0,0,0,0,0,0)), #Matrix with the thresholds mxMatrix( type="Full", nrow=2, ncol=6, free=c(T,T,T,T,T,T,T,T,T,T,T,T), values=c(-1 , -0.55 , -1.1 , -0.53 , -0.91 , -0.39 , 0.4 , 1.45 , 0.2 , 1.01 , 1.26, 0.7), dimnames=list(c(),c("x1","x2","x3","y1","y2","y3")), byrow=TRUE, name="thresh" ), mxRAMObjective(A="A", S="S", F="F" , M="M", thresholds="thresh") ) test1<-mxRun(model_mx) ############## Output ############# ausgabe<-summary(test1) result[,i]<-ausgabe$parameters$Estimate result_std[,i]<-ausgabe$parameters$Std.Estimate loglike[,i]<-ausgabe$Minus2LogLikelihood }