[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