Excellent!
So to provide the final summary for everyone's sake, based on Dr.
Ripley's revisions, and now a revised logexposure "family" function to
use with the improved link function (until 2.4.0 makes it simpler) :
# Logistical Exposure Link Function
# See Shaffer, T. 2004. A unifying approach to analyzing nest success.
# Auk 121(2): 526-540.
logexp <- function(days = 1)
{
linkfun <- function(mu) qlogis(mu^(1/days))
linkinv <- function(eta) plogis(eta)^days
mu.eta <- function(eta)
days*.Call("logit_mu_eta", eta,
PACKAGE = "stats")*plogis(eta)^(days-1)
valideta <- function(eta) TRUE
link <- paste("logexp(", days, ")", sep="")
structure(list(linkfun = linkfun, linkinv = linkinv,
mu.eta = mu.eta, valideta = valideta, name = link),
class = "link-glm")
}
# Modified binomial family function (that allows logexp link function)
logexposure<-function (link="logexp",ExposureDays) {
variance <- function(mu) mu * (1 - mu)
validmu <- function(mu) all(mu > 0) && all(mu < 1)
dev.resids <- function(y, mu, wt) .Call("binomial_dev_resids",
y, mu, wt, PACKAGE = "stats")
aic <- function(y, n, mu, wt, dev) {
m <- if (any(n > 1))
n
else wt
-2 * sum(ifelse(m > 0, (wt/m), 0) * dbinom(round(m *
y), round(m), mu, log = TRUE))
}
initialize <- expression({
if (NCOL(y) == 1) {
if (is.factor(y)) y <- y != levels(y)[1]
n <- rep.int(1, nobs)
if (any(y < 0 | y > 1)) stop("y values must be 0 <= y <= 1")
mustart <- (weights * y + 0.5)/(weights + 1)
m <- weights * y
if (any(abs(m - round(m)) > 0.001))
warning("non-integer successes in a binomial glm!")
} else if (NCOL(y) == 2) {
if (any(abs(y - round(y)) > 0.001))
warning("non-integer counts in a binomial glm!")
n <- y[, 1] + y[, 2]
y <- ifelse(n == 0, 0, y[, 1]/n)
weights <- weights * n
mustart <- (n * y + 0.5)/(n + 1)
} else stop("for the binomial family,",
" y must be a vector of 0 and 1's\n",
"or a 2 column matrix where col 1 is",
" no. successes and col 2 is no. failures")
})
structure(list(family="binomial", link=logexp(ExposureDays),
linkfun=logexp(ExposureDays)$linkfun,
linkinv=logexp(ExposureDays)$linkinv, variance=variance,
dev.resids=dev.resids, aic=aic,
mu.eta=logexp(ExposureDays)$mu.eta, initialize=initialize,
validmu=validmu, valideta=logexp$valideta), class = "family")
}
#Example
nestdata<-read.table("http://www.branta.org/ShafferChatNestData/chat.txt";)
chat.glm.logexp<-glm(survive/trials~parastat+nest_ht*patsize,family=logexposure(ExposureDays=nestdata$expos),data=nestdata)
# if you have MASS installed
library(MASS)
chat.step<-stepAIC(chat.glm.logexp,scope=list(upper=~parastat+nest_ht*patsize,lower=~1))
chat.step$anova
summary(chat.step)
Mark Herzog, Ph.D.
Program Leader, San Francisco Bay Research
Wetland Division, PRBO Conservation Science
4990 Shoreline Highway 1
Stinson Beach, CA 94970
(415) 893-7677 x308
[EMAIL PROTECTED]
Prof Brian Ripley wrote:
> A couple more comments:
>
> The null deviance is wrong here, as the code assumes that the link
> function maps constant vectors to constant vectors, which it does not
> here. You can circumvent that by setting an offset.
>
> Even setting dispersion = 1 I get slightly different se's.
>
> Here's a more robust way to specify the link:
>
> logexp <- function(days = 1)
> {
> linkfun <- function(mu) qlogis(mu^(1/days))
> linkinv <- function(eta) plogis(eta)^days
> mu.eta <- function(eta)
> days*.Call("logit_mu_eta", eta, PACKAGE =
> "stats")*plogis(eta)^(days-1)
> valideta <- function(eta) TRUE
> link <- paste("logexp(", days, ")", sep="")
> structure(list(linkfun = linkfun, linkinv = linkinv,
> mu.eta = mu.eta, valideta = valideta, name = link),
> class = "link-glm")
> }
>
> and in 2.4.0 you will able to use binomial(logexp(nestdata$exposure)).
>
______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
