exploring the properties of truncated distributions with 'simple' data

trunc_phenom.R

library(readxl)

my_sum <- function(f) {
    C <- contr.sum(length(levels(f)))
    colnames(C) <- paste0(".sum_",levels(f)[-length(levels(f))])
    return(C)
}

dd <- read_excel("Truncated.compois.problem.xlsx")
dd$pop <- as.factor(dd$pop)
contrasts(dd$pop) <- my_sum(dd$pop)
dd$morph <- as.factor(dd$morph)
contrasts(dd$morph) <- my_sum(dd$morph)

## options(contrasts = c("contr.sum","contr.poly"))

library(glmmTMB)
trunc_compois <- glmmTMB(mates ~ pop*morph,
                         family='truncated_compois', data=dd)
summary(trunc_compois)

pois <- update(trunc_compois, family = poisson)

## diagnose this: simpler
trunc_pois <- update(trunc_compois, family = truncated_poisson)

with(dd, table(mates,pop, morph))

dd_ELS <- subset(dd, pop == "ELS")
trunc_pois_2 <- update(trunc_pois,
                       formula = . ~ morph,
                       data = dd_ELS)
summary(trunc_pois_2)
with(dd_ELS, table(morph, mates))

dd_ELS_L <- subset(dd_ELS, morph == "L")
dd_ELS_S <- subset(dd_ELS, morph == "S")

trunc_pois_3 <- update(trunc_pois_2,
                       formula = ~ 1,
                       data = dd_ELS_L)
## false convergence warning
summary(trunc_pois_3)
with(dd_ELS_L, table(mates))     

## compare with bbmle
library(bbmle)
dtruncatedpois <- function(x, lambda, log) {
    res <- dpois(x, lambda, log = log)
    res[x==0] <- -Inf
    res <- res - ppois(0, lambda, lower.tail = FALSE, log.p = TRUE)
    if (log) res else exp(res)
}
m1 <- mle2(mates ~ dtruncatedpois(exp(log_lambda)),
           start = list(log_lambda = 0),
           data = dd_ELS_L)

nllfun <- function(log_lambda) {
    -sum(dtruncatedpois(dd$mates, exp(log_lambda), log = TRUE))
}

pvec <- seq(-30, 0, length=101)
lvec <- sapply(pvec, nllfun)

plot(pvec, lvec)

## truncated poisson of x ==1
tpl <- function(lambda) (log(lambda) - lambda - log(1-exp(-lambda)))
tpl2 <- function(lambda) (log(lambda) - lambda -
                          ppois(0, lambda, lower.tail = FALSE, log.p = TRUE))

## for lambda << 1, 1-exp(-lambda) → lambda
## negative log likelihood → lambda

tpl(2)
dtruncatedpois(1, 2, TRUE)

D(expression(log(lambda) - lambda - log(1-exp(-lambda))), "lambda")
## 1/lambda - 1 - 1/(exp(lambda) - 1)
## == 0

curve(-1*tpl(x), from = 1e-9, 1, log = "xy")
curve(-1*tpl2(x), from = 1e-12, 1, log = "xy")