Nothing
#' AE model using nuclear family trios
#'
#' @param model a linear mixed model formula, see example below.
#' @param random random effect, see exampe below.
#' @param data data to be analyzed.
#' @param seed random number seed.
#' @param n.sim number of simulations.
#' @param verbose a flag for printing out results.
#'
#' @details
#' This function is adapted from example 7.1 of Rabe-Hesketh et al. (2008). It
#' also procides heritability estimate and confidence intervals.
#'
#' @export
#' @return
#' The returned value is a list containing:
#' - lme.result the linear mixed model result.
#' - h2 the heritability estimate.
#' - CL confidence intervals.
#'
#' @references
#' \insertRef{rh08}{gap}
#'
#' @examples
#' \dontrun{
#' require(gap.datasets)
#' AE3(bwt ~ male + first + midage + highage + birthyr,
#' list(familyid = pdIdent(~var1 + var2 + var3 -1)), mfblong)
#' }
#'
#' @author Jing Hua Zhao
#' @note Adapted from f.mbf.R from the paper.
#' @keywords models
#' @keywords htest
AE3 <- function(model, random, data, seed=1234, n.sim=50000, verbose=TRUE)
{
for(p in c("MASS", "nlme")) {
if (length(grep(paste("^package:", p, "$", sep=""), search())) == 0) {
if (!requireNamespace(p, quietly = TRUE))
warning(paste("AE3 needs package `", p, "' to be fully functional; please install", sep=""))
}
}
res <- nlme::lme(model, random = random, data = data, method = "ML")
lns2 <- attr(res$apVar, "Pars")[1]
lns1 <- attr(res$apVar, "Pars")[2]
var22 <- res$apVar[1, 1]
cov21 <- res$apVar[2, 1]
var11 <- res$apVar[2, 2]
h2 <- exp(2*lns2)/(exp(2*lns1) + exp(2*lns2))
deriv1 <- -2*exp(2*lns1)*exp(2*lns2)/(exp(2*lns1) + exp(2*lns2))^2
deriv2 <- -2*exp(4*lns2)/(exp(2*lns1) + exp(2*lns2))^2 + 2*exp(2*lns2)/(exp(2*lns1) + exp(2*lns2))
se <- sqrt(var11*deriv1^2 + var22*deriv2^2 + 2*cov21*deriv1*deriv2)
lcl95.no <- h2 - qnorm(.975)*se
ucl95.no <- h2 + qnorm(.975)*se
se.qnorm <- se/dnorm(qnorm(h2))
lcl95.qnorm <- pnorm(qnorm(h2) - qnorm(.975)*se.qnorm)
ucl95.qnorm <- pnorm(qnorm(h2) + qnorm(.975)*se.qnorm)
diff <- lns1 - lns2
se.diff <- sqrt(var11 + var22 - 2*cov21)
diff.lower <- diff - 1.96*se.diff
diff.upper <- diff + 1.96*se.diff
h2.upper <- 1/(exp(2*diff.lower) + 1)
h2.lower <- 1/(exp(2*diff.upper) + 1)
set.seed(seed)
mu = c(lns1, lns2)
Sigma = matrix(c(var11, cov21, cov21, var22), ncol = 2)
samp <- exp(2*MASS::mvrnorm(n.sim, mu, Sigma, empirical = TRUE))
samp.h2 <- samp[,2]/(samp[,1] + samp[,2])
CI <- quantile(samp.h2, c(0.025, 0.975))
CL <- matrix(c(lcl95.no,ucl95.no, lcl95.qnorm,ucl95.qnorm, h2.lower,h2.upper, CI[1],CI[2]),ncol=2,byrow=TRUE)
if(verbose)
{
print(summary(res))
cat("\nHeritability:", h2,
"\n95% CI (no transformation): ", CL[1,1], "-", CL[1,2],
"\n95% CI (probit transformation): ", CL[2,1], "-", CL[2,2],
"\n95% CI (simpler transformation):", CL[3,1], "-", CL[3,2],
"\n95% CI (simulated): ", CL[4,1], "-", CL[4,2], "\n")
}
invisible(list(lme.result=res,h2=h2,CL=CL))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.