[R] second try; writing user-defined GLM link function
Jessi Brown
jessilbrown at gmail.com
Fri Apr 21 03:55:13 CEST 2006
An update for all:
Using the combined contributions from Mark and Dr. Ripley, I've been
(apparently) successfully formulating both GLM's and GLMM's (using
the MASS function glmmPQL) analyzing my nest success data. The beta
parameter estimates look reasonable and the top models resemble those
from earlier analyses using a different nest survival analysis
approach.
However, I've now run into problems when trying to predict the daily
survival rates from fitted models. For example, for a model
considering nest height (NestHtZ) and nest age effects (MeanAge and
related terms; there is an overall cubic time trend in this model), I
tried to predict the daily survival rate for each day out of a 67 day
nest cycle (so MeanAge is a vector of 1 to 67) with mean nest height
(also a vector 67 rows in length; both comprise the matrix "nestday").
Here's what happens:
> summary(glm.24)
Call:
glm(formula = Success ~ NestHtZ + MeanAge + I(MeanAge^2) + I(MeanAge^3),
family = logexposure(ExposureDays = vc.apfa$days), data = vc.apfa)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.3264 -1.2341 0.6712 0.8905 1.5569
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 6.5742015 1.7767487 3.700 0.000215 ***
NestHtZ 0.6205444 0.2484583 2.498 0.012504 *
MeanAge -0.6018978 0.2983656 -2.017 0.043662 *
I(MeanAge^2) 0.0380521 0.0152053 2.503 0.012330 *
I(MeanAge^3) -0.0006349 0.0002358 -2.693 0.007091 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 174.86 on 136 degrees of freedom
Residual deviance: 233.82 on 132 degrees of freedom
AIC: 243.82
Number of Fisher Scoring iterations: 13
> glm.24.pred<-predict(glm.24,newdata=nestday, type="response", SE.fit=T)
Warning message:
longer object length
is not a multiple of shorter object length in: plogis(eta)^days
Can anyone tell me what I'm doing wrong?
cheers, Jessi Brown
More information about the R-help
mailing list