[R] logistic regression

Sundar Dorai-Raj sundar.dorai-raj at PDF.COM
Mon Oct 11 16:52:31 CEST 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)))
> 
> 
> 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). How is ist possible to do this in R?
> 
> Thank you in advance, Heike
> 

This is because your newdata=data.frame(sun=xsun) does not contain a 
variable called "qsun". Try the following instead:

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


--sundar




More information about the R-help mailing list