# ----- NOT R code MZ twin data English1 Math1 SocSci1 NatSci1 Vocab1 English2 Math2 SocSci2 NatSci2 Vocab2 22.5236 16.0187 14.8236 15.5180 15.6883 17.0398 14.5703 13.6112 14.8395 14.1761 16.0187 39.0364 19.0094 23.9329 16.8449 14.5703 29.0442 17.5001 21.9457 15.6345 14.8236 19.0094 23.8313 18.6608 18.4180 13.6112 17.5001 18.3289 17.6795 17.0395 15.5180 23.9329 18.6608 33.5830 16.1748 14.8395 21.9457 17.6795 23.0796 14.9928 15.6883 16.8449 18.4180 16.1748 24.1408 14.1761 15.6345 17.0395 14.9928 20.7278 17.0398 14.5703 13.6112 14.8395 14.1761 22.5236 16.0187 14.8236 15.5180 15.6883 14.5703 29.0442 17.5001 21.9457 15.6345 16.0187 39.0364 19.0094 23.9329 16.8449 13.6112 17.5001 18.3289 17.6795 17.0395 14.8236 19.0094 23.8313 18.6608 18.4180 14.8395 21.9457 17.6795 23.0796 14.9928 15.5180 23.9329 18.6608 33.5830 16.1748 14.1761 15.6345 17.0395 14.9928 20.7278 15.6883 16.8449 18.4180 16.1748 24.1408 DZ twin data English1 Math1 SocSci1 NatSci1 Vocab1 English2 Math2 SocSci2 NatSci2 Vocab2 20.4282 15.7278 13.6828 15.1247 13.7292 11.5891 9.7854 9.0182 10.5496 9.8035 15.7278 40.0221 18.1005 23.3223 16.4543 9.7854 19.3924 12.0687 15.5302 11.1979 13.6828 18.1005 23.7431 18.7302 17.2567 9.0182 12.0687 12.4381 13.8789 11.9365 15.1247 23.3223 18.7302 31.8816 16.8254 10.5496 15.5302 13.8789 17.6237 12.2275 13.7292 16.4543 17.2567 16.8254 20.9661 9.8035 11.1979 11.9365 12.2275 13.4137 11.5891 9.7854 9.0182 10.5496 9.8035 20.4282 15.7278 13.6828 15.1247 13.7292 9.7854 19.3924 12.0687 15.5302 11.1979 15.7278 40.0221 18.1005 23.3223 16.4543 9.0182 12.0687 12.4381 13.8789 11.9365 13.6828 18.1005 23.7431 18.7302 17.2567 10.5496 15.5302 13.8789 17.6237 12.2275 15.1247 23.3223 18.7302 31.8816 16.8254 9.8035 11.1979 11.9365 12.2275 13.4137 13.7292 16.4543 17.2567 16.8254 20.9661 # --- R Code starts here # --- mxModel: NMTwins # --- mxModel: National Merit Twins: General Model require(OpenMx) thisOMxModel <- mxModel(model='thisOMxModel', name='thisOMxModel') # --- Macro Variables # nPheno: Number of phenotypes nPheno <- 5 VAstart <- 4 sqrtStdEstart <- 3 # --- matrices # VA: Additive-genetic covariance # StdE: Diagonal (standard deviations of the environmental variables) # Re: within-inidividual correlatiobns for environment variables # phi: cross-sib correlations for environment variables thisOMxModel <- mxModel(thisOMxModel,mxMatrix(name='VA', type='Symm', free=TRUE, nrow=nPheno, ncol=nPheno)) thisOMxModel <- mxModel(thisOMxModel,mxMatrix(name='sqrtStdE', type='Diag', free=TRUE, nrow=nPheno, ncol=nPheno)) thisOMxModel <- mxModel(thisOMxModel,mxMatrix(name='Re', type='Stand', free=TRUE, nrow=nPheno, ncol=nPheno)) thisOMxModel <- mxModel(thisOMxModel,mxMatrix(name='phi', type='Symm', free=TRUE, nrow=nPheno, ncol=nPheno)) # --- algebra thisOMxModel <- mxModel(thisOMxModel, mxAlgebra(sqrtStdE * sqrtStdE, name='StdE')) thisOMxModel <- mxModel(thisOMxModel, mxAlgebra(StdE %*% Re %*% StdE, name='VE')) thisOMxModel <- mxModel(thisOMxModel, mxAlgebra(StdE %*% phi %*% StdE, name='S')) thisOMxModel <- mxModel(thisOMxModel, mxAlgebra(VA + VE, name='VP')) thisOMxModel <- mxModel(thisOMxModel, mxAlgebra(VA + S, name='R12MZ')) thisOMxModel <- mxModel(thisOMxModel, mxAlgebra(.5 %x% VA + S, name='R12DZ')) thisOMxModel <- mxModel(thisOMxModel, mxAlgebra(cbind(rbind(VP, R12MZ), rbind(R12MZ,VP)), name='PreCovMZ', dimnames=dimnames(ObsCovMZ))) thisOMxModel <- mxModel(thisOMxModel, mxAlgebra(cbind(rbind(VP, R12DZ), rbind(R12DZ,VP)), name='PreCovDZ', dimnames=dimnames(ObsCovDZ))) # --- data objects # MZ: MZ Twins # DZ: DZ Twins thisOMxModel <- mxModel(thisOMxModel, mxModel(name='MZ', mxData(observed=ObsCovMZ, type='cov', numObs=509), mxMLObjective(covariance='thisOMxModel.PreCovMZ'))) thisOMxModel <- mxModel(thisOMxModel, mxModel(name='DZ', mxData(observed=ObsCovDZ, type='cov', numObs=330), mxMLObjective(covariance='thisOMxModel.PreCovDZ'))) # start values diag(thisOMxModel@matrices$VA@values) <- VAstart diag(thisOMxModel@matrices$sqrtStdE@values) <- sqrtStdEstart # this does not work #test <- mxRun(thisOMxModel) #summary(test) # does hermine's code fix it? thisOMxModel2 <- mxModel(thisOMxModel, mxAlgebra(MZ.objective + DZ.objective, name='fval'), mxAlgebraObjective('fval')) test2 <- mxRun(thisOMxModel2) test2 <- mxRun(test2) summary(test2) # which formula for -2Log(L)? fmz1 <- 509*(log(det(as.matrix(test2@algebras$PreCovMZ@result))) + sum(diag(as.matrix(ObsCovMZ) %*% solve(as.matrix(test2@algebras$PreCovMZ@result))))) fmz1 fmz2 <- 509*(log(det(as.matrix(test2@algebras$PreCovMZ@result)))) + 508*(sum(diag(as.matrix(ObsCovMZ) %*% solve(as.matrix(test2@algebras$PreCovMZ@result))))) fmz2 fmz3 <- 508*(log(det(as.matrix(test2@algebras$PreCovMZ@result)))) + 508*(sum(diag(as.matrix(ObsCovMZ) %*% solve(as.matrix(test2@algebras$PreCovMZ@result))))) fmz3 fdz1 <- 330*(log(det(as.matrix(test2@algebras$PreCovDZ@result))) + sum(diag(as.matrix(ObsCovDZ) %*% solve(as.matrix(test2@algebras$PreCovDZ@result))))) fdz1 fdz2 <- 330*(log(det(as.matrix(test2@algebras$PreCovDZ@result)))) + 329*(sum(diag(as.matrix(ObsCovDZ) %*% solve(as.matrix(test2@algebras$PreCovDZ@result))))) fdz2 fdz3 <- 329*(log(det(as.matrix(test2@algebras$PreCovDZ@result)))) + 329*(sum(diag(as.matrix(ObsCovDZ) %*% solve(as.matrix(test2@algebras$PreCovDZ@result))))) fdz3 diff1 <- test2@output$Minus2LogLikelihood - (fmz1 + fdz1) diff2 <- test2@output$Minus2LogLikelihood - (fmz2 + fdz2) diff3 <- test2@output$Minus2LogLikelihood - (fmz3 + fdz3) diff1 diff2 diff3 thisTrace <- sum(diag(as.matrix(ObsCovMZ) %*% solve(as.matrix(test2@algebras$PreCovMZ@result)))) thisFDet <- log(det(as.matrix(test2@algebras$PreCovMZ@result))) - log(det(as.matrix(ObsCovMZ))) chimz <- 508*(thisFDet + thisTrace - 10) thisTrace <- sum(diag(as.matrix(ObsCovDZ) %*% solve(as.matrix(test2@algebras$PreCovDZ@result)))) thisFDet <- log(det(as.matrix(test2@algebras$PreCovDZ@result))) - log(det(as.matrix(ObsCovDZ))) chidz <- 329*(thisFDet + thisTrace - 10) chicheck <- (509/508)*chimz + (330/329)*chidz # chicheck should be close to zero zero <- chicheck - 21.483 zero