library(MASS)
#--------------------------------------------------
#
R script for "Using non-normal SEM to resolve the
#
ACDE model in the classical twin design"
#
(Ozaki, K., Toyoda, H., Iwama, N., Kubo, S.,
& Ando, J. (2011).
# Behavior Genetics.)
#
This enables you to estimate A, C, D, and E effects using non-normal SEM
#
Author: Koken Ozaki, 18 Augst
2016 (version5)
#
Affiliation: Tsukuba University, in Japan.
#
E-mail: koken@gssm.otsuka.tsukuba.ac.jp
#
Special Thanks to Dr.Taylor @Univ of British Columbia
#--------------------------------------------------
################################ Usage
###################################
#You have to specify dataset, usevariables, variable names.
#You should not change the script outside
of "From here" to "Up to here."
#The missing value flag has to be
"."
#estimates and fit results are provided.
###########################################################################
############################## From here
##################################
#data input
data<-read.table("bgdata.txt",
na.strings=".") #you have to change
#specification of usevariables
usevariables<-cbind(data[,1],data[,2],data[,3])
#you have to change
# which indicates that 1st(zygosity),
# 2nd(variable 1 of twin 1), 3rd(variable 1
of twin 2) column of the data
# are used as variables.
# Note that first variable must indicate zygosity.
usevariables<-as.data.frame(usevariables)
#specify MZ data and DZ data
data1<-subset(usevariables,
usevariables[1]==1) #MZ #you may have to change
data2<-subset(usevariables,
usevariables[1]==2) #DZ #you may have to change
# the 1st column of the usevariables
is zygosity
# (in this case MZ=1 and DZ=2)
#variable names
v1.name<-"v1"; #you can change
v2.name<-"v2"; #you can change
#The number of initial value sets in
optimization
initialrep<-2000
############################# Up to here
##################################
library(MASS)
xtm<-data1[2] #variable 1 of twin 1 of MZ
xcm<-data1[3] #variable 1 of twin 2 of MZ
xtd<-data2[2] #variable 1 of twin 1 of DZ
xcd<-data2[3] #variable 1 of twin 2 of DZ
M1<-as.matrix(data.frame(xtm,xcm))
D1<-as.matrix(data.frame(xtd,xcd))
M2<-na.omit(M1)
D2<-na.omit(D1)
xtm<-M2[,1]
xcm<-M2[,2]
xtd<-D2[,1]
xcd<-D2[,2]
#standardization
xtm <- xtm-mean(xtm);
xtm<-xtm/sqrt(mean(xtm*xtm))
xcm <- xcm-mean(xcm);
xcm<-xcm/sqrt(mean(xcm*xcm))
xtd <- xtd-mean(xtd);
xtd<-xtd/sqrt(mean(xtd*xtd))
xcd <- xcd-mean(xcd);
xcd<-xcd/sqrt(mean(xcd*xcd))
M<-as.matrix(data.frame(xtm,xcm))
D<-as.matrix(data.frame(xtd,xcd))
MZrep<-nrow(M);
DZrep<-nrow(D);
final.est1<-matrix(0, ncol=9, nrow=1)
final.fit1<-matrix(0, ncol=1, nrow=1)
final.se21<-matrix(0, ncol=9, nrow=1)
final.est2<-matrix(0, ncol=10, nrow=1)
final.fit2<-matrix(0, ncol=1, nrow=1)
final.se22<-matrix(0, ncol=10, nrow=1)
final.est3<-matrix(0, ncol=10, nrow=1)
final.fit3<-matrix(0, ncol=1, nrow=1)
final.se23<-matrix(0, ncol=10, nrow=1)
#2nd order moments of MZs
m1.2<-mean(M[,1]*M[,1])
m11<-mean(M[,1]*M[,2])
m2.2<-mean(M[,2]*M[,2])
#3rd order moments of MZs
m1.3<-mean(M[,1]*M[,1]*M[,1])
m21<-mean(M[,1]*M[,1]*M[,2])
m12<-mean(M[,1]*M[,2]*M[,2])
m2.3<-mean(M[,2]*M[,2]*M[,2])
#2nd order moments of DZs
d1.2<-mean(D[,1]*D[,1])
d11<-mean(D[,1]*D[,2])
d2.2<-mean(D[,2]*D[,2])
#3rd order moments of DZs
d1.3<-mean(D[,1]*D[,1]*D[,1])
d21<-mean(D[,1]*D[,1]*D[,2])
d12<-mean(D[,1]*D[,2]*D[,2])
d2.3<-mean(D[,2]*D[,2]*D[,2])
#weight matrix
weight.matrix <- function(X){
p<- ncol(X)
X<-t(t(X)-colMeans(X))
############# w22 ##################
uni2 <- matrix(0,p,p)
uni2[lower.tri(uni2,diag=TRUE)]<-1
w22index <- w22index_2 <- which(uni2==TRUE, arr.ind=TRUE)
w22nrow <- nrow(w22index)
w22nele <-
w22nrow*(w22nrow+1)/2
w22index1 <- matrix(0,w22nele,w22nrow)
for(i in 1:w22nrow){
w22index1[(1+i*(i-1)/2):(
(i*(i-1)/2)+i ), i]<- rep(1, i)
}
w22index2_1 <- diag(w22nrow)
w22index2 <- matrix(0, w22nele, w22nrow)
w22cumindex <- cumsum(1:w22nrow)
w22index2[1,1]<-1
for (i in 2: w22nrow){
w22index2[(w22cumindex[i-1]+1):(w22cumindex[i]), ] <- w22index2_1[1:i,]
}
w22index <- cbind(w22index1%*%w22index,w22index2%*%w22index)
w22 <- matrix(0, w22nrow,
w22nrow)
w22ele <- numeric(w22nele)
for (o in 1:w22nele){
i <- w22index[o,1]
j
<- w22index[o,2]
k
<- w22index[o,3]
l
<- w22index[o,4]
w22ele[o]<-(mean(X[,i]*X[,j]*X[,k]*X[,l])
-mean(X[,i]*X[,j])*mean(X[,k]*X[,l]))
}
w22[upper.tri(w22, diag=TRUE)]<-w22ele
############# w33 ##################
uni3 <- array(0,c(p,p,p))
uni3[, , 1][lower.tri(uni3[,,1], diag=TRUE)]<-1
for(i in 2:p ){
uni3[i:p, i:p, i] <- uni3[1:((p+1)-i), 1:((p+1)-i),1]
}
w33index <- w33index_2 <- which(uni3==TRUE, arr.ind=TRUE)
w33nrow <- nrow(w33index)
w33nele <-
w33nrow*(w33nrow+1)/2
w33index1 <- matrix(0,w33nele,w33nrow)
for(i in 1:w33nrow){
w33index1[(1+i*(i-1)/2):(
(i*(i-1)/2)+i ), i]<- rep(1, i)
}
w33index2_1 <- diag(w33nrow)
w33index2 <- matrix(0, w33nele, w33nrow)
w33cumindex <- cumsum(1:w33nrow)
w33index2[1,1]<-1
for (i in 2: w33nrow){
w33index2[(w33cumindex[i-1]+1):(w33cumindex[i]), ] <- w33index2_1[1:i,]
}
w33index <- cbind(w33index1%*%w33index,w33index2%*%w33index)
w33 <- matrix(0, w33nrow,
w33nrow)
w33ele <- numeric(w33nele)
for (o in 1:w33nele){
i <- w33index[o,1]
j
<- w33index[o,2]
k
<- w33index[o,3]
l
<- w33index[o,4]
m
<- w33index[o,5]
n
<- w33index[o,6]
w33ele[o]<-(mean(X[,i]*X[,j]*X[,k]*X[,l]*X[,m]*X[,n])
-mean(X[,m]*X[,n])*mean(X[,i]*X[,j]*X[,k]*X[,l])
-mean(X[,l]*X[,n])*mean(X[,i]*X[,j]*X[,k]*X[,m])
-mean(X[,l]*X[,m])*mean(X[,i]*X[,j]*X[,k]*X[,n])
-mean(X[,i]*X[,j]*X[,k])*mean(X[,l]*X[,m]*X[,n])
-mean(X[,j]*X[,k])*mean(X[,i]*X[,l]*X[,m]*X[,n])
+mean(X[,j]*X[,k])*mean(X[,m]*X[,n])*mean(X[,i]*X[,l])
+mean(X[,j]*X[,k])*mean(X[,l]*X[,n])*mean(X[,i]*X[,m])
+mean(X[,j]*X[,k])*mean(X[,l]*X[,m])*mean(X[,i]*X[,n])
-mean(X[,i]*X[,k])*mean(X[,j]*X[,l]*X[,m]*X[,n])
+mean(X[,i]*X[,k])*mean(X[,m]*X[,n])*mean(X[,j]*X[,l])
+mean(X[,i]*X[,k])*mean(X[,l]*X[,n])*mean(X[,j]*X[,m])
+mean(X[,i]*X[,k])*mean(X[,l]*X[,m])*mean(X[,j]*X[,n])
-mean(X[,i]*X[,j])*mean(X[,k]*X[,l]*X[,m]*X[,n])
+mean(X[,i]*X[,j])*mean(X[,m]*X[,n])*mean(X[,k]*X[,l])
+mean(X[,i]*X[,j])*mean(X[,l]*X[,n])*mean(X[,k]*X[,m])
+mean(X[,i]*X[,j])*mean(X[,l]*X[,m])*mean(X[,k]*X[,n]))
}
w33[upper.tri(w33, diag=TRUE)]<-w33ele
############# w23 ##################
w23index1 <- matrix(0,w33nrow*w22nrow,w33nrow)
for(i in 1:w33nrow){
w23index1[((w22nrow)*(i-1)+1):(w22nrow*i), i]<- rep(1, w22nrow)
}
w23index2 <- matrix(0,w33nrow*w22nrow,w22nrow)
for(i in 1:w33nrow){
w23index2[((w22nrow)*(i-1)+1):(w22nrow*i), (1:w22nrow)]<- diag(w22nrow)
}
w23index <- cbind(w23index1%*%w33index_2,
w23index2%*%w22index_2)
w23 <- matrix(0, w22nrow,
w33nrow)
for (o in 1:(w33nrow*w22nrow)){
i <- w23index[o,1]
j
<- w23index[o,2]
k
<- w23index[o,3]
l
<- w23index[o,4]
m
<- w23index[o,5]
w23[o]<-(mean(X[,i]*X[,j]*X[,k]*X[,l]*X[,m])
-mean(X[,i]*X[,j]*X[,k])*mean(X[,l]*X[,m])
-mean(X[,j]*X[,k])*mean(X[,i]*X[,l]*X[,m])
-mean(X[,i]*X[,k])*mean(X[,j]*X[,l]*X[,m])
-mean(X[,i]*X[,j])*mean(X[,k]*X[,l]*X[,m]))
}
################ combine #######################
w32 <- matrix(0, w33nrow, w22nrow)
w <- rbind(cbind(w22,w23),
cbind(w32,w33))
w <- w+t(w)-diag(diag(w))
iw<-try(solve(w))
if (class(iw)=='try-error')
{
iw<-ginv(w)
}
return(iw)
}
WM<-weight.matrix(M)
WD<-weight.matrix(D)
#Model specification and estimation
#C and E are non-normal
fr1<-function(p){
a<-p[1];
c<-p[2];
d<-p[3];
e<-p[4];
sig<-p[5];
sm2<-a^2+c^2+e^2+d^2;
sm12<-a^2+c^2+d^2;
sd2<-a^2+c^2+e^2+d^2;
sd12<-0.5*a^2+c^2+0.25*d^2;
sm3<-(c^3)*sig+(e^3)*sig;
sm21<-(c^3)*sig;
sd3<-(c^3)*sig+(e^3)*sig;
sd21<-(c^3)*sig;
M_calc_mom<-c(m1.2, m11, m2.2, m1.3, m21, m12, m2.3)
M_sig_mom<-c(sm2, sm12, sm2, sm3, sm21, sm21, sm3)
D_calc_mom<-c(d1.2, d11, d2.2, d1.3, d21, d12, d2.3)
D_sig_mom<-c(sd2, sd12, sd2, sd3, sd21, sd21, sd3)
#GLS estimation
GLS<-nrow(M)*t(M_calc_mom-M_sig_mom)%*%(WM)%*%(M_calc_mom-M_sig_mom)+nrow(D)*t(D_calc_mom-D_sig_mom)%*%(WD)%*%(D_calc_mom-D_sig_mom)
GLS<-GLS*1000000
return(GLS)
}
initial<-matrix(0, ncol=3,
nrow=initialrep)
for(i in
1:nrow(initial))
{
for(j in 1:ncol(initial))
{
initial[i,j]<-runif(1)
}
}
initialsig<-matrix(0, ncol=1, nrow=initialrep)
for(i in 1:nrow(initialsig))
{
for(j in 1:ncol(initialsig))
{
initialsig[i,j]<-rchisq(1,1)
}
}
fit1<-1000000000000000000
est1<-c(0,0,0,0,0)
se1<-c(0,0,0,0,0)
for(i in
1:nrow(initial))
{
initial.a<-abs(initial[i,1]);
initial.c<-abs(initial[i,2]);
initial.e<-sqrt(m1.2-m11);
initial.d<-abs(initial[i,3]);
initial.sig<-abs(initialsig[i,1]);
b0<-c(initial.a,
initial.c, initial.d, initial.e, initial.sig)
res1<-optim(b0,
fr1, method="BFGS", hessian=TRUE);
if (res1$convergence<1){
if (res1$value<fit1)
{fit1<-res1$value; est1<-res1$par; se1<-sqrt(diag(ginv(res1$hessian/1000000)));}}
}
#ACE model and C and E are non-normal
fr10<-function(p){
a<-p[1];
c<-p[2];
d<-0;
e<-p[3];
sig<-p[4];
sm2<-a^2+c^2+e^2+d^2;
sm12<-a^2+c^2+d^2;
sd2<-a^2+c^2+e^2+d^2;
sd12<-0.5*a^2+c^2+0.25*d^2;
sm3<-(c^3)*sig+(e^3)*sig;
sm21<-(c^3)*sig;
sd3<-(c^3)*sig+(e^3)*sig;
sd21<-(c^3)*sig;
M_calc_mom<-c(m1.2, m11, m2.2, m1.3, m21, m12, m2.3)
M_sig_mom<-c(sm2, sm12, sm2, sm3, sm21, sm21, sm3)
D_calc_mom<-c(d1.2, d11, d2.2, d1.3, d21, d12, d2.3)
D_sig_mom<-c(sd2, sd12, sd2, sd3, sd21, sd21, sd3)
#GLS estimation
GLS<-nrow(M)*t(M_calc_mom-M_sig_mom)%*%(WM)%*%(M_calc_mom-M_sig_mom)+nrow(D)*t(D_calc_mom-D_sig_mom)%*%(WD)%*%(D_calc_mom-D_sig_mom)
GLS<-GLS*1000000
return(GLS)
}
initial<-matrix(0, ncol=3,
nrow=initialrep)
for(i in
1:nrow(initial))
{
for(j in 1:ncol(initial))
{
initial[i,j]<-runif(1)
}
}
initialsig<-matrix(0, ncol=1, nrow=initialrep)
for(i in 1:nrow(initialsig))
{
for(j in 1:ncol(initialsig))
{
initialsig[i,j]<-rchisq(1,1)
}
}
fit10<-1000000000000000000
est10<-c(0,0,0,0)
se10<-c(0,0,0,0)
for(i in
1:nrow(initial))
{
initial.a<-abs(initial[i,1]);
initial.c<-abs(initial[i,2]);
initial.e<-sqrt(m1.2-m11);
initial.sig<-abs(initialsig[i,1]);
b0<-c(initial.a,
initial.c, initial.e, initial.sig)
res10<-optim(b0,
fr10, method="BFGS", hessian=TRUE);
if (res10$convergence<1){
if (res10$value<fit10)
{fit10<-res10$value; est10<-res10$par; se10<-sqrt(diag(ginv(res10$hessian/1000000)));}}
}
#D and E are non-normal
fr2<-function(p){
a<-p[1];
c<-p[2];
d<-p[3];
e<-p[4];
sig<-p[5];
signui<-p[6];
sm2<-a^2+c^2+e^2+d^2;
sm12<-a^2+c^2+d^2;
sd2<-a^2+c^2+e^2+d^2;
sd12<-0.5*a^2+c^2+0.25*d^2;
sm3<-(d^3+e^3)*sig;
sm21<-(d^3)*sig;
sd3<-(d^3+e^3)*sig;
sd21<-(d^3)*signui;
M_calc_mom<-c(m1.2, m11, m2.2, m1.3, m21, m12, m2.3)
M_sig_mom<-c(sm2, sm12, sm2, sm3, sm21, sm21, sm3)
D_calc_mom<-c(d1.2, d11, d2.2, d1.3, d21, d12, d2.3)
D_sig_mom<-c(sd2, sd12, sd2, sd3, sd21, sd21, sd3)
#GLS estimation
GLS<-nrow(M)*t(M_calc_mom-M_sig_mom)%*%(WM)%*%(M_calc_mom-M_sig_mom)+nrow(D)*t(D_calc_mom-D_sig_mom)%*%(WD)%*%(D_calc_mom-D_sig_mom)
GLS<-GLS*1000000
return(GLS)
}
initial<-matrix(0, ncol=3,
nrow=initialrep)
for(i in
1:nrow(initial))
{
for(j in 1:ncol(initial))
{
initial[i,j]<-runif(1)
}
}
initialsig<-matrix(0, ncol=2, nrow=initialrep)
for(i in 1:nrow(initialsig))
{
for(j in 1:ncol(initialsig))
{
initialsig[i,j]<-rchisq(1,1)
}
}
fit2<-1000000000000000000
est2<-c(0,0,0,0,0,0)
se2<-c(0,0,0,0,0,0)
for(i in
1:nrow(initial))
{
initial.a<-abs(initial[i,1]);
initial.c<-abs(initial[i,2]);
initial.e<-sqrt(m1.2-m11);
initial.d<-abs(initial[i,3]);
initial.sig<-abs(initialsig[i,1]);
initial.signui<-abs(initialsig[i,2]);
b0<-c(initial.a,
initial.c, initial.d, initial.e, initial.sig, initial.signui)
res2<-optim(b0,
fr2, method="BFGS", hessian=TRUE);
if (res2$convergence<1){
if (res2$value<fit2)
{fit2<-res2$value; est2<-res2$par; se2<-sqrt(diag(ginv(res2$hessian/1000000)));}}
}
#C, D, and E are
non-normal
fr3<-function(p){
a<-p[1];
c<-p[2];
d<-p[3];
e<-p[4];
sig<-p[5];
signui<-p[6];
sm2<-a^2+c^2+e^2+d^2;
sm12<-a^2+c^2+d^2;
sd2<-a^2+c^2+e^2+d^2;
sd12<-0.5*a^2+c^2+0.25*d^2;
sm3<-(d^3+c^3+e^3)*sig;
sm21<-(d^3+c^3)*sig;
sd3<-(d^3+c^3+e^3)*sig;
sd21<-(d^3)*signui+(c^3)*sig;
M_calc_mom<-c(m1.2, m11, m2.2, m1.3, m21, m12, m2.3)
M_sig_mom<-c(sm2, sm12, sm2, sm3, sm21, sm21, sm3)
D_calc_mom<-c(d1.2, d11, d2.2, d1.3, d21, d12, d2.3)
D_sig_mom<-c(sd2, sd12, sd2, sd3, sd21, sd21, sd3)
#GLS estimation
GLS<-nrow(M)*t(M_calc_mom-M_sig_mom)%*%(WM)%*%(M_calc_mom-M_sig_mom)+nrow(D)*t(D_calc_mom-D_sig_mom)%*%(WD)%*%(D_calc_mom-D_sig_mom)
GLS<-GLS*1000000
return(GLS)
}
initial<-matrix(0, ncol=3,
nrow=initialrep)
for(i in
1:nrow(initial))
{
for(j in 1:ncol(initial))
{
initial[i,j]<-runif(1)
}
}
initialsig<-matrix(0, ncol=2, nrow=initialrep)
for(i in 1:nrow(initialsig))
{
for(j in 1:ncol(initialsig))
{
initialsig[i,j]<-rchisq(1,1)
}
}
fit3<-1000000000000000000
est3<-c(0,0,0,0,0,0)
se3<-c(0,0,0,0,0,0)
for(i in
1:nrow(initial))
{
initial.a<-abs(initial[i,1]);
initial.c<-abs(initial[i,2]);
initial.e<-sqrt(m1.2-m11);
initial.d<-abs(initial[i,3]);
initial.sig<-abs(initialsig[i,1]);
initial.signui<-abs(initialsig[i,2]);
b0<-c(initial.a,
initial.c, initial.d, initial.e, initial.sig, initial.signui)
res3<-optim(b0,
fr3, method="BFGS", hessian=TRUE);
if (res3$convergence<1){
if (res3$value<fit3)
{fit3<-res3$value; est3<-res3$par; se3<-sqrt(diag(ginv(res3$hessian/1000000)));}}
}
RMSEA1<-sqrt(max((res1$value/1000000)/((14-5)*(nrow(M)+nrow(D)))-(1/((nrow(M)+nrow(D))-1)),0));
RMSEA2<-sqrt(max((res2$value/1000000)/((14-6)*(nrow(M)+nrow(D)))-(1/((nrow(M)+nrow(D))-1)),0));
RMSEA3<-sqrt(max((res3$value/1000000)/((14-6)*(nrow(M)+nrow(D)))-(1/((nrow(M)+nrow(D))-1)),0));
RMSEA10<-sqrt(max((res10$value/1000000)/((14-4)*(nrow(M)+nrow(D)))-(1/((nrow(M)+nrow(D))-1)),0));
CHI1<-res1$value/1000000;
CHI2<-res2$value/1000000;
CHI3<-res3$value/1000000;
CHI10<-res10$value/1000000;
AIC1<-res1$value/1000000-2*(14-5);
AIC2<-res2$value/1000000-2*(14-6);
AIC3<-res3$value/1000000-2*(14-6);
AIC10<-res10$value/1000000-2*(14-4);
BIC1<-res1$value/1000000-(log(nrow(M)+nrow(D)))*(14-5);
BIC2<-res2$value/1000000-(log(nrow(M)+nrow(D)))*(14-6);
BIC3<-res3$value/1000000-(log(nrow(M)+nrow(D)))*(14-6);
BIC10<-res10$value/1000000-(log(nrow(M)+nrow(D)))*(14-4);
est.results<-matrix(NA,8,13);
rownames(est.results)<-c("Model 1(C and E
are non-normal)", "SE of Model 1", "Model 2(D and E are
non-normal)", "SE of Model 2", "Model 3(C, D, and E are non-normal)", "SE of Model 3",
"Model 10(ACE model using 3rd)", "SE of Model 10");
colnames(est.results)<-c("a","c","d","e","RMSEA","AIC","BIC","Chisquared","df",
"skewness of C", "skewness of D", "skewness of
E", "E[D1^2*D2]");
est.results[1,1]<-est1[1]^2/(est1[1]^2+est1[2]^2+est1[3]^2+est1[4]^2);
est.results[1,2]<-est1[2]^2/(est1[1]^2+est1[2]^2+est1[3]^2+est1[4]^2);
est.results[1,3]<-est1[3]^2/(est1[1]^2+est1[2]^2+est1[3]^2+est1[4]^2);
est.results[1,4]<-est1[4]^2/(est1[1]^2+est1[2]^2+est1[3]^2+est1[4]^2);
est.results[1,5]<-RMSEA1;
est.results[1,6]<-AIC1;
est.results[1,7]<-BIC1;
est.results[1,8]<-CHI1;
est.results[1,9]<-9;
est.results[1,10]<-est1[5];
est.results[1,11]<-est1[5];
est.results[1,12]<-0;
est.results[1,13]<-NA;
est.results[2,1]<-sqrt((sqrt(est.results[1,1])/est1[1])^2*(se1[1]^2)*(2*est1[1])^2)
est.results[2,2]<-sqrt((sqrt(est.results[1,2])/est1[2])^2*(se1[2]^2)*(2*est1[2])^2)
est.results[2,3]<-sqrt((sqrt(est.results[1,3])/est1[3])^2*(se1[3]^2)*(2*est1[3])^2)
est.results[2,4]<-sqrt((sqrt(est.results[1,4])/est1[4])^2*(se1[4]^2)*(2*est1[4])^2)
est.results[2,10]<-se1[5]
est.results[2,11]<-se1[5]
est.results[3,1]<-est2[1]^2/(est2[1]^2+est2[2]^2+est2[3]^2+est2[4]^2);
est.results[3,2]<-est2[2]^2/(est2[1]^2+est2[2]^2+est2[3]^2+est2[4]^2);
est.results[3,3]<-est2[3]^2/(est2[1]^2+est2[2]^2+est2[3]^2+est2[4]^2);
est.results[3,4]<-est2[4]^2/(est2[1]^2+est2[2]^2+est2[3]^2+est2[4]^2);
est.results[3,5]<-RMSEA2;
est.results[3,6]<-AIC2;
est.results[3,7]<-BIC2;
est.results[3,8]<-CHI2;
est.results[3,9]<-8;
est.results[3,10]<-0;
est.results[3,11]<-est2[5];
est.results[3,12]<-est2[5];
est.results[3,13]<-est2[6];
est.results[4,1]<-sqrt((sqrt(est.results[3,1])/est2[1])^2*(se2[1]^2)*(2*est2[1])^2)
est.results[4,2]<-sqrt((sqrt(est.results[3,2])/est2[2])^2*(se2[2]^2)*(2*est2[2])^2)
est.results[4,3]<-sqrt((sqrt(est.results[3,3])/est2[3])^2*(se2[3]^2)*(2*est2[3])^2)
est.results[4,4]<-sqrt((sqrt(est.results[3,4])/est2[4])^2*(se2[4]^2)*(2*est2[4])^2)
est.results[4,11]<-se2[5]
est.results[4,12]<-se2[5]
est.results[4,13]<-se2[6]
est.results[5,1]<-est3[1]^2/(est3[1]^2+est3[2]^2+est3[3]^2+est3[4]^2);
est.results[5,2]<-est3[2]^2/(est3[1]^2+est3[2]^2+est3[3]^2+est3[4]^2);
est.results[5,3]<-est3[3]^2/(est3[1]^2+est3[2]^2+est3[3]^2+est3[4]^2);
est.results[5,4]<-est3[4]^2/(est3[1]^2+est3[2]^2+est3[3]^2+est3[4]^2);
est.results[5,5]<-RMSEA3;
est.results[5,6]<-AIC3;
est.results[5,7]<-BIC3;
est.results[5,8]<-CHI3;
est.results[5,9]<-8;
est.results[5,10]<-est3[5];
est.results[5,11]<-est3[5];
est.results[5,12]<-est3[5];
est.results[5,13]<-est3[6];
est.results[6,1]<-sqrt((sqrt(est.results[5,1])/est3[1])^2*(se3[1]^2)*(2*est3[1])^2)
est.results[6,2]<-sqrt((sqrt(est.results[5,2])/est3[2])^2*(se3[2]^2)*(2*est3[2])^2)
est.results[6,3]<-sqrt((sqrt(est.results[5,3])/est3[3])^2*(se3[3]^2)*(2*est3[3])^2)
est.results[6,4]<-sqrt((sqrt(est.results[5,4])/est3[4])^2*(se3[4]^2)*(2*est3[4])^2)
est.results[6,10]<-se3[5]
est.results[6,11]<-se3[5]
est.results[6,12]<-se3[5]
est.results[6,13]<-se3[6]
est.results[7,1]<-est10[1]^2/(est10[1]^2+est10[2]^2+est10[3]^2);
est.results[7,2]<-est10[2]^2/(est10[1]^2+est10[2]^2+est10[3]^2);
est.results[7,3]<-NA
est.results[7,4]<-est10[3]^2/(est10[1]^2+est10[2]^2+est10[3]^2);
est.results[7,5]<-RMSEA10;
est.results[7,6]<-AIC10;
est.results[7,7]<-BIC10;
est.results[7,8]<-CHI10;
est.results[7,9]<-10;
est.results[7,10]<-est10[4];
est.results[7,11]<-0;
est.results[7,12]<-est10[4];
est.results[7,13]<-NA;
est.results[8,1]<-sqrt((sqrt(est.results[7,1])/est10[1])^2*(se10[1]^2)*(2*est10[1])^2)
est.results[8,2]<-sqrt((sqrt(est.results[7,2])/est10[2])^2*(se10[2]^2)*(2*est10[2])^2)
est.results[8,3]<-NA
est.results[8,4]<-sqrt((sqrt(est.results[7,4])/est10[3])^2*(se10[3]^2)*(2*est10[3])^2)
est.results[8,10]<-se10[4]
est.results[8,12]<-se10[4]
sample.statistics.of.MZ.pairs<-matrix(NA,1,7);
rownames(sample.statistics.of.MZ.pairs)<-c("statictics");
colnames(sample.statistics.of.MZ.pairs)<-c("V(T1)","Cor(T1,T2)","V(T2)","Skew(T1)","Mean
of T1*T1*T2","Mean of T1*T2*T2","Skew(T2)");
sample.statistics.of.MZ.pairs[1,1]<-m1.2;
sample.statistics.of.MZ.pairs[1,2]<-m11;
sample.statistics.of.MZ.pairs[1,3]<-m2.2;
sample.statistics.of.MZ.pairs[1,4]<-m1.3;
sample.statistics.of.MZ.pairs[1,5]<-m21;
sample.statistics.of.MZ.pairs[1,6]<-m12;
sample.statistics.of.MZ.pairs[1,7]<-m2.3;
sample.statistics.of.DZ.pairs<-matrix(NA,1,7);
rownames(sample.statistics.of.DZ.pairs)<-c("statictics");
colnames(sample.statistics.of.DZ.pairs)<-c("V(T1)","Cor(T1,T2)","V(T2)","Skew(T1)","Mean
of T1*T1*T2","Mean of T1*T2*T2","Skew(T2)");
sample.statistics.of.DZ.pairs[1,1]<-d1.2;
sample.statistics.of.DZ.pairs[1,2]<-d11;
sample.statistics.of.DZ.pairs[1,3]<-d2.2;
sample.statistics.of.DZ.pairs[1,4]<-d1.3;
sample.statistics.of.DZ.pairs[1,5]<-d21;
sample.statistics.of.DZ.pairs[1,6]<-d12;
sample.statistics.of.DZ.pairs[1,7]<-d2.3;
round(est.results,
digits=3)
round(sample.statistics.of.MZ.pairs,
digits=3)
round(sample.statistics.of.DZ.pairs,
digits=3)
#Standardized variables are analyzed