cpglm.R
#######################################################
## compound Poisson GLM ##
## Author: Wayne Zhang, actuary_zhang@hotmail.com ##
#######################################################
cpglm <- function(formula, link = "log", data, weights, offset,
subset, na.action = NULL, contrasts = NULL,
control = list(), chunksize = 0, ...) {
call <- match.call()
if (missing(data))
data <- environment(formula)
# use bigglm for big data sets
if (chunksize) {
ans <- cpglm.profile.bigglm(call, data)
} else {
fr <- cpglm.mf(call, contrasts)
link.power <- make.link.power(link)
control <- do.call("cpglm.control", control)
ans <- cpglm.profile(fr, link.power, control)
}
ans@formula <- formula
ans@call <- call
ans@na.action <- na.action
ans@contrasts <- contrasts
return(ans)
}
# function to implement the profile likelihood approach
# Use cpglm.profile.bigglm for big data sets
cpglm.profile <- function(fr, link.power = 0, control = list()){
# return profile loglikelihood
# log.phi indicates whether phi (as in parm[1]) is on the log scale
llik_profile <- function(parm, log.phi = TRUE){
phi <- ifelse(log.phi, exp(parm[1]), parm[1])
p <- parm[2]
fit2 <- cpglm.fit(fr, p = p, link.power)
- sum(log(dtweedie.series(fr$Y, p, fit2$fitted.values, phi)))
}
# generate starting values for phi
pstart <- 1.5
fit <- cpglm.fit(fr, pstart, link.power)
phistart <- sum((fit$weights * fit$residuals^2)) / fit$df.residual
parm <- c(log(phistart), pstart)
# optimize the profile loglikelihood
opt_ans <- nlminb(parm,llik_profile, gradient = NULL,
hessian = NULL, log.phi = TRUE,
lower = c(-Inf, control$bound.p[1]),
upper = c(Inf, control$bound.p[2]),
control = list(trace = control$trace,
iter.max = control$max.iter))
if (opt_ans$convergence) warning(opt_ans$message)
# parameter values at optima
p.max <- opt_ans$par[2]
phi.max <- exp(opt_ans$par[1])
# fit glm using the optimized index parameter
fit <- cpglm.fit(fr, p.max, link.power)
class(fit) <- "glm"
# compute vcov for p and phi
pm <- c(phi.max, p.max)
hs <- hess(pm,llik_profile, log.phi = FALSE)
dimnames(hs) <- list(c("phi","p"),c("phi","p"))
vc <- vcov(fit)
attr(vc,"phi_p") <- svd.inv(hs)
# return results
out <- new("cpglm",
coefficients = fit$coefficients, residuals = fit$residuals,
fitted.values = fit$fitted.values, weights = fit$weights,
linear.predictors = fit$linear.predictors,
df.residual = as.integer(fit$df.residual),
deviance = fit$deviance, call = call("foo"),
aic = dtweedie.nlogl(fr$Y, fitted(fit), phi.max, p.max) + 2 * fit$rank,
formula = ~ 1, control = control, contrasts = NULL,
p = p.max, phi = phi.max, iter = fit$iter, converged = fit$converged,
link.power= link.power, model.frame = fr$mf, na.action = NULL,
offset = fr$offset, prior.weights = fit$prior.weights, y = fr$Y,
inits = NULL, vcov = vc)
return(out)
}
# function to implement the profile likelihood approach for big data sets
cpglm.profile.bigglm <- function(call, data){
# restruct the call
if (is.null(control <- call$control))
control <- list()
if (is.null(link <- call$link))
link <- "log"
control <- do.call("cpglm.control", eval(control))
link.power <- make.link.power(eval(link))
mc <- match(c("link", "control", "family"), names(call), 0L)
call <- call[-c(1, mc[mc>0])]
pstart <- 1.5
call$family <- tweedie(var.power = pstart, link.power = link.power)
# default maxit in bigglm: the default seems not to be enough
call$maxit <- ifelse(is.null(call$maxit), 50, call$maxit)
Y <- eval(call$formula[[2]], envir = data)
# profile loglikelihood
llik_profile <- function(parm, log.phi = TRUE){
phi <- ifelse(log.phi, exp(parm[1]), parm[1])
p <- parm[2]
call$family <- tweedie(var.power = p, link.power = link.power)
fit2 <- do.call("bigglm", as.list(call))
- sum(log(dtweedie.series(Y, p, fitted(fit2, data)$fitted.values, phi)))
}
# generate starting values for phi
fit <- do.call("bigglm", as.list(call))
phistart <- fit$qr$ss / fit$df.resid
parm <- c(log(phistart), pstart)
# optimize the profiled loglikelihood
opt_ans <- nlminb(parm, llik_profile, gradient = NULL,
hessian = NULL, log.phi = TRUE,
lower = c(-Inf, control$bound.p[1]),
upper = c(Inf, control$bound.p[2]),
control = list(trace = control$trace,
iter.max = control$max.iter))
if (opt_ans$convergence) warning(opt_ans$message)
p.max <- opt_ans$par[2]
phi.max <- exp(opt_ans$par[1])
# fit glm using the optimized index parameter
call$family <- tweedie(var.power = p.max, link.power = link.power)
fit <- do.call("bigglm", as.list(call))
fs <- fitted(fit, data)
# compute vcov for p and phi
pm <- c(phi.max, p.max)
hs <- hess(pm, llik_profile, log.phi = FALSE)
dimnames(hs) <- list(c("phi","p"), c("phi","p"))
vc <- vcov(fit)
attr(vc,"phi_p") <- svd.inv(hs)
# return results
out <- new("cpglm", coefficients = coef(fit),
residuals = fs$residuals, fitted.values = fs$fitted.values,
linear.predictors = fs$linear.predictors,
weights = fs$weights, df.residual = as.integer(fit$df.resid),
aic = dtweedie.nlogl(Y, fs$fitted.values, phi.max, p.max) +
2 * (fit$n - fit$df.resid),
deviance = fit$deviance, call = call("foo"),
formula = ~ 1, control = control,
contrasts = NULL, p = p.max, phi = phi.max,
iter = fit$iterations, converged = fit$converged,
link.power= link.power, model.frame = as.data.frame(NULL),
na.action = NULL, offset = fs$offset,
prior.weights = fs$prior.weights, y = Y,
inits = NULL, vcov = vc)
return(out)
}