How to include heteroscedasticity in the predict function of a mixed model in lme

Question

Helloy everybody,

I plotted a mixed model from "nlme" with "ggplot2" and I would like to consider different variance in my category variable in the model (weights=varIdent(form = ~1 | category)). When I checked the plot I realised that the model prediction in the plot did no include the heteroscedasticity, altough the model does. How can I adjust the predict function, so that includes the heteroscedasticity?

create some data:

no <- 1:20
type <- factor("A","B")
set.seed(1)
baseline <- rnorm(20, 30, 10)
set.seed(2)
event <- rnorm(20, 60, 30)
set.seed(3)
postevent1 <- rnorm(20, 40, 20)
set.seed(4)
postevent2 <- rnorm(20, 30, 10)
set.seed(5)
postevent3 <- rnorm(20, 20, 5)
data <- data.frame(no, baseline, event, postevent1, postevent2, postevent3)
data$type[data$no %% 2 == 0] <- "A"
data$type[data$no %% 2 != 0] <- "B"
data$event[data$type == "A"] <- data$event[data$type == "A"] -2
data$no <- factor(data$no)
data$type <- factor(data$type)

transform the data to the long format

library(dplyr)

long <- data %>% gather(key = category, value = measure, -no, -type)

create the model

library(nlme)
model <- lme(measure ~ category*type, random = ~ 1|no, data= long,   
weights=varIdent(form = ~1 | category), method = "REML")
new <- long %>% select(-measure)
long$pred <- predict(model, newdata = new)

plot it

library(ggplot2)
ggplot(long, aes(x=category, y =  measure, colour = type,    
group=interaction(type,category))) +
geom_point() +
facet_grid(~ type) +
geom_line(aes(y=pred), size=0.8, colour = "black") +
theme(axis.text.x = element_text(angle=30, hjust=1, vjust=1))

Answer 1

This is my solution using the "AICcmodavg" package. If anyone finds a better solution, please let me know.

predicting the values

library(AICcmodavg)
pred <- predictSE.lme(model, newdata = new, print.matrix = TRUE)
pred <- data.frame(pred)
newlong <- data.frame(long, pred)

plotting the mean with the 95% confidence interval from the prediction

ggplot(newlong, aes(x= category, y= measure)) + 
geom_point(aes(x = category, y = fit)) + 
facet_grid(~ type) + 
geom_point(aes(colour = type), alpha = 0.5, size = 3)  + 
geom_errorbar(aes(ymax=fit + se.fit*1.96, ymin=fit - se.fit*1.96), 
position=position_dodge(0.9), width=0.25) +
theme(axis.text.x = element_text(angle=30, hjust=1, vjust=1))