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 > set.seed(1001) > > ## test 1 > x <- 0:10 > y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8) > d <- data.frame(x,y) > suppressWarnings(m1 <- mle2(y~dpois(lambda=ymax/(1+x/xhalf)), + parameters=list(ymax~1,xhalf~1), + start=list(ymax=1,xhalf=1),data=d)) > > suppressWarnings(p1 <- profile(m1)) > > suppressWarnings(m2 <- mle2(y~dpois(lambda=ymax/(1+x/xhalf)), + start=list(ymax=1,xhalf=1),data=d)) > > ## should be able to omit parameters (?) or > ## have them taken from > ## test 2: > > ReedfrogSizepred <- + structure(list(TBL = as.integer(c(9, 9, 9, 12, 12, 12, 21, 21, + 21, 25, 25, 25, 37, 37, 37)), Kill = as.integer(c(0, 2, 1, 3, + 4, 5, 0, 0, 0, 0, 1, 0, 0, 0, 0))), .Names = c("TBL", "Kill"), class = "data.frame", row.names = c("1", + "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", + "14", "15")) > > VBlogist <- function(x,sizep1,sizep2,sizep3) { + exp(sizep1*(sizep3-x))/(1+exp(sizep2*sizep1*(sizep3-x))) + } > startp <- list(sizep1=0,sizep2=1,sizep3=12) > mle2(Kill~dbinom(prob=VBlogist(TBL,sizep1,sizep2,sizep3),size=10), + start=startp, + method="Nelder-Mead", + data=ReedfrogSizepred) Call: mle2(minuslogl = Kill ~ dbinom(prob = VBlogist(TBL, sizep1, sizep2, sizep3), size = 10), start = startp, method = "Nelder-Mead", data = ReedfrogSizepred) Coefficients: sizep1 sizep2 sizep3 -0.5944408 1.6799300 12.9078275 Log-likelihood: -12.15 > > ## test 3: > f <- factor(rep(1:2,each=20)) > xhalf <- c(5,10) > ymax <- 10 > x <- rep(0:19,2) > y <- rpois(40,ymax/(1+x/xhalf[f])) > d <- data.frame(x,y) > ## plot(x,y,col=as.numeric(f)) > > m3 <- mle2(y~dpois(lambda=ymax/(1+x/xhalf)), + parameters=list(xhalf~f), + start=list(ymax=1,xhalf=1),data=d) > > m4 <- mle2(y~dpois(lambda=ymax/(1+x/xhalf)), + parameters=list(ymax~f,xhalf~f), + start=list(ymax=1,xhalf=1),data=d) Warning messages: 1: In dpois(x = c(16L, 8L, 6L, 6L, 8L, 0L, 2L, 3L, 5L, 3L, 1L, 5L, : NaNs produced 2: In dpois(x = c(16L, 8L, 6L, 6L, 8L, 0L, 2L, 3L, 5L, 3L, 1L, 5L, : NaNs produced > > suppressWarnings(m5 <- mle2(y~dpois(lambda=ymax/(1+x/xhalf)), + parameters=list(ymax~f), + start=list(ymax=1,xhalf=1),data=d)) > > anova(m2,m3,m4) Likelihood Ratio Tests Model 1: m2, y~dpois(lambda=ymax/(1+x/xhalf)) Model 2: m3, y~dpois(lambda=ymax/(1+x/xhalf)): xhalf~f Model 3: m4, y~dpois(lambda=ymax/(1+x/xhalf)): ymax~f, xhalf~f Tot Df Deviance Chisq Df Pr(>Chisq) 1 2 57.208 2 3 173.004 115.7960 1 <2e-16 *** 3 4 172.415 0.5894 1 0.4427 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > anova(m2,m5,m4) Likelihood Ratio Tests Model 1: m2, y~dpois(lambda=ymax/(1+x/xhalf)) Model 2: m5, y~dpois(lambda=ymax/(1+x/xhalf)): ymax~f Model 3: m4, y~dpois(lambda=ymax/(1+x/xhalf)): ymax~f, xhalf~f Tot Df Deviance Chisq Df Pr(>Chisq) 1 2 57.208 2 3 177.101 119.8930 1 <2e-16 *** 3 4 172.415 4.6864 1 0.0304 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > AICtab(m2,m3,m4,m5) dAIC df m2 0.0 2 m3 117.8 3 m4 119.2 4 m5 121.9 3 > > GobySurvival <- + structure(list(exper = as.integer(c(1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5)), year = as.integer(c(2000, 2000, 2000, 2000, 2000, + 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, + 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, + 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, + 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, + 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, + 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, + 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, + 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, + 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, + 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, + 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, + 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, + 2000, 2000, 2000, 2000, 2000, 2000, 2001, 2001, 2001, 2001, 2001, + 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, + 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, + 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, + 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, + 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, + 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, + 2001, 2001, 2001, 2001, 2001, 2002, 2002, 2002, 2002, 2002, 2002, + 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, + 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, + 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, + 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, + 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, + 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, + 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, + 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, + 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, + 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, + 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, + 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, + 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, + 2002)), site = structure(as.integer(c(2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 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, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1)), .Label = c("backreef", "patchreef"), class = "factor"), + head = structure(as.integer(c(15, 15, 15, 15, 15, 15, 15, + 15, 15, 15, 15, 4, 4, 4, 19, 19, 24, 24, 24, 24, 24, 24, + 6, 6, 6, 6, 6, 6, 6, 6, 9, 9, 9, 10, 10, 10, 10, 10, 10, + 10, 10, 10, 10, 10, 13, 13, 13, 13, 3, 3, 3, 3, 3, 3, 3, + 3, 2, 2, 2, 2, 5, 5, 5, 5, 12, 12, 12, 12, 7, 7, 7, 11, 11, + 11, 11, 11, 11, 11, 11, 11, 14, 14, 14, 23, 23, 23, 23, 23, + 23, 23, 23, 23, 22, 22, 22, 8, 8, 8, 8, 8, 8, 8, 8, 8, 20, + 20, 20, 20, 20, 20, 20, 20, 21, 21, 21, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 17, 17, 17, 17, 17, 17, 17, 16, 16, 16, 16, 16, + 16, 16, 16, 18, 18, 18, 26, 26, 26, 55, 55, 55, 57, 57, 41, + 41, 41, 45, 45, 47, 47, 48, 48, 58, 58, 34, 34, 34, 34, 35, + 35, 35, 35, 50, 50, 50, 32, 32, 32, 25, 25, 25, 25, 25, 33, + 33, 33, 28, 28, 31, 31, 31, 36, 36, 36, 44, 44, 44, 44, 29, + 29, 29, 27, 27, 27, 40, 40, 40, 46, 46, 46, 46, 46, 39, 39, + 39, 39, 30, 30, 30, 30, 30, 51, 51, 51, 51, 51, 51, 56, 56, + 56, 56, 56, 56, 52, 52, 52, 52, 52, 52, 55, 55, 55, 53, 53, + 53, 57, 57, 57, 57, 57, 57, 35, 35, 35, 35, 35, 35, 33, 33, + 33, 33, 33, 33, 29, 29, 29, 45, 45, 45, 45, 45, 45, 38, 38, + 38, 38, 38, 38, 27, 27, 27, 27, 27, 27, 59, 59, 59, 59, 59, + 59, 54, 54, 54, 54, 54, 54, 39, 39, 39, 39, 39, 39, 42, 42, + 42, 41, 41, 41, 41, 41, 41, 49, 49, 49, 46, 46, 46, 46, 46, + 46, 47, 47, 47, 47, 47, 47, 37, 37, 37, 43, 43, 43, 43, 43, + 43, 40, 40, 40, 40, 40, 40, 48, 48, 48, 48, 48, 48, 51, 51, + 51, 45, 45, 45, 41, 41, 41, 47, 47, 47, 37, 37, 37, 49, 49, + 49, 34, 34, 34, 25, 25, 25)), .Label = c("p1", "p10", "p11", + "p12", "p13", "p14", "p15", "p16", "p17", "p18", "p19", "p2", + "p20", "p21", "p3", "p4", "p42", "p5", "p51", "p6", "p7", + "p70", "p8", "p9", "r10", "r11", "r13", "r14", "r15", "r17", + "r18", "r19", "r2", "r20", "r21", "r22", "r23", "r24", "r25", + "r26", "r27", "r28", "r29", "r3", "r30", "r33", "r34", "r35", + "r36", "r37", "r41", "r45", "r47", "r48", "r5", "r6", "r7", + "r8", "r9"), class = "factor"), density = as.integer(c(11, + 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 3, 3, 3, 2, 2, 6, + 6, 6, 6, 6, 6, 8, 8, 8, 8, 8, 8, 8, 8, 3, 3, 3, 11, 11, 11, + 11, 11, 11, 11, 11, 11, 11, 11, 4, 4, 4, 4, 8, 8, 8, 8, 8, + 8, 8, 8, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 9, + 9, 9, 9, 9, 9, 9, 9, 9, 3, 3, 3, 9, 9, 9, 9, 9, 9, 9, 9, + 9, 3, 3, 3, 9, 9, 9, 9, 9, 9, 9, 9, 9, 8, 8, 8, 8, 8, 8, + 8, 8, 3, 3, 3, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 7, + 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 2, 2, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 4, 4, + 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 5, 5, 5, 5, 5, 3, 3, + 3, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 5, 5, 5, 5, 5, 4, 4, 4, 4, 5, 5, 5, 5, 5, 11, 11, + 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 6, 6, 6, 6, 6, 6, + 3, 3, 3, 3, 3, 3, 11, 11, 11, 11, 11, 11, 6, 6, 6, 6, 6, + 6, 11, 11, 11, 11, 11, 11, 3, 3, 3, 6, 6, 6, 6, 6, 6, 11, + 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, + 11, 11, 6, 6, 6, 6, 6, 6, 11, 11, 11, 11, 11, 11, 3, 3, 3, + 11, 11, 11, 11, 11, 11, 3, 3, 3, 6, 6, 6, 6, 6, 6, 11, 11, + 11, 11, 11, 11, 3, 3, 3, 11, 11, 11, 11, 11, 11, 6, 6, 6, + 6, 6, 6, 11, 11, 11, 11, 11, 11, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3)), qual = as.integer(c(1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, + 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, + 9, 9, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, + 10, 10, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, + 14, 14, 14, 14, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, + 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, + 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 9, + 9, 9, 9, 9, 9, 11, 11, 11, 11, 11, 12, 12, 12, 12, 18, 18, + 18, 18, 18, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, + 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, + 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, + 10, 10, 10, 10, 10, 12, 12, 12, 12, 12, 12, 16, 16, 16, 16, + 16, 16, 2, 2, 2, 5, 5, 5, 8, 8, 8, 9, 9, 9, 10, 10, 10, 9, + 9, 9, 4, 4, 4, 3, 3, 3)), d1 = as.integer(c(1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 4, 1, 1, 11, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 8, 8, 4, 8, 11, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 11, 11, 1, 1, 1, 4, 4, 11, 11, 11, 4, 8, 11, 11, + 1, 1, 1, 11, 1, 1, 8, 11, 1, 1, 11, 1, 1, 1, 1, 1, 1, 1, + 11, 11, 1, 8, 11, 4, 8, 8, 8, 11, 11, 11, 11, 11, 1, 1, 8, + 1, 1, 1, 1, 1, 1, 1, 4, 8, 1, 1, 1, 1, 1, 1, 4, 11, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 11, 11, 1, 1, 1, 1, 1, 1, 8, 1, + 1, 1, 1, 1, 8, 11, 11, 1, 4, 11, 1, 1, 3, 1, 1, 1, 1, 1, + 1, 1, 4, 2, 12, 2, 12, 3, 12, 2, 12, 1, 1, 1, 1, 1, 1, 1, + 12, 1, 1, 1, 1, 1, 4, 1, 1, 1, 2, 4, 1, 1, 12, 1, 1, 1, 1, + 4, 1, 1, 12, 1, 1, 3, 8, 1, 2, 12, 1, 1, 1, 1, 1, 8, 1, 1, + 3, 3, 12, 1, 1, 2, 12, 1, 2, 4, 8, 8, 1, 2, 3, 1, 1, 1, 1, + 1, 1, 1, 3, 3, 1, 1, 3, 1, 1, 1, 1, 1, 1, 1, 1, 12, 1, 1, + 8, 1, 2, 10, 1, 1, 12, 1, 1, 3, 1, 1, 1, 1, 2, 2, 1, 4, 6, + 3, 3, 4, 1, 4, 12, 1, 1, 3, 1, 1, 3, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 3, 6, 1, 1, 1, 1, 1, 1, 1, 1, 12, 1, 1, 12, 3, 6, + 10, 1, 1, 12, 1, 1, 8, 1, 2, 12, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 12, 2, 2, 12, 1, 12, 12, 4, 4, 4, 1, 1, 2, 1, 1, + 1, 1, 1, 8, 1, 1, 2, 1, 1, 4, 1, 1, 12, 1, 1, 12, 1, 3, 12, + 2, 4, 12, 2, 10, 12, 1, 1, 8, 1, 1, 8)), d2 = as.integer(c(4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 8, 4, 4, 70, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 11, 11, 8, 11, 70, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 70, 70, 4, 4, 4, 8, 8, 70, 70, 70, + 8, 11, 70, 70, 4, 4, 4, 70, 4, 4, 11, 70, 4, 4, 70, 4, 4, + 4, 4, 4, 4, 4, 70, 70, 4, 11, 70, 8, 11, 11, 11, 70, 70, + 70, 70, 70, 4, 4, 11, 4, 4, 4, 4, 4, 4, 4, 8, 11, 4, 4, 4, + 4, 4, 4, 8, 70, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 70, 70, + 4, 4, 4, 4, 4, 4, 11, 4, 4, 4, 4, 4, 11, 70, 70, 4, 8, 70, + 2, 2, 4, 2, 2, 2, 2, 2, 2, 2, 8, 3, 70, 3, 70, 4, 70, 3, + 70, 2, 2, 2, 2, 2, 2, 2, 70, 2, 2, 2, 2, 2, 8, 2, 2, 2, 3, + 8, 2, 2, 70, 2, 2, 2, 2, 8, 2, 2, 70, 2, 2, 4, 12, 2, 3, + 70, 2, 2, 2, 2, 2, 12, 2, 2, 4, 4, 70, 2, 2, 3, 70, 2, 3, + 8, 12, 12, 2, 3, 4, 2, 2, 2, 2, 2, 2, 2, 4, 4, 2, 2, 4, 2, + 2, 2, 2, 2, 2, 2, 2, 70, 2, 2, 10, 2, 3, 12, 2, 2, 70, 2, + 2, 4, 2, 2, 2, 2, 3, 3, 2, 6, 8, 4, 4, 6, 2, 6, 70, 2, 2, + 4, 2, 2, 4, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 4, 8, 2, 2, 2, + 2, 2, 2, 2, 2, 70, 2, 2, 70, 4, 8, 12, 2, 2, 70, 2, 2, 10, + 2, 3, 70, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 70, 3, 3, 70, + 2, 70, 70, 6, 6, 6, 2, 2, 3, 2, 2, 2, 2, 2, 10, 2, 2, 3, + 2, 2, 6, 2, 2, 70, 2, 2, 70, 2, 4, 70, 3, 6, 70, 3, 12, 70, + 2, 2, 10, 2, 2, 10))), .Names = c("exper", "year", "site", + "head", "density", "qual", "d1", "d2"), 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", "87", "88", "89", "90", + "91", "92", "93", "94", "95", "96", "97", "98", "99", "100", + "101", "102", "103", "104", "105", "106", "107", "108", "109", + "110", "111", "112", "113", "114", "115", "116", "117", "118", + "119", "120", "121", "122", "123", "124", "125", "126", "127", + "128", "129", "130", "131", "132", "133", "134", "135", "136", + "137", "138", "139", "140", "141", "142", "143", "144", "145", + "146", "147", "148", "149", "150", "151", "152", "153", "154", + "155", "156", "157", "158", "159", "160", "161", "162", "163", + "164", "165", "166", "167", "168", "169", "170", "171", "172", + "173", "174", "175", "176", "177", "178", "179", "180", "181", + "182", "183", "184", "185", "186", "187", "188", "189", "190", + "191", "192", "193", "194", "195", "196", "197", "198", "199", + "200", "201", "202", "203", "204", "205", "206", "207", "208", + "209", "210", "211", "212", "213", "214", "215", "216", "217", + "218", "219", "220", "221", "222", "223", "224", "225", "226", + "227", "228", "229", "230", "231", "232", "233", "234", "235", + "236", "237", "238", "239", "240", "241", "242", "243", "244", + "245", "246", "247", "248", "249", "250", "251", "252", "253", + "254", "255", "256", "257", "258", "259", "260", "261", "262", + "263", "264", "265", "266", "267", "268", "269", "270", "271", + "272", "273", "274", "275", "276", "277", "278", "279", "280", + "281", "282", "283", "284", "285", "286", "287", "288", "289", + "290", "291", "292", "293", "294", "295", "296", "297", "298", + "299", "300", "301", "302", "303", "304", "305", "306", "307", + "308", "309", "310", "311", "312", "313", "314", "315", "316", + "317", "318", "319", "320", "321", "322", "323", "324", "325", + "326", "327", "328", "329", "330", "331", "332", "333", "334", + "335", "336", "337", "338", "339", "340", "341", "342", "343", + "344", "345", "346", "347", "348", "349", "350", "351", "352", + "353", "354", "355", "356", "357", "358", "359", "360", "361", + "362", "363", "364", "365", "366", "367", "368", "369")) > > dicweib <- function(x,shape,scale,log=FALSE) { + if (is.matrix(x)) { + day1 <- x[,1] + day2 <- x[,2] + } else { + day1 <- x[1] + day2 <- x[2] + } + v <- log(pweibull(day2,shape,scale)-pweibull(day1,shape,scale)) + if (log) v else exp(v) + } > > GS2 <- transform(GobySurvival, + day1 = d1-1, + day2 = ifelse(d2==70,Inf,d2-1), + fexper=factor(exper)) > totmeansurv <- with(GS2,mean((d1+d2)/2)) > > mle2(cbind(day1,day2)~dicweib(exp(shape),exp(scale)), + parameters=list(scale~fexper+qual*density), + start=list(scale=log(totmeansurv),shape=0),data=GS2) Call: mle2(minuslogl = cbind(day1, day2) ~ dicweib(exp(shape), exp(scale)), start = list(scale = log(totmeansurv), shape = 0), data = GS2, parameters = list(scale ~ fexper + qual * density)) Coefficients: scale.(Intercept) scale.fexper2 scale.fexper3 scale.fexper4 1.950601011 -1.070739935 -0.767760213 -0.131513595 scale.fexper5 scale.qual scale.density scale.qual:density 0.004852567 -0.013727672 -0.219867981 0.012638159 shape -1.001618792 Log-likelihood: -443.06 There were 14 warnings (use warnings() to see them) > > ## GH 8 > set.seed(1001) > lymax <- c(0,2) > lhalf <- 0 > x <- sort(runif(200)) > g <- factor(sample(c("a","b"),200,replace=TRUE)) > y <- rnbinom(200,mu=exp(lymax[g])/(1+x/exp(lhalf)),size=2) > d2 <- data.frame(x,g,y) > fit3b <- mle2(y~dnbinom(mu=exp(lymax)/(1+x/exp(lhalf)),size=exp(logk)), + parameters=list(lhalf~1,lymax~g),data=d2, + start=list(lhalf=0,lymax=0,logk=0)) > coef(fit3b) lhalf lymax.(Intercept) lymax.gb logk -0.03968142 -0.10104397 2.12733072 0.40478459 > stopifnot(all.equal(names(coef(fit3b)), + c("lhalf", "lymax.(Intercept)", "lymax.gb", "logk"))) > > proc.time() user system elapsed 2.09 0.23 2.31