首页 > 解决方案 > 将误差带添加到线性混合模型 (lme4)

问题描述

我有一个带有分类(物种)和连续(温度)预测器的混合模型(lme4)。我在不同的地块对数据进行了采样,分布在 3 个样带上。我正在寻找一种方法来围绕我的模型预测绘制一些误差带。在此处输入图像描述

这里有一些示例数据(我是模拟数据的新手):

library(lme4)
library(tidyverse)    

Transect <- rep(c("A", "B", "C"), each = 15, length.out = 135)
Plot <- rep(1:45, each = 3)
Species <- rep(c("D", "E", "F"), each = 45)
    
Temp1 <- (seq(5,10, length.out = 15)) + rnorm(n=15)
Temp2 <- (seq(5,10, length.out = 15)) + rnorm(n=15)
Temp3 <- (seq(5,10, length.out = 15)) + rnorm(n=15)
Temp4 <- (seq(5,10, length.out = 15)) + rnorm(n=15)
Temp5 <- (seq(5,10, length.out = 15)) + rnorm(n=15)
Temp6 <- (seq(5,10, length.out = 15)) + rnorm(n=15)
Temp7 <- (seq(5,10, length.out = 15)) + rnorm(n=15)
Temp8 <- (seq(5,10, length.out = 15)) + rnorm(n=15)
Temp9 <- (seq(5,10, length.out = 15)) + rnorm(n=15)
Temp <- c(Temp1, Temp2, Temp3, Temp4, Temp5, Temp6, Temp7, Temp8, Temp9)
    
Onset1 <- (seq(90,160, length.out = 15)) + rnorm(n=15, sd = 20)
Onset2 <- (seq(90,160, length.out = 15)) + rnorm(n=15, sd = 20)
Onset3 <- (seq(90,160, length.out = 15)) + rnorm(n=15, sd = 20)
Onset4 <- (seq(90,160, length.out = 15)) + rnorm(n=15, sd = 20) + 20
Onset5 <- (seq(90,160, length.out = 15)) + rnorm(n=15, sd = 20) + 20
Onset6 <- (seq(90,160, length.out = 15)) + rnorm(n=15, sd = 20) + 20
Onset7 <- (seq(90,160, length.out = 15)) + rnorm(n=15, sd = 20) + 40
Onset8 <- (seq(90,160, length.out = 15)) + rnorm(n=15, sd = 20) + 40
Onset9 <- (seq(90,160, length.out = 15)) + rnorm(n=15, sd = 20) + 40
Onset <- c(Onset1, Onset2, Onset3, Onset4, Onset5, Onset6, Onset7, Onset8, Onset9)
    
data <- data.frame(Onset, Transect, Plot, Species, Temp)

我拟合了一个线性混合模型:

model <- lmer(Onset ~ Temp * Species + (1|Transect/Plot), data = data)

我想为每个物种都有一条回归线,所以我写了一个小函数来提取每个物种的系数,因为我不确定是否有其他方法可以做到这一点:

extract.coef <- function(model){
  coefficients.model <<- as.numeric(fixef(model))
  intercept.D <<- coefficients.model [1]
  intercept.E <<- intercept.D + coefficients.model[3]
  intercept.F <<- intercept.D + coefficients.model[4]
  slope.D <<- coefficients.model [2]
  slope.E <<- slope.D + coefficients.model[5]
  slope.F <<- slope.D + coefficients.model[6]
  start.D <<- as.numeric(data %>% filter(Species == "D") %>% arrange(Temp) %>% slice(n = 1) %>% select ("Temp"))
  start.E <<- as.numeric(data %>% filter(Species == "E") %>% arrange(Temp) %>% slice(n = 1) %>% select ("Temp"))
  start.F <<- as.numeric(data %>% filter(Species == "F") %>% arrange(Temp) %>% slice(n = 1) %>% select ("Temp"))
  End.D <<- as.numeric(data %>% filter(Species == "D") %>% arrange(desc(Temp)) %>% slice(n = 1) %>% select ("Temp"))
  End.E <<- as.numeric(data %>% filter(Species == "E") %>% arrange(desc(Temp)) %>% slice(n = 1) %>% select ("Temp"))
  End.F <<- as.numeric(data %>% filter(Species == "F") %>% arrange(desc(Temp)) %>% slice(n = 1) %>% select ("Temp"))
}

我提取了系数:

extract.coef(model)

并绘制模型:

ggplot(data, aes(Temp, Onset, colour = Species, shape = Transect)) +
  geom_point(size = 0.9, position = position_jitter (width = 0.1, height = 2)) +
  geom_segment(aes(x = start.D, xend = End.D, y = intercept.D + slope.D * start.D, yend = intercept.D + slope.D * End.D), colour = "red", size = 0.4) +
  geom_segment(aes(x = start.E, xend = End.E, y = intercept.E + slope.E * start.E, yend = intercept.E + slope.E * End.E), colour = "green", size = 0.4) +
  geom_segment(aes(x = start.F, xend = End.F, y = intercept.F + slope.F * start.F, yend = intercept.F + slope.F * End.F), colour = "blue", size = 0.4) +
  theme_bw() +
  scale_shape_manual(labels = c("A", "B", "C"), values = c(16,17,15), name = "Transect:") +
  scale_colour_manual(values = c( "red", "green", "blue"), labels = c("D", "E", "F"), name = "Species:") +
  labs (x = "Temperature [°C]", y = "DOY", colour = "Species", shape = "Transect") +
  theme(axis.title.x = element_text(size = 10), axis.title.y = element_text(size = 10), axis.text.x = element_text(size =8), axis.text.y = element_text(size =8)) +
coord_cartesian(xlim = c(3, 12), expand = TRUE)

有没有办法在我的回归线周围绘制误差带?

标签: rggplot2lme4

解决方案


推荐阅读