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