混合效应MLE2模型R

jgwigjjp  于 2023-07-31  发布在  其他
关注(0)|答案(1)|浏览(112)

我需要拟合一个Ricker模型来拟合我的数据，该模型对每个区域都有随机效应，然后绘制输出。我使用How do I plot a mle2 fit of a model in ggplot2, along with the data?的答案来编写自己的模型：

dat_HG_Ricker=data[,c("EST.1.100m","Redds.100m","Clean")]

rickerfun <- function(x,a,b) {
    a*x*exp(-b*x)
}

m1 <- mle2(EST.1.100m ~ dlnorm(meanlog=log(rickerfun(Redds.100m,a,b)),
                         sdlog=exp(logsdlog)),
       method="L-BFGS-B",
       lower=c(a=1e-2,b=1e-2,logsdlog=-10),
       start=list(a=55.80234,b=0.8058173,logsdlog=0),
       data=dat_HG_Ricker)

slnorm <- function(meanlog, sdlog) {
   list(title="Log-normal",
        median=exp(meanlog),
        mean=exp(meanlog+sdlog^2/2))
}

m2 <- lm(log(EST.1.100m) ~ Redds.100m+ offset(log(Redds.100m)),
          data=dat_HG_Ricker)

m3 <- glm(EST.1.100m ~ Redds.100m+ offset(log(Redds.100m)),
      data=dat_HG_Ricker,
      family=Gamma(link="log"))

m4 <- nls(EST.1.100m ~ a*Redds.100m*exp(-b*Redds.100m),
      start=rsv,
      data=dat_HG_Ricker)

predframe <- data.frame(Redds.100m=seq(0,5.5,length=51))
predframe$mle2 <- predict(m1,newdata=predframe)
predframe$mle_med <- predict(m1,newdata=predframe,location="median")
predframe$loglm <- exp(predict(m2,newdata=predframe))
predframe$glm <- predict(m3,newdata=predframe,type="response")
predframe$nls <- predict(m4,newdata=predframe)

predframe_m <- reshape2::melt(predframe,id.var="Redds.100m",
                              variable.name="model",
                              value.name="EST.1.100m")

library(ggplot2)
ggplot(dat_HG_Ricker,aes(Redds.100m,EST.1.100m))+ geom_point() +
    geom_smooth(method="glm",
                formula=y~x + offset(log(x)),
                method.args=list(family=Gamma(link="log")))+
    geom_point(data=predframe_m,aes(colour=model,shape=model))

字符串
在代码中

rsv
$a
[1] 55.80234

$b
[1] 0.8058173

型
因此，模型的随机效应将针对变量“清洁”。我不想使用mle2（）。理想情况下，我希望图看起来像这样，每个区域（“干净”）都有自己的线来显示拟合，另外还有一条模型平均线：

的数据
数据来源：

