R:将组和单个多项式趋势线添加到 GMM 图

R: Adding group and individual polynomial trend lines to a GMM plot

我正在为如何将个人和群体趋势线添加到我的地块而苦苦挣扎。 (R 并使用 ggplot2)。

这是我正在使用的代码:

MensHG.fm2=lmer(HGNewtons~Temperature+QuadTemp+Run+(1|Subject),MenstrualData) #model

plot.hg<-data.frame(MensHG.fm2@frame,fitted.re=fitted(MensHG.fm2))

g1<-ggplot(plot.hg,aes(x=Temperature,y=HGNewtons))+geom_point()

g2<-g1+facet_wrap(~Subject, nrow=6)+ylab(bquote('HG MVF (N)'))+xlab(bquote('Hand ' ~T[sk] ~(degree*C)))

g3<-g2+geom_smooth(method="glm", formula=y~ploy(x,2), se=FALSE) #This gives me my individual trendlines

现在我想将数据的 g1 部分的趋势线(即总体趋势)放在我的每个单独的地块上 - 执行此操作的最佳方法是什么?如果我使用代码,我可以看到趋势:

g5=g1+geom_smooth(method="glm", formula=y~poly(x,2), se=FALSE)

但是,一旦我执行了 facet-wrap(我得到与 g3 相同的输出),这条趋势线就消失了

使用g4<-g3+geom_smooth(data=MensHG.fm2)

似乎无法解决问题

没有你的数据的最小工作示例,我使用了内置的 iris 数据。为了演示,我在这里假装 物种 是不同的主题。

library(lme4)
library(ggplot2)

fit.iris <- lmer(Sepal.Width ~ Sepal.Length + I(Sepal.Length^2) + (1|Species), data = iris)

为了简单起见,我还使用了两个额外的包,broomdplyrbroom 中的 augment 与上面 ..., fitted.re=fitted(MensHG.fm2) 所做的相同,但有一些额外的功能。我也使用 dplyr::select,但这并不是严格需要的,具体取决于您想要的输出(图 2 与图 3 之间的差异)。

library(broom)
library(dplyr)   
augment(fit.iris)
# output here truncated for simplicity
  Sepal.Width Sepal.Length I.Sepal.Length.2. Species  .fitted       .resid   .fixed ...
1         3.5          5.1             26.01  setosa 3.501175 -0.001175181 2.756738 
2         3.0          4.9             24.01  setosa 3.371194 -0.371193601 2.626757 
3         3.2          4.7             22.09  setosa 3.230650 -0.030649983 2.486213 
4         3.1          4.6             21.16  setosa 3.156417 -0.056417409 2.411981 
5         3.6          5.0                25  setosa 3.437505  0.162495354 2.693068 
6         3.9          5.4             29.16  setosa 3.676344  0.223656271 2.931907 
ggplot(augment(fit.iris), 
       aes(Sepal.Length, Sepal.Width)) + 
  geom_line(#data = augment(fit.iris) %>% select(-Species), 
            aes(y = .fixed, color = "population"), size = 2) +
  geom_line(aes(y = .fitted, color = "subject", group = Species), size = 2) +
  geom_point() + 
  #facet_wrap(~Species, ncol = 2) +
  theme(legend.position = c(0.75,0.25))

请注意,我 # 注释了两个语句:data = ...facet_wrap(...)。注释掉这些行后,您将得到如下输出:

您在整个范围内对总体进行了平滑处理(.fixed 对于固定效应),然后您进行了显示拟合模型值 (.fitted) 的组平滑处理,同时考虑到主题级截取。

然后您可以通过删除代码段中的第二个 #-注释标记来在方面显示它:

这是相同的,但由于拟合值仅存在于每个主题级别面板的原始数据范围内,因此总体平滑被截断到该范围。

为了解决这个问题,我们可以删除第一个 #-注释标记: