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

[Prof Brian Ripley]

> glm.24.pred<-predict(glm.24,newdata=nestday, type="response", SE.fit=T)

What is SE.fit?  The help says se.fit.  That _may_ be a problem.

However, I think the real problem is that the link function argument 
includes a reference to vc.apfa$days that is appropriate for fitting, not 
prediction. One way out might be (untested)

attach(va.apfa)
glm.24 <- glm(formula = Success ~ NestHtZ + MeanAge + I(MeanAge^2) + 
I(MeanAge^3), family = logexposure(ExposureDays = days), data = vc.apfa)
detach()
attach(nestday)
glm.24.pred<-predict(glm.24,newdata=nestday, type="response", SE.fit=T)
detach()

so that 'days' refers to the appropriate dataset both when fitting and 
predicting.

(This is bending glm() to do something it was not designed to do, so some 
hoop-jumping is needed.)


On Thu, 20 Apr 2006, Jessi Brown wrote:

> 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
>
> ______________________________________________
> R-help at 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
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595