[R] Poisson regression with gamma random effects
Daniel A. Powers
dpowers at mail.la.utexas.edu
Fri Feb 16 19:42:41 CET 2001
Goran -
(A convoluted solution...)
I coded a problem like this in gauss using an S-plus program by Naryan
Sastry at Rand (given below). It is set up for piecewise constant hazard
models, but can easily be tricked to modeling counts as it calls a poisson
regression (glm) in the M-step. Note that I have not tested this in R,
but the gauss code gives the same results as stata's xtpois on McCullough
and Nelder's ship data.
Cheers,
Dan
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Dan Powers
Associate Professor, Sociology
University of Texas at Austin
------------------------------------------------------------------------------
repwe <- function(formula, data, offset1=exposure, clust1=cluster,
interv1=interval, alpha.start=1, eps=1e-07, maxit=250, maxit1=20) {
# S-plus program to estimate a hazard model with a single gamma-distributed
# random effect using the E-M algorithm
# Author: Narayan Sastry
# Make sure formula includes no constant or offset, and that the first G
# variables describe the baseline hazard. Exposure should be in months.
# Sample call:
# model.out <- repwe(event~-1+factor(interval)+urb,data=ps.dat,clust1=cl,
# alpha.start=.5)
if(!is.data.frame(data)) stop("Data must be a data frame")
off.col <- match(deparse(substitute(offset1)),names(data))
exposure <- names(data[off.col])
clus.col <- match(deparse(substitute(clust1)),names(data))
cluster <- data[,clus.col]
inter.col <- match(deparse(substitute(interv1)),names(data))
interval <- data[,inter.col]
form1 <- as.character(formula)
formula2 <- paste(form1[2],"~",form1[3],"+offset(log(",exposure,"/12))")
event.col <- match(form1[2],names(data))
event <- data[,event.col]
tot.event <- tot.integ.haz <- rep(0,nrow(data))
N <- length(unique(cluster))
mu.hat1 <- mu.hat <- xi.hat <- t.ev.cl <- t.i.h.cl <- rep(0,N)
G <- max(interval)
off.set <- data[,off.col]/12
for(i in 1:N){
cl.num <- unique(cluster)[i]
tot.event[cluster==cl.num] <- t.ev.cl[i] <- sum(event[cluster==cl.num])
}
iter <- 0
phi.dif <- 1
cat("Iteration: ",format(iter)," alpha: ",format(alpha.start),"\n")
while((abs(phi.dif) > eps*10 | iter < 31) & iter <= maxit1) {
if(iter>0) {
old.coef <- haz.mod$coef
data$mu.hat1 <- (alpha.new + tot.event)/(alpha.new + tot.integ.haz) }
alpha <- if(iter==0) alpha.start else alpha.new
if(iter==1) formula2 <- paste(form1[2],"~",form1[3]," +
offset(log(",exposure,"/12) + log(mu.hat1))")
haz.mod <- glm(formula2,poisson,data,x=T)
log.L <- N*(alpha*log(alpha)-log(gamma(alpha))) + sum(event*predict(haz.mod)) +
sum(log(gamma(alpha+t.ev.cl))) - sum((alpha+t.ev.cl)*log(alpha+t.i.h.cl))
cat("logL = ",format(log.L),"\n")
integ.haz <- exp(predict(haz.mod)) * off.set
for(i in 1:N){
cl.num <- unique(cluster)[i]
tot.integ.haz[cluster==cl.num] <- t.i.h.cl[i] <-
sum(integ.haz[cluster==cl.num]) }
xi.hat <- digamma(alpha + t.ev.cl) - log(alpha + t.i.h.cl)
mu.hat <- (alpha + t.ev.cl)/(alpha + t.i.h.cl)
xi.mu.dif <- sum(xi.hat - mu.hat)
alpha.new <- newt.raph(alpha,xi.mu.dif,N,eps,maxit)
phi.dif3 <- 1
iter3 <- 0
while(abs(phi.dif3) > eps*10 & iter3 <= maxit) {
xi.hat <- digamma(alpha.new + t.ev.cl) - log(alpha.new + t.i.h.cl)
mu.hat <- (alpha.new + t.ev.cl)/(alpha.new + t.i.h.cl)
xi.mu.dif <- sum(xi.hat - mu.hat)
alpha.old <- alpha.new
alpha.new <- newt.raph(alpha.old,xi.mu.dif,N,eps,maxit)
phi.dif3 <- 1/alpha.new - 1/alpha.old
iter3 <- iter3 + 1
}
phi.dif <- 1/alpha - 1/alpha.new
iter <- iter + 1
cat("Iteration: ",format(iter)," alpha: ",format(alpha.new),"\n")
}
xi.hat <- digamma(alpha + t.ev.cl) - log(alpha + t.i.h.cl)
mu.hat <- (alpha + t.ev.cl)/(alpha + t.i.h.cl)
d2Q.dtheta2 <- t(haz.mod$x * haz.mod$weights) %*% haz.mod$x
p <- nrow(d2Q.dtheta2)
E.w.D <- cbind(rbind(d2Q.dtheta2,rep(0,p)),rep(0,p+1))
E.w.D[p+1,p+1] <- -N * (1/alpha - trigamma(alpha))
var.w <- mu.hat / (alpha + t.i.h.cl)
var.logw <- trigamma(alpha + t.ev.cl)
cov.w.logw <- mu.hat * (digamma(alpha + t.ev.cl + 1) -
digamma(alpha + t.ev.cl))
var.logww <- var.logw + var.w - 2 * cov.w.logw
cov.w.logww <- cov.w.logw - var.w
A <- rep(0,p)
if(p-G!=0) dHAZ.dbeta <- haz.mod$x[,(G+1):p] * integ.haz
dHAZ.dgam <- integ.haz
for(i in 1:N) {
cl.num <- unique(cluster)[i]
for(j in 1:p) {
A[j] <- if(j<=G) { sum(dHAZ.dgam[cluster==cl.num & interval==j]) }
else { if(p-G>1) sum(dHAZ.dbeta[cluster==cl.num,(j-G)])
else if(p-G==1) sum(dHAZ.dbeta[cluster==cl.num]) } }
T.i <- rbind(cbind(-A,rep(0,p)),c(0,1))
Var.W.i <- matrix(c(var.w[i],rep(cov.w.logww[i],2),var.logww[i]),2)
Var.U.i <- T.i %*% Var.W.i %*% t(T.i)
Var.U <- if(i==1) Var.U.i else Var.U + Var.U.i }
SVD <- svd (E.w.D - Var.U)
variances <- sqrt(diag( t(SVD$u %*% (t(SVD$v) * (1/SVD$d))) ))
log.L <- N*(alpha*log(alpha)-log(gamma(alpha))) + sum(event*predict(haz.mod)) +
sum(log(gamma(alpha+t.ev.cl))) - sum((alpha+t.ev.cl)*log(alpha+t.i.h.cl))
cat("logL = ",format(log.L),"\n")
final.mod <- haz.mod
final.mod$alpha <- alpha
cat("1/alpha: ",format(1/alpha),"\n")
cat("t(1/alpha):",format((1/alpha)/sqrt(alpha^(-4) * variances[p+1]^2)),"\n")
final.mod$alpha.se <- variances[p+1]
cat("s.e.(1/alpha): ",format(sqrt(alpha^(-4) * variances[p+1]^2)),"\n")
final.mod$phi.dif <- phi.dif
final.mod$iter <- iter
final.mod$std.errs <- variances[1:p]
cat("New t-values: ","\n")
cat(format(abs(haz.mod$coef/final.mod$std.errs)),"\n")
return(final.mod)
}
newt.raph <- function(alpha.old,xi.mu.dif,n,epsilon,maxit) {
iterations <- 0
phi.dif <- 1
alpha <- alpha.old
while(abs(phi.dif) > epsilon & iterations <= maxit) {
dQ.dalpha <- xi.mu.dif + n*(1 + log(alpha) - digamma(alpha))
d2Q.dalpha2 <- n*(1 / alpha - trigamma(alpha))
alpha.new <- alpha - dQ.dalpha / d2Q.dalpha2
phi.dif <- 1/alpha - 1/alpha.new
iterations <- iterations + 1
alpha <- alpha.new
}
return(alpha)
}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list