[R] logistic regression

[Prof Brian Ripley]

On Mon, 11 Oct 2004, Heike Zimmermann wrote:

> Hello,
> 
> I have a problem concerning logistic regressions.  When I add a quadratic
> term to my linear model, I cannot draw the line through my scatterplot
> anymore, which is no problem without the quadratic term.
> In this example my binary response variable is "incidence", the explanatory
> variable is "sun":
> > model0<-glm(incidence~1,binomial)
> > model1<-glm(incidence~sun,binomial)
> > anova(model0,model1,test="Chi")
> Analysis of Deviance Table
> 
> Model 1: incidence ~ 1
> Model 2: incidence ~ sun
>   Resid. Df Resid. Dev  Df Deviance P(>|Chi|)
> 1       299     332.94                       
> 2       298     287.19   1    45.74 1.349e-11
> > qsun<-sun^2
> > model2<-glm(incidence~sun+qsun,binomial)
> > anova(model1,model2,test="Chi")
> Analysis of Deviance Table
> 
> Model 1: incidence ~ sun
> Model 2: incidence ~ sun + qsun
>   Resid. Df Resid. Dev  Df Deviance P(>|Chi|)
> 1       298    287.194                       
> 2       297    280.623   1    6.571     0.010
> 
> So the second, non-linear, model explains more than the first. 
> Now to create a graph I write:
> 
> > plot(sun,incidence)
> > min(sun)
> [1] 0
> > max(sun)
> [1] 90
> > xsun<-seq(0,90,1)
> 
> >lines(xsun,predict(model2,type="response",data.frame(sun=xsun)))

Why are you predicting using a new data frame containing just one of the 
two predictors?

> 
> Error in model.frame(formula, rownames, variables, varnames, extras,
> extranames,  : 
>         variable lengths differ
> > 
> 
> So this is the message I receive everytime I try to draw the fitted values
> of my model. I know for sure that exactly the same command works in S-Plus
> (with the same data). 

If so, you have been trapped by a bug in S-PLUS, for you have at no point
supplied the quadratic term needed for prediction.  I suspect more likely
that the command is accepted but gives incorrect answers.  See the
warnings in MASS, for example.

> How is ist possible to do this in R?

NB: note the use of spaces for readability.

mydata <- data.frame(sun)
model2 <- glm(incidence ~ sun + I(sun^2), binomial, data=mydata)
plot(sun, incidence)
lines(xsun, predict(model2, type="response", data.frame(sun=0:90)))

PLEASE use a meaningful subject line.

-- 
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