R Under development (unstable) (2023-12-07 r85661 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(bbmle) Loading required package: stats4 > > ## goby data in dump format > > x <- structure(list(indiv = structure(as.integer(c(20, 77, 79, 21, + 33, 40, 11, 28, 43, 85, 56, 49, 29, 37, 57, 36, 66, 65, 19, 69, + 47, 60, 23, 25, 39, 84, 12, 5, 76, 55, 32, 10, 75, 4, 78, 80, + 86, 48, 54, 22, 18, 61, 41, 74, 68, 14, 53, 45, 30, 17, 62, 3, + 7, 50, 34, 82, 8, 70, 38, 52, 2, 63, 81, 15, 44, 58, 13, 26, + 73, 83, 59, 42, 72, 67, 35, 16, 1, 46, 27, 64, 51, 24, 71, 6, + 9, 31)), .Label = c("f10al1", "f10al2", "f10al3", "f10r1", "f10r2", + "f11al1", "f11al2", "f11al3", "f11al4", "f11r1", "f11r2", "f11r3", + "f12al1", "f12al2", "f12al3", "f12al4", "f12al5", "f12r1", "f12r2", + "f12r3", "f12r4", "f12r5", "f12r6", "f13al1", "f13r1", "f14al1", + "f14al2", "f14r1", "f14r2", "f15al1", "f15al2", "f15r1", "f15r2", + "f18al1", "f18al2", "f18r1", "f18r2", "f19al1", "f19r1", "f19r2", + "f1al1", "f1al2", "f1r1", "f20al1", "f20al2", "f20al3", "f20r1", + "f20r2", "f20r3", "f2al1", "f2al2", "f2al3", "f2al4", "f2r1", + "f2r2", "f2r3", "f2r4", "f3al1", "f3al2", "f3r1", "f3r2", "f4al1", + "f5al1", "f5al2", "f5r1", "f5r2", "f6al1", "f6al2", "f6r1", "f7al1", + "f7al2", "f7al3", "f7al4", "f7al5", "f7r1", "f7r2", "f7r3", "f7r4", + "f7r5", "f7r6", "f9al1", "f9al2", "f9al4", "f9r1", "f9r2", "f9r3" + ), class = "factor"), group = structure(as.integer(c(5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)), .Label = c("AL", + "AL-Rat5th", "AL-RatOv", "R", "R-ALat5th"), class = "factor"), + lifespan = as.integer(c(391, 370, 346, 341, 334, 320, 319, + 317, 314, 307, 295, 260, 30, 10, 397, 380, 364, 355, 352, + 341, 340, 339, 336, 320, 314, 312, 308, 302, 296, 290, 284, + 267, 263, 263, 255, 253, 242, 222, 220, 181, 64, 36, 192, + 192, 189, 186, 183, 181, 180, 176, 173, 171, 170, 169, 166, + 11, 247, 235, 234, 233, 232, 224, 221, 220, 215, 210, 210, + 204, 202, 17, 13, 301, 300, 296, 281, 271, 253, 250, 241, + 239, 232, 221, 220, 214, 33, 30))), .Names = c("indiv", "group", + "lifespan"), class = "data.frame", row.names = c("1", "2", "3", + "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", + "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", + "27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", + "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", "48", + "49", "50", "51", "52", "53", "54", "55", "56", "57", "58", "59", + "60", "61", "62", "63", "64", "65", "66", "67", "68", "69", "70", + "71", "72", "73", "74", "75", "76", "77", "78", "79", "80", "81", + "82", "83", "84", "85", "86")) > > mlife <- log(mean(x$lifespan)) > Bm0w <- mle2(lifespan~dweibull(scale=exp(llambda),shape=alpha), + start=list(llambda=mlife,alpha=1), + data=x) > Bm1w <- mle2(lifespan~dweibull(scale=exp(llambda),shape=alpha), + start=list(llambda=mlife,alpha=1), + parameters=list(llambda~group), + data=x) Warning message: In dweibull(x = c(391L, 370L, 346L, 341L, 334L, 320L, 319L, 317L, : NaNs produced > Bm2w <- mle2(lifespan~dweibull(scale=exp(llambda),shape=alpha), + start=list(llambda=mlife,alpha=1), + parameters=list(llambda~group,alpha~group), + data=x) Warning messages: 1: In dweibull(x = c(391L, 370L, 346L, 341L, 334L, 320L, 319L, 317L, : NaNs produced 2: In dweibull(x = c(391L, 370L, 346L, 341L, 334L, 320L, 319L, 317L, : NaNs produced > Bm3w <- mle2(lifespan~dweibull(scale=exp(llambda),shape=alpha), + start=list(llambda=mlife,alpha=3), + parameters=list(alpha~group), + data=x) Warning messages: 1: In dweibull(x = c(391L, 370L, 346L, 341L, 334L, 320L, 319L, 317L, : NaNs produced 2: In dweibull(x = c(391L, 370L, 346L, 341L, 334L, 320L, 319L, 317L, : NaNs produced 3: In dweibull(x = c(391L, 370L, 346L, 341L, 334L, 320L, 319L, 317L, : NaNs produced 4: In dweibull(x = c(391L, 370L, 346L, 341L, 334L, 320L, 319L, 317L, : NaNs produced > anova(Bm0w,Bm1w) Likelihood Ratio Tests Model 1: Bm0w, lifespan~dweibull(scale=exp(llambda),shape=alpha) Model 2: Bm1w, lifespan~dweibull(scale=exp(llambda),shape=alpha): llambda~group Tot Df Deviance Chisq Df Pr(>Chisq) 1 2 1043.5 2 6 1015.5 27.945 4 1.28e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > anova(Bm0w,Bm1w,Bm2w) Likelihood Ratio Tests Model 1: Bm0w, lifespan~dweibull(scale=exp(llambda),shape=alpha) Model 2: Bm1w, lifespan~dweibull(scale=exp(llambda),shape=alpha): llambda~group Model 3: Bm2w, lifespan~dweibull(scale=exp(llambda),shape=alpha): llambda~group, alpha~group Tot Df Deviance Chisq Df Pr(>Chisq) 1 2 1043.5 2 6 1015.5 27.945 4 1.28e-05 *** 3 10 1008.8 6.736 4 0.1505 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > anova(Bm0w,Bm3w,Bm2w) Likelihood Ratio Tests Model 1: Bm0w, lifespan~dweibull(scale=exp(llambda),shape=alpha) Model 2: Bm3w, lifespan~dweibull(scale=exp(llambda),shape=alpha): alpha~group Model 3: Bm2w, lifespan~dweibull(scale=exp(llambda),shape=alpha): llambda~group, alpha~group Tot Df Deviance Chisq Df Pr(>Chisq) 1 2 1043.5 2 6 1038.5 4.9434 4 0.2932 3 10 1008.8 29.7377 4 5.535e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > AICctab(Bm0w,Bm1w,Bm2w,Bm3w,sort=TRUE,nobs=nrow(x),delta=TRUE) dAICc df Bm1w 0.0 6 Bm2w 3.1 10 Bm0w 19.0 2 Bm3w 23.0 6 > > > proc.time() user system elapsed 1.46 0.14 1.59