rm(list=ls(all=T)) mzmData <- read.table("mzm.dat",sep="",na.strings=".", header=T) dzmData <- read.table("dzm.dat",sep="",na.strings=".", header=T) mzfData <- read.table("mzf.dat",sep="",na.strings=".", header=T) dzfData <- read.table("dzf.dat",sep="",na.strings=".", header=T) dzoData <- read.table("dzo.dat",sep="",na.strings=".", header=T) #---------------------------------------------------------------------------------------------- # Start of the Model in which a SNP is included as definition variable affecting the common factor #---------------------------------------------------------------------------------------------- require(OpenMx) manifestVars <- c("mo3t1","fa3t1","mo7t1","fa7t1","mo10t1","fa10t1","mo12t1","fa12t1","mo3t2","fa3t2","mo7t2","fa7t2","mo10t2","fa10t2","mo12t2","fa12t2") # Means, Freely estimated, starting values derived from estimation of the model in the entire (both genotyped and non genotyped) dataset AllMean <- mxMatrix( type="Full", nrow=1, ncol=8, free=T, values=0, labels=c("meanm3","meanf3","meanm7","meanf7","meanm10","meanf10","meanm12","meanf12"), name="Mean") # Regressions of individual items on sex , freely estimated, starting values based on entire dataset AllBeta <- mxMatrix( type="Full", nrow=1, ncol=8, free=T, values=c(-0.1,-0.1,rep(-0.7,6)), labels=c("beta1","beta2","beta3","beta4","beta5","beta6","beta7","beta8"), name="Beta") # Regression of the common factor on the SNP value, this is the main parameter of interest! ############################## ASSIGNMENT ####################################### # Below (in ALLSnpBeta) we need to adjust the number of columns to accomodate 5 SNPs instead of one # It is now marked to a questionmark, it should be easy to figure out the required value. ############################## ASSIGNMENT ####################################### AllSnpBeta <- mxMatrix( type="Full", nrow=1, ncol=?, free=T, values= rep(.01, 5), labels=c("SNP1", "SNP2", "SNP3", "SNP4", "SNP5"), name="BetaSNP") # Factor loadings of the items on the common factor, fixed, based on the entire dataset AllLoadings <- mxMatrix( type="Full", nrow=8, ncol=1, free=F, values=c(1.2337,1.2310,2.3359,2.0936,2.5046,2.2797,2.3230,2.1039), name="facL") # Correlation between the common factors of twin 1 and twin 2, freely estimated (different accross zygosities) AllMZCor <- mxMatrix( type="Full", nrow=1, ncol=1, free=T, values=.6, labels="MZcor", name="facCorMZ") AllDZCor <- mxMatrix( type="Full", nrow=1, ncol=1, free=T, values=.1, labels="DZcor", name="facCorDZ") # Factor loadings of the items on the (item-specific) residual factors, fixed at values derived from the entire sample AllResidual <- mxMatrix( type="Diag", nrow=8, ncol=8, free=F, values=c(1.7753, 1.7186, 1.8036, 1.7091, 1.7403, 1.6928, 1.7810, 1.7563), name="res") # (Diagonal matrix of) covariances between the residual factors of twin 1 - twin 2, fixed at values derived from the entire sample AllResidualMZCov <- mxMatrix( type="Diag", nrow=8, ncol=8, free=F, values=c(0.6465, 0.6334, 0.5994, 0.6365, 0.5652, 0.6110, 0.6231, 0.6472), name="resCorMZ") AllResidualDZCov <- mxMatrix( type="Diag", nrow=8, ncol=8, free=F, values=c(0.1539, 0.1807, 0.3259, 0.3872, 0.3108, 0.3767, 0.3063, 0.4101), name="resCorDZ") # Matrix Algebra # Variance/covariance between the items (within twin 1 [twin 2]) due to the common factor AllFactorVariance <- mxAlgebra( expression= facL%*%t(facL), name="facVar") # Twin 1 - twin 2 covariance due to the common factor AllFactorCovMZ <- mxAlgebra( expression= facL%*%facCorMZ%*%t(facL), name="facCovMZ") AllFactorCovDZ <- mxAlgebra( expression= facL%*%facCorDZ%*%t(facL), name="facCovDZ") # Residual item variance within twin 1 (twin 2) AllResVar <- mxAlgebra( expression= res%*%t(res), name="resVar") # Residual covariance between the items of twin 1 - twin 2 AllResCovMZ <- mxAlgebra( expression= res%*%resCorMZ%*%t(res), name="resCovMZ") AllResCovDZ <- mxAlgebra( expression= res%*%resCorDZ%*%t(res), name="resCovDZ") # The expected 16 x 16 phenotypic covariance matrix AllExpCovMZ <- mxAlgebra( expression= rbind(cbind(facVar+resVar, facCovMZ+resCovMZ), cbind(facCovMZ+resCovMZ, facVar+resVar)), name="covMZ") AllExpCovDZ <-mxAlgebra( expression= rbind(cbind(facVar+resVar, facCovDZ+resCovDZ), cbind(facCovDZ+resCovDZ, facVar+resVar)), name="covDZ") AllModel <- mxModel("All",AllMean,AllBeta,AllSnpBeta,AllLoadings,AllMZCor,AllDZCor,AllResidual, AllResidualMZCov,AllResidualDZCov,AllFactorVariance,AllFactorCovMZ,AllFactorCovDZ, AllResVar,AllResCovMZ,AllResCovDZ,AllExpCovMZ,AllExpCovDZ) #---------------------------------------------------------------------------------------------- MZMData <- mxData(mzmData, type="raw") MZMSex <- mxMatrix( type="Full", nrow=1, ncol=2, free=F, labels=c("data.SEXt1","data.SEXt2"), name="sex") ######################################### ASSIGNMENT ######################################### # We are now reading in 5 SNPs at the same time! Adapt the dimentions of the matrix to accomodate # 5 SNPs for 2 twins. If you nee a hint look back at the 1SNP script! adapt from there! # also you need to add labels for all SNPs i included the labels statement below add it to the # MZMSnp object directly below! #labels=c("data.rs6265t1", "data.rs4680t1", "data.rs6314t1", "data.rs1007023t1", # "data.rs12231356t1", "data.rs6265t2", "data.rs4680t2", # "data.rs6314t2", "data.rs1007023t2", "data.rs12231356t2") # ######################################### ASSIGNMENT ######################################### MZMSnp <- mxMatrix( type="Full", nrow=?, ncol=?, free=F, , name="snp") MZMMean <- mxAlgebra( expression= cbind(All.Mean,All.Mean) + sex%x%All.Beta + (All.BetaSNP%*%snp)%x%t(All.facL), name="mean" ) MZMObj <-mxFIMLObjective( covariance="All.covMZ", means="mean", dimnames=manifestVars) MZMModel <- mxModel("MZM",MZMData,MZMSex,MZMSnp,MZMMean,MZMObj) #---------------------------------------------------------------------------------------------- DZMData <- mxData(dzmData, type="raw") DZMSex <- mxMatrix( type="Full", nrow=1, ncol=2, free=F, labels=c("data.SEXt1","data.SEXt2"), name="sex") DZMSnp <- mxMatrix( type="Full", nrow=5, ncol=2, free=F, labels=c("data.rs6265t1", "data.rs4680t1", "data.rs6314t1", "data.rs1007023t1", "data.rs12231356t1", "data.rs6265t2", "data.rs4680t2", "data.rs6314t2", "data.rs1007023t2", "data.rs12231356t2"), name="snp") DZMMean <- mxAlgebra( expression= cbind(All.Mean,All.Mean) + sex%x%All.Beta + (All.BetaSNP%*%snp)%x%t(All.facL), name="mean" ) DZMObj <- mxFIMLObjective( covariance="All.covDZ", means="mean", dimnames=manifestVars) DZMModel <- mxModel("DZM",DZMData,DZMSex,DZMSnp,DZMMean,DZMObj ) #---------------------------------------------------------------------------------------------- MZFData <- mxData(mzfData, type="raw") MZFSex <- mxMatrix( type="Full", nrow=1, ncol=2, free=F, labels=c("data.SEXt1","data.SEXt2"), name="sex") MZFSnp <- mxMatrix( type="Full", nrow=5, ncol=2, free=F, labels=c("data.rs6265t1", "data.rs4680t1", "data.rs6314t1", "data.rs1007023t1", "data.rs12231356t1", "data.rs6265t2", "data.rs4680t2", "data.rs6314t2", "data.rs1007023t2", "data.rs12231356t2"), name="snp") MZFMean <- mxAlgebra( expression= cbind(All.Mean,All.Mean) + sex%x%All.Beta + (All.BetaSNP%*%snp)%x%t(All.facL), name="mean" ) MZFObj <- mxFIMLObjective( covariance="All.covMZ", means="mean", dimnames=manifestVars) MZFModel <- mxModel("MZF",MZFData,MZFSex,MZFSnp,MZFMean,MZFObj) #---------------------------------------------------------------------------------------------- DZFData <- mxData(dzfData, type="raw") DZFSex <- mxMatrix( type="Full", nrow=1, ncol=2, free=F, labels=c("data.SEXt1","data.SEXt2"), name="sex") DZFSnp <- mxMatrix( type="Full", nrow=5, ncol=2, free=F, labels=c("data.rs6265t1", "data.rs4680t1", "data.rs6314t1", "data.rs1007023t1", "data.rs12231356t1", "data.rs6265t2", "data.rs4680t2", "data.rs6314t2", "data.rs1007023t2", "data.rs12231356t2"), name="snp") DZFMean <- mxAlgebra( expression= cbind(All.Mean,All.Mean) + sex%x%All.Beta + (All.BetaSNP%*%snp)%x%t(All.facL), name="mean" ) DZFObj <- mxFIMLObjective( covariance="All.covDZ", means="mean", dimnames=manifestVars) DZFModel <- mxModel("DZF",DZFData,DZFSex,DZFSnp,DZFMean,DZFObj) #---------------------------------------------------------------------------------------------- DOSData <- mxData(dzoData, type="raw") DOSSex <- mxMatrix( type="Full", nrow=1, ncol=2, free=F, labels=c("data.SEXt1","data.SEXt2"), name="sex") DOSSnp <- mxMatrix( type="Full", nrow=5, ncol=2, free=F, labels=c("data.rs6265t1", "data.rs4680t1", "data.rs6314t1", "data.rs1007023t1", "data.rs12231356t1", "data.rs6265t2", "data.rs4680t2", "data.rs6314t2", "data.rs1007023t2", "data.rs12231356t2"), name="snp") DOSMean <- mxAlgebra( expression= cbind(All.Mean,All.Mean) + sex%x%All.Beta + (All.BetaSNP%*%snp)%x%t(All.facL), name="mean" ) DOSObj <- mxFIMLObjective( covariance="All.covDZ", means="mean", dimnames=manifestVars) DOSModel <- mxModel("DOS",DOSData,DOSSex,DOSSnp, DOSSnp,DOSMean,DOSObj) #---------------------------------------------------------------------------------------------- SUMObj <- mxAlgebra( MZM.objective + DZM.objective + MZF.objective + DZF.objective + DOS.objective, name="h1") SUMAlg <- mxAlgebraObjective("h1") factor5SNPsModel <- mxModel("factorSNPs",AllModel,MZMModel,DZMModel,MZFModel,DZFModel,DOSModel,SUMObj,SUMAlg) factor5SNPsFit <- mxRun(factor5SNPsModel) factor5SNPsSumm <- summary(factor5SNPsFit) factor5SNPsSumm #---------------------------------------------------------------------------------------------- # Drop the SNP effect: factorNoSNPsModel <- factor5SNPsModel ############ Assignment, create the new values for the BetaSNP matrices ############ factorNoSNPsModel$All@matrices$BetaSNP@free[] = ?????? factorNoSNPsModel$All@matrices$BetaSNP@values[] = ? factorNoSNPsFit = mxRun(factorNoSNPsModel) factorNoSNPsSumm <- summary(factorNoSNPsFit) factorNoSNPsSumm mxCompare(factor5SNPsFit,factorNoSNPsFit)