https://github.com/cran/bbmle
Tip revision: 13982dbec4de7e8f454e60e5769fa597e6de56aa authored by Ben Bolker on 11 May 2022, 07:20:02 UTC
version 1.0.25
version 1.0.25
Tip revision: 13982db
formulatest.Rout.save
R Under development (unstable) (2019-06-19 r76722) -- "Unsuffered Consequences"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
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.
Natural language support but running in an English locale
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.488 0.132 2.757