dput(dat_HG_Ricker)
structure(list(EST.1.100m = c(50.03436426, 29.6619718333333, 
21.3333333333333, 17.5644444466667, 10.9090909066667, 1.33333333333333, 
3.926589776, 9.74298464, 12.6666666666667, 25.6666666666667, 
18.73469388, 16.11965812, 28.6349206333333, 58.4074074066667, 
36.7659574466667, 30.1363636333333, 32.37777778, 10, 13.93939394, 
32.6210045666667, 46.7801418466667, 44.12658228, 36.9047619066667, 
37.5257731933333, 19.53488372, 23.8095238066667, 4.8, 1.10168078533333, 
12.6, 50.08695652, 7.33333333333333, 7.33333333333333, 15.05376344, 
29.7721519, 19.25, 25.2470588266667, 29.5909090933333, 21.09722222, 
26.57777778, 10.98039216, 30.68992248, 22.2189054733333, 24.9263157866667, 
16.4833333333333, 19.0188679266667, 37.5, 20.73429952, 39.79288026, 
20.1165048533333, 30.3003663, 21.81818182, 21.3918128666667, 
19.6756756733333, 20.6666666666667, 2.836471968, 13.0888888866667, 
12.6666666666667, 8.15514054, 10, 4.41666666666667, 6.96, 8.5106383, 
21.3333333333333, 32.4133333333333, 14.875, 10.7916666666667, 
8.08080808, 10.2702702733333, 25.14285714, 7.56989247333333, 
14.9787234066667, 15.44827586, 19.1020408133333, 7.55555555333333, 
7.33333333333333, 7.98039216, 18.60465116, 16.1333333333333, 
2.66666666666667, 6.66666666666667, 4.66666666666667, 5.33333333333333, 
24, 20.6666666666667, 30.4411764733333, 23.2962962933333, 63.05084746, 
20.6666666666667, 58.5347985333333, 49.18996416, 35.4430379733333, 
48.32982456, 57.1358024666667, 26.1333333333333, 59.9925925933333, 
77.0833333333333, 65.3195876266667, 62.58064516, 52.5821596266667, 
105, 32.3037974666667, 17.90804598, 25.31073446, 9.52, 24.92307692, 
9.33333333333333, 20.3333333333333, 13, 11, 27.90277778, 15.3333333333333, 
25.1700680266667, 14.31884058, 16.6666666666667, 16, 41.95767196, 
24.5333333333333, 13.8461538466667, 13.0666666666667, 15.3333333333333, 
11.3333333333333, 20.5050505066667, 18.75, 9.33333333333333, 
14.71264368, 83.8666666666667, 50.9390681, 27.2333333333333, 
50.30420712, 57.11442786, 75.8137254666667, 46.8841607533333, 
64.5882352933333, 3.401360544, 62.84126984, 12.4, 57.92063492, 
1.50617283933333, 57.7073170733333, 5.33333333333333, 20.5333333333333, 
25.7627118666667, 23.4347826066667, 19.7727272733333, 16.1, 27.1515151533333, 
45.6038647333333, 3.33333333333333, 0.666666666666667, 5.33333333333333, 
6.588235294, 33.25170068, 32.3076923066667, 3.33333333333333, 
4.24036956666667, 6, 5.30975916533333, 9.62962962666667, 4.66666666666667, 
36.8609271533333, 25.3255813933333, 30.0487804866667, 4.925925926, 
4.48697225266667, 2.66666666666667, 10, 8.04597701333333, 11.3333333333333, 
10, 12, 2.78260869533333, 2.66666666666667, 8.28571428666667, 
21.2105263133333, 31.1538461533333, 13.2173913066667, 7.43589743333333, 
14.6666666666667, 15.3333333333333, 20.8333333333333, 25.1234567933333, 
25.81560284, 15.3333333333333, 19.54285714, 10.7466666666667, 
9.33333333333333, 8, 16, 19.9689922466667, 10, 15.12, 2.19047619066667, 
12.0888888866667, 19.3333333333333, 2.823529412, 13.3333333333333, 
9.16594128, 19.8333333333333, 22.6285714266667, 10.9090909066667, 
16.2095238066667, 29.2307692333333, 15.5833333333333, 14.98245614, 
40.18079096, 36.8547008533333, 31.9444444466667, 9.89333333333333, 
4.336925954, 69.5338346, 13.87755102, 14.9333333333333, 26.08955224, 
47.37078652, 22.8194444466667, 25.5163398666667, 3.43434343466667, 
11.6097561, 15.5555555533333, 10.4347826066667, 18.12121212, 
17.3333333333333, 16.7272727266667, 51.24919094, 19.1489361733333, 
27.9545454533333, 9.82222222, 13.5714285733333, 9.62962962666667, 
6.02898550733333), Redds.100m = c(0.328497125650151, 0.342184505885574, 
0.437996167533534, 0.218998083766767, 0.0273747604708459, 0.150561182589652, 
0.191623323295921, 0.0821242814125376, 0.177935943060498, 0.410621407062688, 
0.46537092800438, 0.437996167533534, 0.301122365179305, 0.164248562825075, 
0.164248562825075, 0.150561182589652, 0.205310703531344, 0.506433068710649, 
0.547495209416918, 0.287434984943882, 0.314809745414728, 0.109499041883384, 
0.123186422118806, 0.0821242814125376, 0.109499041883384, 0.287434984943882, 
0.0547495209416918, 0.218998083766767, 0.218998083766767, 0.0547495209416918, 
0.191623323295921, 0.396934026827265, 3.3997441052824, 2.48583439956132, 
3.56424785231219, 2.86967647596417, 3.3997441052824, 4.73405227563517, 
3.01590202887955, 2.77828550539207, 3.07073661122281, 2.90623286419302, 
3.25351855236703, 3.41802229939682, 3.19868397002376, 2.44927801133248, 
2.12027051727289, 1.66331566441236, 1.55364649972583, 1.53536830561141, 
1.37086455858161, 1.40742094681046, 0.73112776457686, 0.511789435203802, 
0.639736794004752, 0.603180405775909, 0.44018643190057, 0.349559813568099, 
0.401346452615225, 0.854479544277576, 0.414293112377007, 1.10046607975142, 
0.919212843086484, 0.595546349041947, 1.17814603832211, 0.699119627136199, 
1.03573278094252, 0.867426204039358, 0.841532884515795, 0.80269290523045, 
0.97099948213361, 0.893319523562921, 0.932159502848265, 0.737959606421543, 
0.168306576903159, 0.349559813568099, 0.245986535473848, 0.479026411185914, 
0.168306576903159, 0.297773174520974, 0.323666494044537, 0.103573278094252, 
0.257653834495302, 0.257653834495302, 0.121248863291907, 0.242497726583813, 
0.409214913610185, 0.36374658987572, 0.621400424371022, 0.500151561079115, 
0.682024856016975, 0.409214913610185, 0.560775992725068, 0.439527129433162, 
0.575932100636557, 0.36374658987572, 0.469839345256138, 0.424371021521673, 
0.212185510760837, 0.45468323734465, 0.670933113129645, 0.867052023121387, 
1.75474814203138, 1.96118909991742, 1.45540875309661, 1.3934764657308, 
1.36251032204789, 1.05284888521883, 0.949628406275805, 1.20767960363336, 
0.774153592072667, 1.4037985136251, 1.4037985136251, 0.949628406275805, 
1.03220478943022, 0.609000825763832, 0.815441783649876, 0.433526011560694, 
0.990916597853014, 1.20767960363336, 0.712221304706854, 0.639966969446738, 
0.660611065235343, 0.268373245251858, 0.309661436829067, 2.33064014916097, 
2.67246737103791, 2.08203853325047, 3.07644499689248, 1.61591050341827, 
1.52268489745183, 1.55376009944065, 1.95773772529521, 0.74580484773151, 
1.39838408949658, 0.52827843380982, 0.435052827843381, 0.652579241765072, 
0.870105655686762, 0.279676817899316, 1.2119328775637, 1.08763206960845, 
0.932256059664388, 1.55376009944065, 1.24300807955252, 0.932256059664388, 
0.683654443753884, 0.310752019888129, 0.341827221876942, 0.652579241765072, 
1.33623368551896, 1.6780609073959, 1.61591050341827, 1.6780609073959, 
1.95773772529521, 0.994406463642014, 1.24300807955252, 1.2119328775637, 
0.652579241765072, 0.510783200908059, 0.638479001135074, 0.297956867196368, 
0.312145289443814, 0.368898978433598, 0.482406356413167, 0.411464245175936, 
0.297956867196368, 0.297956867196368, 0.227014755959137, 0.439841089670829, 
0.411464245175936, 0.624290578887628, 0.482406356413167, 0.454029511918275, 
0.92485549132948, 0.130057803468208, 0.968208092485549, 2.03757225433526, 
0.852601156069364, 1.64739884393064, 1.76300578034682, 1.22832369942197, 
1.35838150289017, 1.48843930635838, 1.63294797687861, 1.77745664739884, 
0.852601156069364, 0.867052023121387, 0.520231213872832, 0.563583815028902, 
0.794797687861272, 0.794797687861272, 0.679190751445087, 0.505780346820809, 
0.260115606936416, 0.823699421965318, 0.346820809248555, 0.173410404624277, 
0.534682080924855, 0.317919075144509, 1.11612443600095, 1.13195598828465, 
0.941977360880234, 0.862819599461727, 0.0949893137022085, 0.253304836539223, 
0.22164173197182, 0.277052164964775, 0.134568194411462, 0.316631045674028, 
0.387873030950685, 0.569935882213251, 0.609514762922505, 0.5620201060714, 
0.269136388822924, 0.324546821815879, 0.308715269532178, 0.44328346394364, 
0.213725955829969, 0.284967941106626, 0.34037837409958, 0.245389060397372, 
0.229557508113671, 0.269136388822924, 0.245389060397372, 0.213725955829969, 
0.229557508113671, 0.284967941106626, 0.205810179688118, 0.253304836539223
), Clean = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 
9L, 9L, 9L, 9L, 9L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 
11L, 11L, 11L, 11L, 11L, 11L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 
13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 
13L, 13L, 13L, 13L, 13L, 13L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 
16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 
16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L), .Label = c("Big CR", 
"Coal CR", "Elk CR", "Goat CR", "Granite CR", "Lion CR", "Lower Goat CR", 
"Middle Squeezer CR", "Morrison CR", "NF Coal CR", "Ole CR", 
"SF Coal CR", "Squeezer CR", "Upper Goat CR", "Upper Squeezer CR", 
"Whale CR", "Wounded Buck CR"), class = "factor")), class = "data.frame", row.names = c(11L, 
12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 
25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 37L, 38L, 
39L, 40L, 42L, 43L, 44L, 96L, 97L, 101L, 102L, 104L, 105L, 108L, 
109L, 110L, 112L, 113L, 114L, 115L, 116L, 117L, 119L, 120L, 121L, 
122L, 123L, 125L, 126L, 128L, 129L, 137L, 138L, 139L, 142L, 143L, 
145L, 146L, 147L, 149L, 150L, 151L, 152L, 153L, 154L, 155L, 156L, 
157L, 158L, 159L, 160L, 161L, 162L, 164L, 165L, 167L, 169L, 195L, 
196L, 197L, 198L, 200L, 201L, 202L, 203L, 204L, 205L, 206L, 207L, 
208L, 209L, 210L, 211L, 213L, 214L, 221L, 222L, 228L, 230L, 231L, 
232L, 234L, 235L, 236L, 237L, 238L, 239L, 240L, 241L, 242L, 243L, 
245L, 246L, 247L, 248L, 249L, 250L, 251L, 252L, 254L, 263L, 265L, 
266L, 268L, 269L, 270L, 271L, 272L, 273L, 274L, 275L, 276L, 277L, 
278L, 279L, 280L, 281L, 282L, 283L, 284L, 285L, 286L, 287L, 289L, 
290L, 293L, 296L, 297L, 298L, 299L, 300L, 302L, 303L, 304L, 352L, 
353L, 355L, 364L, 365L, 366L, 367L, 368L, 369L, 371L, 372L, 373L, 
374L, 375L, 376L, 434L, 435L, 436L, 445L, 446L, 448L, 449L, 450L, 
451L, 452L, 453L, 454L, 455L, 456L, 457L, 458L, 459L, 460L, 461L, 
462L, 463L, 464L, 465L, 466L, 468L, 469L, 481L, 485L, 487L, 488L, 
490L, 492L, 493L, 494L, 495L, 496L, 497L, 498L, 499L, 500L, 501L, 
502L, 503L, 504L, 505L, 506L, 507L, 508L, 510L, 511L, 513L, 514L, 
515L, 516L, 517L, 519L))

