Note: This might be more appropriate for Cross Validated, and I'm fine moving it over there if requested, but I thought I would try here first since it could be r-related.
I am comparing growth curve parameter estimates between two groups using nonlinear mixed models with nlme. I am using boot_nlme() from nlraa to bootstrap the confidence intervals for the model parameter estimates and predicted data.
Below is the code I am using but with the penguin.data from FlexParamCurve. It doesn't have the same issue exactly, but provides the code in a usable format. If this isn't enough, I can edit with my own data.
For my personal data, my issue stems from the predicted data's CIs overlapping (i.e. showing no statistical difference) at the asymptote when according to the model and bootstrapped parameter estimate CIs, they do not overlap (i.e. showing statistical difference).
My question is: Why would the bootstrapped CIs of the model parameter estimates dispute the boostrapped predicted data's CIs from the same model? Could this be me not fully understanding the bootstrap method? Or could there be something else I am missing (code-related)?
Any help/insight would be appreciated.
library(FlexParamCurve) #loads dataset and package 'nlme'
library(ggplot2)
library(nlraa)
library(boot)
library(dplyr)
library(car)
VarFunc.Auto<-varPower(form=~fitted(.))
##creating model
richards.func <- function(age, A, Ti, k, d){
A * (1 + (d - 1) * exp(( (-k) * (age - Ti)) / ( d ^ ( d / (1 - d))))) ^ (1 / (1 - d))
}
ggplot(penguin.data, aes (x = ckage, y = weight, color = year)) +
geom_smooth()
#fixed asymptote
mod <- nlme(weight ~ SSlogis(ckage, Asym, R0, lrc),
data = penguin.data,
fixed= list(Asym ~ year,
R0 ~ year,
lrc ~ year),
random = Asym ~ 1,
start = c(Asym = 1000, 0,
R0 = 21, 0,
lrc = 1, 0),
control = list(maxIter = 100),
na.action = na.pass)
summary(mod)
peng.mdiff <- function(x, ckage = seq(0, 80, length.out = 500)){
ndat <- expand.grid(ckage = ckage,
year = c('2000','2002'),
nest = NA,
bandid = NA,
stringsAsFactors = TRUE)
prd <- predict(x, newdata = ndat, level = 0)
}
set.seed(123)
#this takes a few minutes to run on my computer
system.time(peng.bt <- boot_nlme(mod, peng.mdiff, R = 500, cores = 3))
set.seed(123)
system.time(peng.bt.ci <- confint(peng.bt, level = 0.95))
prd.df <- as.data.frame(peng.bt.ci)
prd.df <- prd.df %>%
rename_at('2.5 %', ~'lower.ci') %>%
rename_at('97.5 %', ~'upper.ci')
ndat1 <- with(penguin.data,
expand.grid(
ckage=seq(0, 80 , length.out=500),
year = c('2000','2002'),
nest = NA,
bandid = NA
))
newX <- ndat1 %>%
mutate(prd = predict_nlme(mod, newdata = ndat1, level = 0))
comb.df1 <- cbind(prd.df, ndat1, 'prd' = newX$prd)
ci.plot <- ggplot() +
geom_ribbon(comb.df1, mapping = aes(x = ckage, ymin = lower.ci, ymax = upper.ci, fill = year ), alpha = 0.50) +
geom_line(data = newX, aes(x = ckage, y = prd, group = year), color = 'black', alpha = 0.75)
ci.plot
For my personal data, my issue stems from the predicted data's CIs overlapping (i.e. showing no statistical difference) at the asymptote when according to the model and bootstrapped parameter estimate CIs, they do not overlap (i.e. showing statistical difference).
Overlapping confidence intervals does not imply no statistically significant difference. Confidence intervals can overlap and there can still be a statistically significant difference. That is, 95% confidence intervals can still overlap and the difference be significantly different, at the 5% level.