如何在R中绘制混合效应模型的结果?

k2fxgqgv  于 2023-11-14  发布在  其他
关注(0)|答案(1)|浏览(121)

我想使用以下数据拟合一个混合效应模型,并使用R.绘制建模数据,以显示使用该模型预测值的平均95% CI。

tempEf <- data.frame(
  N = rep(c("Nlow", "Nhigh"), each=300),
  Myc = rep(c("AM", "ECM"), each=150, times=2),
  TRTYEAR = runif(600, 1, 15),
  ID = rep(c("1","2","3","4","5"), each=10, times=12)   #5 sites
)

# Make up some response data
tempEf$r <- 2*tempEf$TRTYEAR +                   
  -8*as.numeric(tempEf$Myc=="ECM") +
  4*as.numeric(tempEf$N=="Nlow") +
  0.1*tempEf$TRTYEAR * as.numeric(tempEf$N=="Nlow") +
  0.2*tempEf$TRTYEAR*as.numeric(tempEf$Myc=="ECM") +
  -11*as.numeric(tempEf$Myc=="ECM")*as.numeric(tempEf$N=="Nlow")+ 
  0.5*tempEf$TRTYEAR*as.numeric(tempEf$Myc=="ECM")*as.numeric(tempEf$N=="Nlow")+ 
  as.numeric(tempEf$ID) +  #Random intercepts; intercepts will increase by 1
  tempEf$TRTYEAR/10*rnorm(600, mean=0, sd=2)    #Add some noise

字符串
enter image description here

6psbrbz9

6psbrbz91#

根据您提供的信息,您只能猜测您想要什么,因为您没有包含任何适合模型的代码。这里有一个建议,您可以如何打印/绘制混合模型的调整预测,希望这对您有帮助。
有几个R软件包可以帮助你计算预测/估计的边际均值/边际效应。在这个特定的例子中,我使用的是ggeffects -package。你可以在链接网站上找到很多例子和文档。

library(lme4)
#> Loading required package: Matrix
library(ggeffects)

tempEf <- data.frame(
  N = rep(c("Nlow", "Nhigh"), each = 300),
  Myc = rep(c("AM", "ECM"), each = 150, times = 2),
  TRTYEAR = runif(600, 1, 15),
  ID = rep(c("1", "2", "3", "4", "5"), each = 10, times = 12) # 5 sites
)

# Make up some response data
tempEf$r <- 2 * tempEf$TRTYEAR +
  -8 * as.numeric(tempEf$Myc == "ECM") +
  4 * as.numeric(tempEf$N == "Nlow") +
  0.1 * tempEf$TRTYEAR * as.numeric(tempEf$N == "Nlow") +
  0.2 * tempEf$TRTYEAR * as.numeric(tempEf$Myc == "ECM") +
  -11 * as.numeric(tempEf$Myc == "ECM") * as.numeric(tempEf$N == "Nlow") +
  0.5 * tempEf$TRTYEAR * as.numeric(tempEf$Myc == "ECM") * as.numeric(tempEf$N == "Nlow") +
  as.numeric(tempEf$ID) + # Random intercepts; intercepts will increase by 1
  tempEf$TRTYEAR / 10 * rnorm(600, mean = 0, sd = 2) # Add some noise

m <- lmer(r ~ TRTYEAR * Myc + N + (1 | ID), data = tempEf)

ggpredict(m, c("TRTYEAR", "Myc"))
#> # Predicted values of r
#> 
#> # Myc = AM
#> 
#> TRTYEAR | Predicted |         95% CI
#> ------------------------------------
#>       0 |      4.11 | [ 2.59,  5.63]
#>       2 |      8.29 | [ 6.82,  9.76]
#>       6 |     16.64 | [15.24, 18.05]
#>       8 |     20.82 | [19.42, 22.22]
#>      10 |     25.00 | [23.59, 26.41]
#>      16 |     37.53 | [36.01, 39.06]
#> 
#> # Myc = ECM
#> 
#> TRTYEAR | Predicted |          95% CI
#> -------------------------------------
#>       0 |     -9.24 | [-10.77, -7.71]
#>       2 |     -4.20 | [ -5.67, -2.72]
#>       6 |      5.89 | [  4.48,  7.30]
#>       8 |     10.93 | [  9.53, 12.33]
#>      10 |     15.97 | [ 14.56, 17.38]
#>      16 |     31.10 | [ 29.58, 32.63]
#> 
#> Adjusted for:
#> *  N = Nhigh
#> * ID = 0 (population-level)
#> 
#> Not all rows are shown in the output. Use `print(..., n = Inf)` to show
#>   all rows.

ggpredict(m, c("TRTYEAR", "Myc")) |> plot()

字符串
x1c 0d1x的数据
创建于2023-11-05使用reprex v2.0.2

相关问题