https://github.com/cran/cplm
Tip revision: 31a5c508df295d9e5b5c8d3ffdc2c8b27b4e819f authored by Wayne Zhang on 26 October 2011, 00:00:00 UTC
version 0.3-1
version 0.3-1
Tip revision: 31a5c50
cpglmm.R
#######################################################
## Compound Poisson GLMM ##
## Author: Wayne Zhang, actuary_zhang@hotmail.com ##
#######################################################
cpglmm <- function(formula, link="log", data, weights, offset,
subset, na.action, betastart=NULL, phistart=NULL,
pstart=NULL, contrasts = NULL, control = list()) {
call <- match.call()
if (missing(data))
data <- environment(formula)
link.power <- make.link.power(link)
# create model frame
fr <- lme4:::lmerFrames(call, formula, contrasts)
offset <- wts <- NULL
if (length(fr$wts))
wts <- fr$wts
if (length(fr$off))
offset <- fr$off
if (!is.null(betastart)){
if (length(betastart) != ncol(fr$X))
stop(gettextf("number of 'betastart' is %d should
equal %d (number of mean parameters)",
length(betastart), ncol(fr$X)), domain = NA)
}
if (!is.null(phistart) && length(phistart)>1)
stop("multiple values specified for 'phistart'")
if (!is.null(phistart) && phistart<=0)
stop("value of 'phistart' should be greater than 0")
if (!is.null(pstart) && length(pstart)>1)
stop("multiple values specified for 'pstart'")
if (!is.null(pstart) && (pstart<=1 || pstart>=2))
stop("value of 'pstart' should be between 1 and 2")
if (is.null(pstart))
pstart <- 1.5
# generate starting values
glmFit <- glm.fit(fr$X, fr$Y, weights = wts, offset = offset,
family = tweedie(var.power=pstart,link.power=link.power),
intercept = attr(attr(fr$mf, "terms"), "intercept") > 0)
if (is.null(betastart))
betastart <- glmFit$coefficients
names(betastart) <- names(fr$fixef)
if (is.null(phistart))
phistart <- sum((fr$Y-glmFit$fitted.values)^2/glmFit$fitted.values^pstart)/
glmFit$df.residual
# get factor list and dims
FL <- lme4:::lmerFactorList(formula, fr, 0L, 0L)
ctr <- do.call(cpglmm.control, control)
FL$dims["mxit"] <- ctr$max.iter
FL$dims["mxfn"] <- ctr$max.fun
dm <- lme4:::mkZt(FL, NULL)
dm$dd["verb"] <- ctr$trace
AGQlist <- .Call(lme4:::lme4_ghq, 1)
M1 <- length(levels(dm$flist[[1]]))
n <- ncol(dm$Zt)
if (M1 >= n) {
msg1 <- "Number of levels of a grouping factor for the random effects\n"
msg3 <- "n, the number of observations"
if (dm$dd["useSc"])
stop(msg1, "must be less than ", msg3)
else if (M1 == n)
message(msg1, "is *equal* to ", msg3)
}
# environment to evaluate tweedie density
tmp.env <- new.env( )
environment(ldtweedie) <- tmp.env
# input cpglmm class for optimization
ans <- new(Class = "cpglmm", env = tmp.env, nlmodel = (~I(x))[[2]],
frame = fr$mf, call = call, flist = dm$flist,
Zt = dm$Zt, X = fr$X, y = fr$Y, pWt = unname(glmFit$prior.weights),
offset = unname(fr$off), Gp = unname(dm$Gp), dims = dm$dd,
ST = dm$ST, A = dm$A, Cm = dm$Cm, Cx = (dm$A)@x, L = dm$L,
deviance = dm$dev, fixef = betastart, ranef = numeric(dm$dd[["q"]]),
u = numeric(dm$dd[["q"]]), eta = unname(glmFit$linear.predictors),
mu = unname(glmFit$fitted.values), muEta = numeric(dm$dd[["n"]]),
var = numeric(dm$dd[["n"]]), resid = unname(glmFit$residuals),
sqrtXWt = as.matrix(numeric(dm$dd[["n"]])), sqrtrWt = numeric(dm$dd[["n"]]),
RZX = matrix(0, dm$dd[["q"]], dm$dd["p"]), RX = matrix(0, dm$dd["p"], dm$dd["p"]),
ghx = AGQlist[[1]], ghw = AGQlist[[2]],
p=pstart, phi=phistart,link.power=as.double(link.power),
bound.p=ctr$bound.p)
# run optimization
invisible(.Call("cpglmm_optimize",ans,ldtweedie))
if (ans@dims[["cvg"]] > 6)
warning(lme4:::convergenceMessage(ans@dims[["cvg"]]))
invisible(.Call("mer_update_RX",ans))
ans
}
cpglmm.control <- function(max.iter=300L,
max.fun=900L,
bound.p=c(1.01,1.99),
trace=0){
if (!is.numeric(max.iter) || max.iter <= 0)
stop("value of 'max.iter' must be > 0")
if (!is.numeric(max.fun) || max.fun <= 0)
stop("value of 'max.fun' must be > 0")
if (!is.numeric(bound.p) || length(bound.p)!=2)
stop("'bound.p' must be of length 2")
if (min(bound.p)<1)
stop("invalid lower bound in 'bound.p'")
if (!is.numeric(trace) && !is.logical(trace))
stop("'trace' must be logical or numeric")
bound.p <- sort(bound.p)
list(max.iter=as.integer(max.iter),
max.fun=as.integer(max.fun),
bound.p=as.numeric(bound.p),
trace=as.integer(trace))
}
# return -2*loglike given a list x
ldtweedie <- function(x){
-2*sum(log(dtweedie(y=x[[1]],mu=x[[2]],phi=x[[3]],p=x[[4]])))
}