英文:
mixed effects mle2 model r
问题
我需要将一个Ricker模型拟合到具有每个地区随机效应的数据中，然后绘制输出。我使用了https://stackoverflow.com/questions/57153916/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
因此，模型的随机效应将应用于变量"Clean"。我并不固守于使用mle2()。理想情况下，我想要一个类似于这样的图，每个地区("Clean")都有自己的拟合线，并且另外有一个模型均值线:
数据:
dput(dat_HG_Ricker)
# 数据的结构
这是您的代码的翻译部分。如有需要，请提出进一步的问题。
英文:
I need to fit a ricker model to my data that has random effects for each area and then plot the output. I used the answer from https://stackoverflow.com/questions/57153916/how-do-i-plot-a-mle2-fit-of-a-model-in-ggplot2-along-with-the-data
to write my own model:
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))
In the code
rsv
$a
[1] 55.80234
$b
[1] 0.8058173
So, the random effect to the model would be for the variable "Clean". I am not married to using mle2(). Ideally, I would like a plot that looks something like this, with each area ("Clean") having its own line to show fit, and additionally a model mean line:
Data:
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))
答案1
得分: 1
抱歉，你提供的内容中包含了一些代码和特殊格式的标记，无法直接进行翻译。以下是你提供的内容中的文本部分的翻译：
"Unfortunately still the case that I am aware of no easy, off-the-shelf, deterministic/frequentist way to fit arbitrary generalized nonlinear mixed models in R. You don't appear to need a generalized NLMM (i.e., you're happy with a log-Normal response, which can be model with a log-transformed Gaussian response). (If you did, your choices seem to be (1) brms (not deterministic/frequentist); (2) TMB (not off-the-shelf), (3) JAGS/Nimble/Greta (neither deterministic/frequentist nor off-the-shelf) ...)"
"很不幸，我仍然没有找到在R中轻松、现成的、确定性/频率派的方式来拟合任意的广义非线性混合模型。你似乎不需要广义NLMM（即，你对对数正态响应感到满意，可以使用对数转换的高斯响应进行建模）。 （如果需要的话，你的选择似乎是（1）brms（不确定性/频率派）；（2）TMB（不是现成的）；（3）JAGS/Nimble/Greta（既不确定性/频率派也不是现成的）...）"
"You could use a nonlinear mixed model: nlme() would work (lme4::nlmer() would theoretically work, but is finickier). However, you can get away with a LMM, which is easier — as I said in the comments,"
"你可以使用非线性混合模型：nlme()可以工作（lme4::nlmer() 理论上 可以工作，但更加挑剔）。不过，你也可以使用线性混合模型（LMM），这更容易—就像我在评论中所说的那样，"
"> ... use the "log-Ricker" trick, i.e. fit the model as log(y) ~ x + offset(log(x)) in the linear mixed model engine of your choice (log(y) = a + log(x) + b*x → y = exp(a)*x*exp(b*x))
"> …使用“对数Ricker”技巧，即在你选择的线性混合模型引擎中将模型拟合为 log(y) ~ x + offset(log(x))（log(y) = a + log(x) + b*x → y = exp(a)*x*exp(b*x)）"
"For these data, I ended up assuming the random intercepts and slopes were uncorrelated (i.e., using ||) because otherwise I got a singular fit (correlation = -1); ideally your full data set is larger and you don't have to do this ..."
"对于这些数据，我最终假设随机截距和斜率是不相关的（即，使用 || ），否则我得到了一个奇异拟合（相关性 = -1）；理想情况下，你的完整数据集更大，不必这样做..."
library(lme4)
m1 <- lmer(log(EST.1.100m) ~ offset(log(Redds.100m)) + Redds.100m +
(1 + Redds.100m || Clean),
dat_HG_Ricker)
library(lme4)
m1 <- lmer(log(EST.1.100m) ~ offset(log(Redds.100m)) + Redds.100m +
(1 + Redds.100m || Clean),
dat_HG_Ricker)
"Construct prediction frames, one group-specific and one general (by hand: you could also use emmeans or ggeffects or ...)"
"构建预测框架，一个是特定于组的，一个是通用的（手动：你也可以使用 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))
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))
"Plot:"
"绘图："
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)
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)
"For graphical purposes, you might want to indicate some more about where the data for each group actually lie. You could:"
"为了图形效果，你可能想要更明确地指示每个组的数据实际位于哪里。你可以："
-
use
ggalt::geom_encircle(aes(fill = Clean), colour = NA, alpha = 0.08)to colour in regions -
restrict each population-level prediction curve to the x-range actually covered for that group (perhaps extending the curve with a dashed line, although that gets pretty fussy ...)"
-
使用 `gg
英文:
It is unfortunately still the case that I am aware of no easy, off-the-shelf, deterministic/frequentist way to fit arbitrary generalized nonlinear mixed models in R. You don't appear to need a generalized NLMM (i.e., you're happy with a log-Normal response, which can be model with a log-transformed Gaussian response). (If you did, your choices seem to be (1) brms (not deterministic/frequentist); (2) TMB (not off-the-shelf), (3) JAGS/Nimble/Greta (neither deterministic/frequentist nor off-the-shelf) ...)
You could use a nonlinear mixed model: nlme() would work (lme4::nlmer() would theoretically work, but is finickier). However, you can get away with a LMM, which is easier — as I said in the comments,
> ... use the "log-Ricker" trick, i.e. fit the model as log(y) ~ x + offset(log(x)) in the linear mixed model engine of your choice (log(y) = a + log(x) + b*x → y = exp(a)*x*exp(b*x))
For these data, I ended up assuming the random intercepts and slopes were uncorrelated (i.e., using ||) because otherwise I got a singular fit (correlation = -1); ideally your full data set is larger and you don't have to do this ...
library(lme4)
m1 <- lmer(log(EST.1.100m) ~ offset(log(Redds.100m)) + Redds.100m +
(1 + Redds.100m || Clean),
dat_HG_Ricker)
Construct prediction frames, one group-specific and one general (by hand: you could also use emmeans or ggeffects or ...)
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))
Plot:
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)
For graphical purposes, you might want to indicate some more about where the data for each group actually lie. You could:
-
use
ggalt::geom_encircle(aes(fill = Clean), colour = NA, alpha = 0.08)to colour in regions -
restrict each population-level prediction curve to the x-range actually covered for that group (perhaps extending the curve with a dashed line, although that gets pretty fussy ...)
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论