型

r

1条答案

按热度按时间
nwnhqdif

nwnhqdif1#

不幸的是 * 仍然 * 的情况下，我知道没有简单的，现成的，* 确定性 / 频率 * 的方式来适应任意广义非线性混合模型在R。您似乎不需要一个 * 广义 * NLMM（即，您对对数正态响应很满意，它可以用对数转换的高斯响应建模）。（如果你这样做了，你的选择似乎是（1）brms（非确定性/频率论）;（2）TMB（非现成），（3）JAGS/Nimble/Greta（既不是确定性/频率论，也不是现成）...）
您可以使用非线性混合模型：nlme()可以工作（lme4::nlmer()理论上可以工作，但更挑剔）。然而，你可以使用LMM，这更容易-正如我在评论中所说的，
...使用“log-Ricker”技巧，即在您选择的线性混合模型引擎中拟合模型为log(y) ~ x + offset(log(x))（log(y) = a + log(x) + b*x → y = exp(a)*x*exp(b*x)）
对于这些数据，我最终假设随机截距和斜率是不相关的（即使用||），因为否则我会得到一个奇异的拟合（相关性= -1）;理想情况下，您的完整数据集更大，您不必这样做...

library(lme4)
m1 <- lmer(log(EST.1.100m) ~ offset(log(Redds.100m)) + Redds.100m +
         (1 + Redds.100m || Clean),
     dat_HG_Ricker)

