复制混合模型的点阵图

Replicating lattice graph for a mixed model

我正在尝试使用 ggplot2 为混合模型复制点阵图。我的 ggplot 图看起来非常相似,但我不确定是否适合黄土线模型。

我的目标是使用 ggplot2 从混合模型中添加黄土线。以下是我的命令示例:

library(nlme)
library(ggplot2)
library(lattice)
library(lme4)

data(MathAchieve)
attach(MathAchieve)
mses <- tapply(SES, School, mean) 
mses[as.character(MathAchSchool$School[1:10])] 

Bryk <- as.data.frame(MathAchieve[, c("School", "SES", "MathAch")])

names(Bryk) <- c("school", "ses", "mathach")
sample20 <- sort(sample(7185, 20)) # 20 randomly sampled students

Bryk$meanses <- mses[as.character(Bryk$school)]

Bryk$cses <- Bryk$ses - Bryk$meanses
sector <- MathAchSchool$Sector
names(sector) <- row.names(MathAchSchool)
Bryk$sector <- sector[as.character(Bryk$school)]

attach(Bryk)

cat <- sample(unique(school[sector=="Catholic"]), 20)
Cat.20 <- groupedData(mathach ~ ses | school,  data=Bryk[is.element(school, cat),])

格子图:

trellis.device(color=T)
xyplot(mathach ~ ses | school, data=Cat.20, main="Catholic", 
       panel=function(x, y) {
         panel.loess(x, y, span=1) 
         panel.xyplot(x, y)
         panel.lmline(x, y, lty=2)
       })

使用 ggplot 绘制图表:

ggplot(Cat.20, aes(x = ses, y =mathach )) + 
  geom_point(size=1, shape=1) + 
  stat_smooth(method="lm",se=F)+
  stat_smooth(, colour="Red",se=F)+
  facet_wrap(school~., scale = "free_y")  

如有任何建议,我们将不胜感激。

前言

在开始解释之前,请允许我向您介绍这个问题:Why is it not advisable to use attach() in R, and what should I use instead?

虽然建议您使问题可重现,但您使用的代码可以进行一些清理。例如:

  1. 不要包含代码中未使用的包(我没有看到需要 lme4 包);
  2. 不需要使用data(...)来加载MathAchieve。有关详细信息,请参阅 ?data 中的 "Good Practices" 部分。
  3. 如上所述,不要使用attach()
  4. 为了完全再现,在任何随机抽样之前使用 set.seed()
  5. 对于最小的示例,不要绘制 20 所学校,而较小的数量就可以。

由于您使用的是 tidyverse 软件包之一进行绘图,我推荐使用其集合中的另一个软件包进行数据处理:

library(nlme)
library(ggplot2)
library(lattice)
library(dplyr)

Bryk <- MathAchieve %>%
  select(School, SES, MathAch) %>%
  group_by(School) %>%
  mutate(meanses = mean(SES),
         cses = SES - meanses) %>%
  ungroup() %>%
  left_join(MathAchSchool %>% select(School, Sector),
            by = "School")
colnames(Bryk) <- tolower(colnames(Bryk))

set.seed(123)
cat <- sample(unique(Bryk$school[Bryk$sector == "Catholic"]), 2)
Cat.2 <- groupedData(mathach ~ ses | school,
                     data = Bryk %>% filter(school %in% cat))

说明

说完了,我们来看看loess的相关函数:

来自 ?panel.loess:

panel.loess(x, y, span = 2/3, degree = 1,
            family = c("symmetric", "gaussian"),
            ... # omitted for space
            )

来自 ?stat_smooth:

stat_smooth(mapping = NULL, data = NULL, geom = "smooth",
  method = "auto", formula = y ~ x, span = 0.75, method.args = list(), 
  ... # omitted for space
  )

其中 method = "auto" 默认为 stats 包中的 loess,用于 <1000 次观察。

来自 ?loess:

loess(formula, data, span = 0.75, degree = 2,
      family = c("gaussian", "symmetric"),
      ... #omitted for space
      )

简而言之,黄土地块的默认参数是 lattice 包的 span = 2/3, degree = 1, family = "symmetric"ggplot2 包的 span = 0.75, degree = 2, family = "gaussian"如果您希望结果图匹配:

,您必须指定匹配参数
xyplot(mathach ~ ses | school, data = Cat.2, main = "Catholic", 
       panel=function(x, y) {
         panel.loess(x, y, span=1, col = "red")  # match ggplot's colours
         panel.xyplot(x, y, col = "black")       # to facilitate comparison
         panel.lmline(x, y, lty=2, col = "blue")
       })

ggplot(Cat.2, aes(x = ses, y = mathach)) + 
  geom_point(size = 2, shape = 1) +
  stat_smooth(method = "lm", se = F)+
  stat_smooth(span = 1,
              method.args = list(degree = 1, family = "symmetric"),
              colour = "red", se = F)+
  facet_wrap(school ~ .) +
  theme_classic() # less cluttered background to facilitate comparison