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

Question

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

Rodlangde 表示根长 Mehlich 是指植物有效磷 Lokalitet 是地方的意思

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

我非常想制作一个图，您可以在其中看到 10 个（我有 10 个位置）不同的线性图，这些线性图具有相同的斜率但不同的截距，该模型由该模型组成。我现在尝试在网上搜索at solution 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

Answer 1

拟合模型（使用显式 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)))

预测“级别 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)