字符串
构建预测框架，一个特定于组，一个一般（手动：您也可以使用emmeans或ggeffects或...）

pframe <- with(dat_HG_Ricker,
               expand.grid(Clean = factor(unique(Clean)),
                           Redds.100m = seq(0, 5, length = 51)))
pframe0 <- unique(pframe["Redds.100m"])
pframe$EST.1.100m <- exp(predict(m1, newdata = pframe))
pframe0$EST.1.100m <- exp(predict(m1,    newdata = pframe0, re.form = NA))

型
图：

library(ggplot2); theme_set(theme_bw())
ggplot(dat_HG_Ricker, aes(Redds.100m, EST.1.100m, colour = Clean)) +
    geom_point() +
    geom_line(data=pframe) +
    geom_line(data=pframe0, colour = "black", lwd = 2)

型
x1c 0d1x的数据
出于图形的目的，您可能需要指示更多关于每个组的数据实际位于何处的信息。您可以：

  • 使用ggalt::geom_encircle(aes(fill = Clean), colour = NA, alpha = 0.08)在区域中着色
  • 将每个人口级别的预测曲线限制为该组实际覆盖的x范围（也许用虚线扩展曲线，尽管这会变得非常繁琐......）
赞(0）分享  回复(0）举报  2023-07-31
我来回答

相关问题

    查看更多

    热门标签

    更多
    JavaquerypythonNode开发语言requestUtil数据库Table后端算法LoggerMessageElementParser

    最新问答

    更多
    • 蜀ICP备2022004421号-2 NULL123 https://www.null123.com © All rights reserved
    • 本站内容来源互联网,如果侵犯您的权益请联系我们删除, 联系方式：448109455@qq.com