4
$\begingroup$

I have a question regarding confidence intervals from a non-linear mixed effects model. I ran the following model on my data:

#Set up model 
ModelFunc <- function(t, A, B, C) { A + B * t ^ C }

ModelGradient <- deriv(
  body(ModelFunc)[[2]], 
  namevec = c("A" ,"B", "C"), 
  function.arg=ModelFunc )

model1 <- nlmer(
  Resp ~ ModelGradient(t=time, A, B, C) ~ (A | ID),  
  data = DF,
  start = c(A=60, B=-0.003, C=3.1)) #not reproducible because I can't post the data, but I'm posting the model output below

It gave me the following fixed effects and covariance matrix:

model1_fixef <- structure(c(60.4587890233239, -0.00172745908090093, 3.07643744499072), .Names = c("A", "B", "C"))

model1_vcov <- new("dpoMatrix", x = c(1.65646896448117, -0.000137861599380447, -0.0228735238880192, 
                       -0.000137861599380447, 6.11544583017694e-07, 0.00011478199667793, 
                       -0.0228735238880192, 0.00011478199667793, 0.021782162709204),
               Dim = c(3L, 3L),
               Dimnames = list(c("A", "B", "C"), c("A", "B", "C")))

To compute the 95% confidence interval, I did the following:

#Resample 10000 parameter estimates from fixed effects and cov matrix:
pars.resamp <- mvrnorm(100000, mu = model1_fixef , Sigma = model1_vcov)

#Calculate y values using resampled parameters
xvals <- seq(1, 30, length.out = 100)
yvals <- apply(pars.resamp, 1, function(x)  ModelFunc(xvals, x[1], x[2], x[3]))

#determine 2.5 and 97.5% quantiles of the 10000 y values at each x value
dfCI <- data.frame(t(apply(yvals, 1, quantile, c(0.025,0.975))) )

#Combine output in data frame
output <- data.frame(x=xvals,
                     y=ModelFunc(xvals, model1_fixef[1], model1_fixef[2], model1_fixef[3]),
                     lower=dfCI$"X2.5.",
                     upper=dfCI$"X97.5.")

#Plot
ggplot(output, aes(x=xvals))+geom_line(aes(y=y))+
  geom_ribbon(aes(ymin=lower, ymax=upper), fill="dodgerblue3", alpha=0.3)

This gives me the following: enter image description here

Published confidence intervals are usually symmetric, but this one is clearly not. The reason for the asymmetry ( I think ) is that the B parameter is really small and becomes bigger than 0 in some of the resampled parameters, resulting in an upward curve. The exponential model gives a better fit than a simple linear model and based on prior literature we expect an exponential relationship between our variables. Is it fine to report the model like this? Or is there a better way to compute the confidence intervals?

$\endgroup$
1
  • 1
    $\begingroup$ (Trying to) get a symmetric confidence interval is just a construction decision. There is no intrinsic reason why confidence intervals should be symmetric in general. And even for symmetric distributions, unsymmetric confidence intervals are well known, e.g. one-sided c.i.s. Deciding to construct a symmetric confidence interval is just a decision to get rid of the ambiguity that for a given distribution often many intervals achieve the requested confidence level. This is reduced to "the" confidence interval by making it symmetric (if possible), as narrow as possible, one-sided, etc. $\endgroup$ May 29 '18 at 8:51
0
+50
$\begingroup$

I think you are simply capturing the effects of the joint variability of the parameter estimates. That is to say, if you plot the different options for $B \pm \text{ME}_b$ and for the different options for $C \pm \text{ME}_c$, then you would obtain unsymmetrical CIs, particularly as the t gets quite large.

Here is some sample code to help demonstrate this:

xvals <- seq(0,30,length.out=100)
MEb <- .0012
MEc <- .10

A <- 60.5
B <- -0.002
C <- 3.076
Blo <- B - MEb
Bhi <- B + MEb
Clo <- C - MEc
Chi <- C + MEc

## add ME to just the parameter B
yvals <- A+B*xvals^C
yvalslo <- A+Blo*xvals^C
yvalshi <- A+Bhi*xvals^C

plot(xvals,yvals,type="l",lwd=3)
lines(xvals,yvalslo,lty=2)
lines(xvals,yvalshi,lty=2)

## add ME to just the parameter C
yvals <- A+B*xvals^C
yvalslo <- A+B*xvals^Clo
yvalshi <- A+B*xvals^Chi

plot(xvals,yvals,type="l",lwd=3)
lines(xvals,yvalslo,lty=2)
lines(xvals,yvalshi,lty=2)

## add ME to both the parameters B & C
yvals <- A+B*xvals^C
yvalslo <- A+Blo*xvals^Clo
yvalshi <- A+Bhi*xvals^Chi

plot(xvals,yvals,type="l",lwd=3)
lines(xvals,yvalslo,lty=2)
lines(xvals,yvalshi,lty=2)

In my humble opinion, I believe the graphic representation you've provided is a very compelling way to report the confidence intervals for a model estimate such as this.

$\endgroup$
2
  • $\begingroup$ Thank you for your answer. An additional question: my colleague used SAS to get 95% confidence intervals for the same model. It seems that SAS used the delta method, which (apparently) gives symmetric CIs. I'm not so familiar with the mathematical details and I'm wondering what to use and why. Also, does such a C.I. give any indication that the model does not fit the data very well? $\endgroup$ May 30 '18 at 18:34
  • 1
    $\begingroup$ I believe the CI's are symmetric for the parameter estimates, but they would not be symmetric when converted back to the original scale. If you notice, the margin of errors in my code are symmetric...but their effect on the final graphs are not symmetric. $\endgroup$
    – Gregg H
    May 30 '18 at 22:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.