如何可视化线性混合效应模型预测

How to visualize linear mixed-effects model predictions

我有一个非常简单的线性混合效应模型:

Rodlangde 表示根长 Mehlich意为植物有效磷 Lokalitet 表示地区

model<-lme(Rodlangde~Mehlich,random=~1|Lokalitet)

我非常想制作一个图,您可以在其中看到 10 个(我有 10 个地点)不同的线性图,它们具有相同的斜率但模型由不同的截距组成。我现在已经尝试在网上搜索解决方案2天了,但是我能找到的那些代码太复杂了我看不懂,或者我找不到我需要哪些包才能使用这些代码。谁能帮我用一个简单的代码来可视化同一个图中的 10 个不同的图形?

亲切的问候

Data:
Lokalitet   pH  Mehlich Kol Rodlangde
Odrup   7.02    0   0.919642857 6.362373845
Odrup   6.87    0   0.875   5.372476457
Odrup   7.09    0   0.868421053 14.23942792
Odrup   6.64    0   0.939393939 4.640122704
Orhoje  6.81    12.13896843 0.83    2.842893319
Orhoje  7.44    7.062027912 0.741666667 4.399108137
Orhoje  7.23    6.915193254 0.917355372 3.597793514
Orhoje  6.73    3.930162033 0.909090909 5.28750758
Melby   5.2 28.20132199 0.669642857 2.541898484
Melby   5.35    14.97459413 0.519685039 2.790724941
Melby   5.04    8.352860756 0.596153846 2.927228501
Melby   5.02    10.51701575 0.596153846 1.538074359
Kallingedal 8.4 17.47092431 0.458646617 8.059178499
Kallingedal 8.33    21.74560339 0.703703704 10.50345245
Kallingedal 8.3 21.34370501 0.762295082 7.610537154
Kallingedal 8.37    25.06114498 0.770491803 11.88896483
Ravrigtigkalk   5.61    5.117349119 0.952380952 9.307948512
Ravrigtigkalk   5.92    3.400217532 0.85046729  12.80110763
Ravrigtigkalk   5.77    3.358878819 0.607476636 14.82758346
Ravrigtigkalk   5.82    2.854552095 0.9375  4.231563699
Karsemose   5.28    0   0.813084112 12.06213863
Karsemose   5.36    1.312479611 0.838095238 7.341806594
Karsemose   5.32    0   0.898148148 10.1038273
Karsemose   5.34    0   0.821782178 8.16704508
Lergraven   8.43    44.62536835 0.847457627 13.48193914
Lergraven   8.41    39.52348256 0.884297521 12.67270404
Lergraven   8.39    43.26503035 0.880597015 21.24738813
Lergraven   8.41    40.8293479  0.770491803 16.12249983
Hvidtjorn   7.98    41.68676311 0.923076923 19.46781449
Hvidtjorn   8.16    43.89098256 0.827868852 14.39349303
Hvidtjorn   8.19    37.35675233 0.942857143 34.98582813
Hvidtjorn   8.2 29.90406084 0.927927928 17.09668084
Ravsurt 5.17    5.061969924 0.821782178 5.956222014
Ravsurt 5.31    9.271879523 0.842975207 12.71456674
Ravsurt 5.47    9.796946179 0.692307692 4.772145446
Ravsurt 5.33    12.27335664 0.852173913 5.802874149
Eskebjerg   5.6 5.866279787 0.805309735 10.41055981
Eskebjerg   5.78    11.59095638 0.961538462 6.981631906
Eskebjerg   5.34    0.381918387 0.789473684 8.218044532
Eskebjerg   5.52    7.376130558 0.942857143 4.018040528

拟合模型(使用显式 data 参数始终是个好主意 - 如果您要对新数据使用 predict,除其他外,这是必要的)

library(nlme)
model <- lme(Rodlangde~Mehlich,random=~1|Lokalitet,data=dd)

创建预测框。 length=51 当模型是线性的(可能只是长度 = 2)时有点矫枉过正,但对非线性或广义线性模型很有用...

pframe <- with(dd,
               expand.grid(Lokalitet=levels(Lokalitet),
                      Mehlich=seq(min(Mehlich),max(Mehlich),length=51)))

预测在"level 1",即地方级别(而不是级别=0,人口级别):

pframe$Rodlangde <- predict(model,newdata=pframe,level=1)

lattice::xyplot作图:

library("lattice")
xyplot(Rodlangde~Mehlich,group=Lokalitet,data=pframe,type="l")

ggplot2:

library("ggplot2")
ggplot(dd,aes(Mehlich,Rodlangde,colour=Lokalitet))+
     geom_point()+
     geom_line(data=pframe)