R/099_explore-weibull.R

In robschick/tangled: Creates data objects and plots for Knowlton, Schick et al. manuscript on entanglement

library(dplyr)
# explore weibull vs exponential
pregnant <- read_csv(here::here('data-raw', '2021-02-25_years_since_pregnancy_1980-2013.csv'))%>% 
  dplyr::select(!starts_with('X')) %>% 
  filter(Year >=1980 & Year <= 2013) %>% 
  drop_na(Pregnant) %>% 
  drop_na(elapsed) %>% 
  filter(elapsed >= 1)

pregnant$Pregnant[pregnant$Pregnant == 'F'] <- '1'
pregnant$Pregnant <- as.numeric(pregnant$Pregnant)

# assemble a dataset like firedata
# contains the unique intervals and the number of unique intervals
pregdata <- pregnant %>% 
  group_by(elapsed) %>% 
  tally()
colnames(pregdata) <- c('years', 'npreg')
allyears <- seq(0, max(pregdata$years), by = 1)
npregyr <- allyears * 0
for(i in 1:length(pregdata$years)){
  npregyr[allyears[pregdata$years[i] == allyears]] <- pregdata$npreg[i]
}
plot(allyears, npregyr, type = 's')

wlik <- function(param){
  lp <- param[1]
  cp <- param[2]
  lik <- log(cp * lp^cp * pregnant$elapsed^(cp - 1)) - (lp * pregnant$elapsed)^cp
  return(-sum(lik))
}

param0 <- c(0.1, 1)
out <- optim(param0, wlik, lower = c(0.01, 0.01), upper = c(10, 10), method = "L-BFGS-B")

# plot weibull density
lw <- out$par[1]
cw <- out$par[2]
time <- seq(0, 40, length = 100)
hazard <- cw*lw^cw*time^(cw-1)
fweib <- hazard*exp(-(lw*time)^cw)
plot(allyears, npregyr / sum(npregyr), type = 's')
lines(time, fweib)
lines(time, hazard, lty = 2)


# LRT
lmle <- 1 / mean(pregnant$elapsed)
like <- -sum(log(lmle) - lmle*pregnant$elapsed)

likw <- out$value
dev <- 2 * (like - likw)
pr <- 1 - pchisq(dev, 1) # difference in number of parameters, e.g. 2-1

## Adding in Covariates
pregnant <- readRDS("/Users/rob/Documents/research/projects/right-whales_PCOD/rightWhaleEntanglement/data/2021-06-15_pregnancy-with-covariates.rds")
xdat <- pregnant[, c("Pregnant", "health_scl", "mon_below_scl", "elapsed", "severity", "num_events", "sev_num", "std_year")]
xdat <- xdat[xdat$Pregnant > 0, ]
xdat <- xdat[xdat$elapsed > 0, ]
xdat$sev_fac <- xdat$sev_num
xdat$sev_fac[xdat$sev_num > 0] <- 1

# likelihood function
wlik <- function(param){
  lp <- param[1]
  b0 <- param[2]
  b1 <- param[3]
  b2 <- param[4]
  cp <- exp(b0 + b1 * xdat$sev_fac + b2 * xdat$mon_below_scl) 
  lik <- log(cp * lp^cp * xdat$elapsed^(cp - 1)) - (lp * xdat$elapsed)^cp
  
  # if(is.infinite(sum(lik))){
  #   browser()
  # }
  return(-sum(lik))
}

param0 <- c(0.1, 1, -1, -1) 
out <- optim(param0, wlik, 
             method = "L-BFGS-B", hessian = TRUE)
out$par

sigma <- solve(out$hessian)
se <- sqrt(diag(sigma))
b1_up <- out$par[3] + 1.96*se[3]
b1_low <- out$par[3] - 1.96*se[3]

b2_up <- out$par[4] + 1.96*se[4]
b2_low <- out$par[4] - 1.96*se[4]

my_pars <- data.frame(cov = c('entanglement', 'health'), mle = c(out$par[3], out$par[4]),
                      low = c(b1_low, b2_low),
                      high = c(b1_up, b2_up))
my_pars