estimate.lmm {LMMstar} | R Documentation |
Delta Method for Mixed Models
Description
Perform a first order delta method
Usage
## S3 method for class 'lmm'
estimate(
x,
f,
df = !is.null(x$df),
robust = FALSE,
type.information = NULL,
level = 0.95,
method.numDeriv = NULL,
average = FALSE,
transform.sigma = NULL,
transform.k = NULL,
transform.rho = NULL,
...
)
Arguments
x |
a |
f |
[function] function of the model coefficient computing the parameter(s) of interest. Can accept extra-arguments. |
df |
[logical] Should degree of freedom, computed using Satterthwaite approximation, for the parameter of interest be output. Can also be a numeric vector providing providing the degrees of freedom relative to each estimate. |
robust |
[logical] Should robust standard errors (aka sandwich estimator) be output instead of the model-based standard errors. |
type.information |
[character] Should the expected information be used (i.e. minus the expected second derivative) or the observed inforamtion (i.e. minus the second derivative). |
level |
[numeric,0-1] the confidence level of the confidence intervals. |
method.numDeriv |
[character] method used to approximate the gradient: either |
average |
[logical] is the estimand the average output of argument |
transform.sigma |
[character] Transformation used on the variance coefficient for the reference level. One of |
transform.k |
[character] Transformation used on the variance coefficients relative to the other levels. One of |
transform.rho |
[character] Transformation used on the correlation coefficients. One of |
... |
extra arguments passed to |
Examples
if(require(lava) && require(nlme)){
#### Random effect ####
set.seed(10)
dL <- sampleRem(1e2, n.times = 3, format = "long")
e.lmm1 <- lmm(Y ~ X1+X2+X3 + (1|id), repetition = ~visit|id, data = dL)
nlme::ranef(e.lmm1, se = TRUE)
e.ranef <- estimate(e.lmm1, f = function(p){nlme::ranef(e.lmm1, p = p)})
e.ranef
if(require(ggplot2)){
df.gg <- cbind(index = 1:NROW(e.ranef), e.ranef)
gg.ranef <- ggplot(df.gg, aes(x = index, y=estimate, ymin=lower, ymax = upper))
gg.ranef + geom_point() + geom_errorbar() + ylab("estimated random effect") + xlab("id")
}
#### ANCOVA via mixed model ####
set.seed(10)
d <- sampleRem(1e2, n.time = 2)
e.ANCOVA1 <- lm(Y2~Y1+X1, data = d)
if(require(reshape2)){
dL2 <- melt(d, id.vars = c("id","X1"), measure.vars = c("Y1","Y2"),
value.name = "Y", variable.name = "time")
dL2$time <- factor(dL2$time, levels = c("Y1","Y2"), labels = c("1","2"))
## estimated treatment effect (no baseline constraint)
e.lmm <- lmm(Y ~ time + time:X1, data = dL2, repetition = ~time|id)
e.delta <- estimate(e.lmm, function(p){
c(Y1 = p["rho(1,2)"]*p["k.2"],
X1 = p["time2:X1"]-p["k.2"]*p["rho(1,2)"]*p["time1:X1"])
}) ## same estimate and similar standard errors.
e.delta ## Degrees of freedom are a bit off though
cbind(summary(e.ANCOVA1)$coef, df = df.residual(e.ANCOVA1))
## estimated treatment effect (baseline constraint)
dL2$time2 <- as.numeric(dL2$time=="2")
e.lmmC <- lmm(Y ~ time2 + time2:X1, data = dL2, repetition = ~time|id)
e.deltaC <- estimate(e.lmmC, function(p){
c(Y1 = p["rho(1,2)"]*p["k.2"],
X1 = p["time2:X1"])
})
e.deltaC ## Degrees of freedom are a bit more accurate
}
}