R coxme: How to get study-specific treatment effect and 95% confidence interval from a mixed effect model?

Question

I am using the coxme to fix a mixed effect Cox model. I have two random effects in the model (intercept [baseline hazard] and slope of the treatment). After running the model, I got the fixed coefficients and random effects output:

Fixed coefficients
                                 coef exp(coef)   se(coef)     z    p
treatment                 0.420463710 1.5226675 0.28313885  1.49 0.14

Random effects
     Group  Variable          Std Dev    Variance   Corr      
     Study  Intercept         1.3979618  1.9542971 -0.6995296
            treatment         0.4815433  0.2318839

I also used the ranef to get the random effect coefficients for each study. Outpout:

         Intercept    treatment
Study1   -1.591497    0.5479255
Study2    1.276003   -0.4392570
Study3   -1.051512    0.3621938
Study4    1.367005   -0.4708623

I was wondering how I can obtain the treatment effect (coefficient) and 95% CI for each study separately. Can I get the point estimate by summing up the overall coefficient of fixed effect and study-specific random effect (e.g. 0.420463710 + 0.5479255 for Study1)? But what about the 95%CI?

Is there any way to get the study-specific treatment effect and 95% CI (and the study's weight), like in ipdforest in Stata?

Thank you very much.

I'll say for the record (*much* belatedly) that there's not an obvious way to get conditional standard deviations of the random effects from `lme`: the mathematics in `vignette("lmer", package = "lme4")` applies, but one would have to figure out how to extract the quantities `V` and `Lambda_{theta}` from an `lme` object ... — Ben Bolker, May 15 '23 at 17:06

score 0 · Answer 1 · answered May 18 '23 at 00:06

Here is a not-completely-worked out answer, with some caveats and commentary:

The general recipe is that if you have a set of linear combinations of parameters (fixed and random) you want to apply — combine that into a matrix M — along with a set of parameters (fixed and random) B and a covariance matrix (for the same parameters) S, then the predictions are M %*% B and the variances of the predictions are diag(M %*% S %*% t(S)) (there are slightly more efficient ways to compute this ...)

The coxme.condsd function below generates what I would call conditional standard deviations for a given set of predictor variables. This isn't quite the same as generating the conditional standard deviations of the slopes; that would involve a different set of linear combinations. But it's at least the beginning of an answer (this is also posted here, and will hopefully be incorporated into coxme at some point in the future ...)

coxme.condsd <- function(fit,
                         data = eval(getCall(fit)$data),
                         grpvar,
                         rform,
                         sd.only = TRUE) {
    ## model.frame() to get a sanitized (NA-free) data set
    vf <- all.vars(fit$formula$fixed)
    ## pick apart random effects terms
    vr <- (fit$formula$random
        |> lapply(function(x) lapply(lme4::findbars(x), all.vars))
        |> unlist(recursive = TRUE)
    )
    mf <- model.frame(reformulate(c(vf, vr)), data)
    ## indicator matrix for grouping variable
    Fmat <- t(Matrix::fac2sparse(mf[[grpvar]]))
    ## Z matrix
    Xr <- model.matrix(rform, mf)
    ## was (special case): Z <- cbind(Fmat, Fmat*mf$ph.ecog)
    Z <- t(Matrix::KhatriRao(t(Xr), t(Fmat)))
    ## remove intercept from RHS of fixed formula
    fform <- update(fit$formula$fixed[-2], ~ 0 + .)
    X <- model.matrix(fform, data = mf)
    ## respect order of vars in cov matrix: RE folowed by FE coefs
    ZX <- cbind(Z,X)
    V <- as.matrix(fit$variance)
    V2 <- ZX %*% V %*% t(ZX)
    if (sd.only) sqrt(diag(V2)) else V2
}

## cc <- c(ranef(fit)$inst[,1], ranef(fit)$inst[,2], coef(fit))
## head(drop(ZX %*% cc))
## head(fit$linear.predictor)
fit <- coxme(Surv(time, status) ~ ph.ecog + age + (1 + ph.ecog|inst), lung)

## this gets conditional SD (including both fixed & random effect uncertainty,
##  but using a plug-in estimate of the RE *variances*) for each observation in
## the original data set
coxme.condsd(fit, grpvar = "inst", rform = ~1 + ph.ecog)

R coxme: How to get study-specific treatment effect and 95% confidence interval from a mixed effect model?

1 Answers1