If you're not tied to glm(), you could use optim() which might allow
you some more flexibility. I modified the code from Venables and
Ripley (2002) pg. 445 to do this same thing a while back. If you use
it you should double check it with a statistician.
#link function
le.link <- function(theta, expos) {log(theta^(1/expos)/(1-theta^
(1/expos)))}
#inverse link
le.inv.link <- function(theta, expos) {((exp(theta)/(1+exp(theta)))
^expos)}
LogExpos <- function(formula, data=NULL, wt=rep(1, length(y)),
start = rep(0, p), expos, ...) {
x <- model.matrix(formula, data=data)
y <- model.frame(formula, data=data)[,1]
fmin <- function(beta, X, y, w) {
p <- le.inv.link(X %*% beta, expos=expos)
-sum(2*w*ifelse(y, log(p), log(1-p)))
}
gmin <- function(beta, X, y, w) {
eta <- X %*% beta; p <- le.inv.link(eta, expos=expos)
as.vector(-2 * (w*dlogis(eta) * ifelse(y, 1/p, -1/(1-
p)))) %*% X
}
if(is.null(dim(x))) dim(x) <- c(length(x), 1)
dn <- dimnames(x)[[2]]
if(!length(dn)) dn <- paste("Var", 1:ncol(x), sep="")
p <- ncol(x)
if(is.factor(y)) y <- (unclass(y) !=1)
fit <- optim(start, fmin, gmin, X = x, y = y, w = wt,
method="BFGS", hessian=T, ...)
names(fit$par) <- dn
fit$se <- sqrt(diag(solve(fit$hessian)))
names(fit$se) <- dn
z <- cbind(fit$par, fit$se); colnames(z) <- c
("Estimate", "std.err.")
cat("\nCoefficients:\n"); print(z)
cat("\nResidual Deviance:", format(fit$value), "\n")
cat("\nConvergence message:", fit$convergence, "\n")
invisible(fit)
}
#Example
set.seed(1)
yy <- rbinom(100, 1, .5)
x1 <- rnorm(100, 1)
x2 <- rnorm(100, 1)
time <- rpois(100, 3)
time <- ifelse(time==0, 1, time)
dat <- data.frame(yy, x1, x2, time)
LogExpos(yy ~ x1 + x2, dat, expos=time)
#Output
Coefficients:
Estimate std.err.
(Intercept) 1.25483347 0.2122538
x1 0.09026530 0.1430969
x2 -0.02418791 0.1195096
Residual Deviance: 158.4612
Convergence message: 0
Good luck,
Richard
--
Richard Chandler, M.S. candidate
Department of Natural Resources Conservation
UMass Amherst
(413)545-1237
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
