Search code examples
rggplot2mixed-modelsmle

mixed effects mle2 model r


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 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: enter image description here

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))

Solution

  • 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*xy = 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)
    

    enter image description here

    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 ...)