tests/mortanal.R

In bbolker/bbmle: Tools for General Maximum Likelihood Estimation

library(bbmle)

## 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)
Bm2w <- mle2(lifespan~dweibull(scale=exp(llambda),shape=alpha),
             start=list(llambda=mlife,alpha=1),
             parameters=list(llambda~group,alpha~group),
             data=x)           
Bm3w <- mle2(lifespan~dweibull(scale=exp(llambda),shape=alpha),
           start=list(llambda=mlife,alpha=3),
           parameters=list(alpha~group),
           data=x)
anova(Bm0w,Bm1w)
anova(Bm0w,Bm1w,Bm2w)
anova(Bm0w,Bm3w,Bm2w)
AICctab(Bm0w,Bm1w,Bm2w,Bm3w,sort=TRUE,nobs=nrow(x),delta=TRUE)