1

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))

enter image description here

Kev
  • 425
  • 3
  • 8
  • 1
    If you're willing to switch to `lme4` instead of `nlme`, the `sim` function might do what you're looking for. – Gregor Thomas Dec 16 '15 at 22:17
  • The variance of A and B is the same, but in the categorie variable the variance differs. The lme model takes this into account, but the predict function from the lme model does not. – Kev Dec 17 '15 at 10:20
  • Dear Gregor, you are right! I have already changed the title. – Kev Dec 17 '15 at 16:50
  • I don't know how to account for heteroscedasticity in lme4. The comment I have found in the internet said that the weight argument is not wotking properly in lmer [link](http://r.789695.n4.nabble.com/weights-argument-in-the-lmer-function-in-lme4-td792717.html) – Kev Dec 17 '15 at 17:30
  • Once again your right. I have just fixed it. – Kev Dec 17 '15 at 22:02

1 Answers1

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))

enter image description here

Kev
  • 425
  • 3
  • 8