[R] second try; writing user-defined GLM link function

[Richard Chandler]

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