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 -- +++ GMX DSL Premiumtarife 3 Monate gratis* + WLAN-Router 0,- EUR* +++ Clevere DSL-Nutzer wechseln jetzt zu GMX: http://www.gmx.net/de/go/dsl
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
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
Hi Heike,
you'd probably want
lines(xsun, predict(model2, type="response",
newdata=data.frame(sun=xsun, qsun=xsun*xsun)))
I hope it helps.
Best,
Dimitris
----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven
Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/16/396887
Fax: +32/16/337015
Web: http://www.med.kuleuven.ac.be/biostat/
http://www.student.kuleuven.ac.be/~m0390867/dimitris.htm
----- Original Message -----
From: "Heike Zimmermann" <heikz at gmx.de>
To: <R-help at stat.math.ethz.ch>
Sent: Monday, October 11, 2004 4:32 PM
Subject: [R] logistic regression
> 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
>
> --
> +++ GMX DSL Premiumtarife 3 Monate gratis* + WLAN-Router 0,- EUR*
> +++
> Clevere DSL-Nutzer wechseln jetzt zu GMX:
> http://www.gmx.net/de/go/dsl
>
> ______________________________________________
> 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
>