#------------------------------------------------------------------------------ # Author: Michael D. Hunter # Date: 2020-03-03 08:11:10 # Filename: threeLevelCholesky.R # Purpose: Create an encapsulated minumal working example of 2- and 3-level # twin models. The later is a twin-neighborhood model. These are based on # fake (simulated) data. #------------------------------------------------------------------------------ # setwd('./Projects/OpenMx/Boulder2020/Multilevel') # setwd('/home/mikeh/2020/Multilevel') #------------------------------------------------------------------------------ # Load packages require(OpenMx) #------------------------------------------------------------------------------ # Load data load('ThreeLevelSmall.RData') # Inspect and comment childData <- gMod1$data$observed familyData <- gMod1$family$data$observed countyData <- gMod1$family$neighborhood$data$observed dim(childData) head(childData) # The 'rel' variable is the coefficient of relatedness. # The 'relu' variable is its complement: 1-rel # Because the definition variable in these models works on the variances # rather than the loadings, rel and relu are used instead of their # square roots. dim(familyData) head(familyData) # Note the inclusion of 'rel' in these data. # This coefficient operates on the family level, not the # child level. dim(countyData) head(countyData) #------------------------------------------------------------------------------ # Set working/structure variables # CP = Conduct Problems # ES = Externalizing # HYP = HYPeractivity # DE = Defiant thePheno <- c('CPFIN')#c('CPFIN', 'ESFIN', 'HYPFIN', 'DEFIN') numPheno <- length(thePheno) mLab <- paste('mean', thePheno, sep='_') bCLatent <- paste('C', 1:numPheno, sep='') bALatent <- paste('AC', 1:numPheno, sep='') ALatent <- paste('A', 1:numPheno, sep='') bLatent <- c(bCLatent, bALatent) wELatent <- paste('E', 1:numPheno, sep='') wALatent <- paste('AU', 1:numPheno, sep='') wLatent <- c(wELatent, wALatent) bvLab <- c(rep(NA, numPheno), rep("data.rel", numPheno)) #between variance labels wvLab <- c(rep(NA, numPheno), rep("data.relu", numPheno)) #within variance labels # A function that take a from and a to, then makes the # Cholesky paths between them. # Optionally, fromLab and toLab can be specified. # These are used to keep the "between" A and "within" A # components equal. # Further named arugments to mxPath can be given in ... cholHelper <- function(from, to, fromLab=from, toLab=to, ...){ lab <- outer(fromLab, toLab, paste, sep='_') lab <- lab[upper.tri(lab, diag=TRUE)] outFrom <- matrix(from, nrow=length(from), ncol=length(to)) outFrom <- outFrom[upper.tri(outFrom, diag=TRUE)] outTo <- matrix(to, nrow=length(from), ncol=length(to), byrow=TRUE) outTo <- outTo[upper.tri(outTo, diag=TRUE)] lbounds <- matrix(NA, nrow=length(from), ncol=length(to)) diag(lbounds) <- 1e-6 lbounds <- lbounds[upper.tri(lbounds, diag=TRUE)] return(mxPath(from=outFrom, to=outTo, labels=lab, lbound=lbounds, ...)) } #------------------------------------------------------------------------------ #------------------------------------------------------------------------------ #------------------------------------------------------------------------------ # Create two-level model #------------------------------------------------------------------------------ #------------------------------------------------------------------------------ #------------------------------------------------------------------------------ #------------------------------------------------------------------------------ # Between Model #Version 2 that is slicker and quicker bModel2 <- mxModel('family', type="RAM", mxData(type="raw", observed=familyData, primaryKey="Family_ID"), latentVars = bLatent, mxPath(bLatent, arrows=2, values=1, free=FALSE, labels=bvLab)) #------------------------------------------------------------------------------ # Within Model #Version 2 that is slicker and quicker wModel2 <- mxModel('child', type="RAM", bModel2, mxData(type="raw", observed=childData, sort=FALSE), manifestVars = thePheno, latentVars = wLatent, mxPath(from="one", to=thePheno, arrows=1, free=TRUE, values=2, labels=mLab), mxPath(wLatent, arrows=2, values=1, free=FALSE, labels=wvLab), cholHelper(wALatent, thePheno, fromLab=ALatent, values=.8), cholHelper(wELatent, thePheno, values=.8), cholHelper(paste0('family.', bCLatent), thePheno, fromLab=bCLatent, joinKey='Family_ID', values=.8), cholHelper(paste0('family.', bALatent), thePheno, fromLab=ALatent, joinKey='Family_ID', values=.8)) #------------------------------------------------------------------------------ # Run 'em wRun2 <- mxRun(wModel2) #------------------------------------------------------------------------------ # Take a look summary(wRun2) p <- coef(wRun2) AChol <- vech2full(p[grep('^A', names(p))]) AChol[upper.tri(AChol)] <- 0 CChol <- vech2full(p[grep('^C', names(p))]) CChol[upper.tri(CChol)] <- 0 EChol <- vech2full(p[grep('^E', names(p))]) EChol[upper.tri(EChol)] <- 0 ACov <- AChol %*% t(AChol) CCov <- CChol %*% t(CChol) ECov <- EChol %*% t(EChol) (ACE <- matrix(c(diag(ACov), diag(CCov), diag(ECov)), ncol=3, dimnames=list(thePheno, c('A', 'C', 'E')))) # A C E #CPFIN 2.808061 0.9066664 1.113353 #ESFIN 1.364529 0.2430282 1.594080 #HYPFIN 2.818496 0.2222725 3.478574 #DEFIN 1.949715 0.2683077 1.109120 # Above are results from real data (ACES <- ACE/rowSums(ACE)) # A C E #CPFIN 0.5816103 0.18779024 0.2305994 #ESFIN 0.4261973 0.07590745 0.4978953 #HYPFIN 0.4323282 0.03409431 0.5335774 #DEFIN 0.5860028 0.08064205 0.3333552 # Above are results from real data rowSums(ACE) sapply(childData[,thePheno, drop=FALSE], var, na.rm=TRUE) # 'Genetic correlations' for the multivariate case (not covered) cov2cor(ACov) cov2cor(CCov) cov2cor(ECov) #------------------------------------------------------------------------------ #------------------------------------------------------------------------------ #------------------------------------------------------------------------------ # Create three-level model #------------------------------------------------------------------------------ #------------------------------------------------------------------------------ #------------------------------------------------------------------------------ #------------------------------------------------------------------------------ # Neighborhood / Site Model cNLatent <- paste('CN', 1:numPheno, sep='') nModel3 <- mxModel('neighborhood', type="RAM", mxData(type="raw", observed=countyData, primaryKey="County_ID"), latentVars = cNLatent, mxPath(cNLatent, arrows=2, values=1, free=FALSE, lbound=1e-6) ) #------------------------------------------------------------------------------ # Between Model bNLatent <- paste('N', 1:numPheno, sep='') fLatent <- c(bLatent, bNLatent) fvLab <- c(bvLab, rep(NA, numPheno)) fModel3 <- mxModel('family', type="RAM", nModel3, mxData(type="raw", observed=familyData, primaryKey="Family_ID"), latentVars = fLatent, mxPath(bLatent, arrows=2, values=1, free=FALSE, labels=bvLab), mxPath(paste0('neighborhood.', cNLatent), bNLatent, values=1, free=FALSE, joinKey='County_ID') ) #------------------------------------------------------------------------------ # Within Model #Version 2 that is slicker and quicker cModel3 <- mxModel('child', type="RAM", fModel3, mxData(type="raw", observed=childData, sort=FALSE), manifestVars = thePheno, latentVars = wLatent, mxPath(from="one", to=thePheno, arrows=1, free=TRUE, values=2, labels=mLab), mxPath(wLatent, arrows=2, values=1, free=FALSE, labels=wvLab), cholHelper(wALatent, thePheno, fromLab=ALatent, values=.8, ubound=10), cholHelper(wELatent, thePheno, values=.8, ubound=10), cholHelper(paste0('family.', bCLatent), thePheno, fromLab=bCLatent, joinKey='Family_ID', values=.8, ubound=10), cholHelper(paste0('family.', bALatent), thePheno, fromLab=ALatent, joinKey='Family_ID', values=.8, ubound=10), cholHelper(paste0('family.', bNLatent), thePheno, fromLab=bNLatent, joinKey='Family_ID', values=.8, ubound=10)) #------------------------------------------------------------------------------ # Run 'em cRun3 <- mxRun(cModel3) #------------------------------------------------------------------------------ # Take a look summary(cRun3) p <- coef(cRun3) AChol <- vech2full(p[grep('^A', names(p))]) AChol[upper.tri(AChol)] <- 0 CChol <- vech2full(p[grep('^C', names(p))]) CChol[upper.tri(CChol)] <- 0 EChol <- vech2full(p[grep('^E', names(p))]) EChol[upper.tri(EChol)] <- 0 NChol <- vech2full(p[grep('^N', names(p))]) NChol[upper.tri(NChol)] <- 0 ACov <- AChol %*% t(AChol) CCov <- CChol %*% t(CChol) ECov <- EChol %*% t(EChol) NCov <- NChol %*% t(NChol) # Variances for each of A, C, E, and N (ACEN <- matrix(c(diag(ACov), diag(CCov), diag(ECov), diag(NCov)), ncol=4, dimnames=list(thePheno, c('A', 'C', 'E', 'N')))) # Proportion of variances for A, C, E, and N (ACENS <- ACEN/rowSums(ACEN)) # Total estimated variance compared to child-level total variance rowSums(ACEN) sapply(childData[,thePheno, drop=FALSE], var, na.rm=TRUE) # 'Genetic correlations' cov2cor(ACov) cov2cor(CCov) cov2cor(ECov) cov2cor(NCov) #------------------------------------------------------------------------------ # Compare # Compare the proportions of variance attributed to A, C, and E in the two level model # to the proportions of variance attributed to A, C, E, and N in the three level model. # What happend to C? # Can you think of a situation where the higher level variance (say is site, or # rater, or data collector, or fMRI scanner, or SNP chip manufacturer) # impacts E instead? What about A? #------------------------------------------------------------------------------ # Optional # Try fitting the same model to Hyperactivity instead. # Does it have the same C -> N phenomena as conduct problems? #------------------------------------------------------------------------------ # Done