Hi,
I need help with how to estimate the A, C,E for multiple variables and store them in a table format.
I have a example data with two values for each twin pair for 5 genes ( A,B,C,D,E).
I started out in following way from the Twin Univariate model script but it's not working.
Data<- read.table("ExampleData.txt",sep="\t",header=T) # Example Data attached.
Varsi<-c("GeneA_","GeneB_","GeneC_","GeneD_","GeneE_")
for(i in 1:5){
nv <- 1
selVars<- paste(Varsi[i],c(rep(1,nv),rep(2,nv)),sep="")
mzDatai<-subset(Data,Type=="MZ",selVars)
dzDatai<-subset(Data,Type=="DZ",selVars)
-- Data Summary------
summary(mzDatai)
colMeans(mzDatai,na.rm=TRUE)
cov(mzDatai,use="complete")
summary(dzDatai)
colMeans(dzDatai,na.rm=TRUE)
cov(dzDatai,use="complete")
---ACE Model----
univACEModel <- mxModel("univACE",
mxModel("ACE",
# Matrices a, c, and e to store a, c, and e path coefficients
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="a11", name="a" ),
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="c11", name="c" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="e11", name="e" ),
# Matrices A, C, and E compute variance components
mxAlgebra( expression=a %% t(a), name="A" ),
mxAlgebra( expression=c %% t(c), name="C" ),
mxAlgebra( expression=e %% t(e), name="E" ),
# Algebra to compute total variances and standard deviations (diagonal only)
mxAlgebra( expression=A+C+E, name="V" ),
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
mxAlgebra( expression=solve(sqrt(IV)), name="sd"),
# Matrix & Algebra for expected means vector
mxMatrix( type="Full", nrow=1, ncol=nv, name="M" ),
mxAlgebra( expression= cbind(Mean,Mean), name="expMean" ),
mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=20, label="mean", name="Mean" ),
mxAlgebra( expression= cbind(Mean,Mean), name="expMean" ),
# Algebra for expected variance/covariance matrix in MZ
mxAlgebra( expression= rbind ( cbind(A+C+E , A+C),
cbind(A+C , A+C+E)), name="expCovMZ" ),
# Algebra for expected variance/covariance matrix in DZ, note use of 0.5, converted to 1*1 matrix
mxAlgebra( expression= rbind ( cbind(A+C+E , 0.5%x%A+C),
cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )
),
mxModel("MZ",
mxData( observed=mzDatai, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovMZ", means="ACE.expMean", dimnames=selVars )
),
mxModel("DZ",
mxData( observed=dzDatai, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovDZ", means="ACE.expMean", dimnames=selVars)
),
mxAlgebra( expression=MZ.objective + DZ.objective, name="sumll" ),
mxAlgebraObjective("sumll"),
mxCI(c('ACE.A', 'ACE.C', 'ACE.E','ACE.V'))
)
}
univACEFit <- mxRun(univACEModel, intervals=T)
univACESumm <- summary(univACEFit)
univACESumm
LL_ACE <- mxEval(objective, univACEFit)
LL_ACE
Generate Table of Parameter Estimates using mxEval
pathEstimatesACE <- print(round(mxEval(cbind(ACE.a,ACE.c,ACE.e), univACEFit),4))
varComponentsACE <- print(round(mxEval(cbind(ACE.A/ACE.V,ACE.C/ACE.V,ACE.E/ACE.V), univACEFit),4))
rownames(pathEstimatesACE) <- 'pathEstimates'
colnames(pathEstimatesACE) <- c('a','c','e')
rownames(varComponentsACE) <- 'varComponentsACE'
colnames(varComponentsACE) <- c('a^2','c^2','e^2')
pathEstimatesACE
varComponentsACE
Any help would be of great help for highthroughput data with ACE estimation for large no.of variables( genes).