Re: [Rd] weird behavior of drop1() for polr models (MASS)

From: John Fox <jfox_at_mcmaster.ca>
Date: Tue, 30 Sep 2008 17:14:09 -0400

Dear Jeroen,

> -----Original Message-----
> From: r-devel-bounces_at_r-project.org [mailto:r-devel-bounces_at_r-project.org]
On
> Behalf Of Jeroen Ooms
> Sent: September-30-08 3:10 PM
> To: r-devel_at_r-project.org
> Subject: Re: [Rd] weird behavior of drop1() for polr models (MASS)
>
>
> Sorry for posting in -dev; I assumed a technical issue. Thanks very much
for
> the responses, I've realized that interaction effects of categorical
> predictors are more complicated than I thought. I've read the contrasts
help
> file and paragraph in the R manual about contrast, however I'm not quite
> sure I fully understand it yet.
>
> I tried Peter's suggestion of fitting "~Infl + Type + Infl:Type" and "Type
+
> Infl:Type", and indeed these result in identical models. For the first
> model, the main effects have 5 parameters and the interaction effect
> contains 6 parameters, for the second model the main effects have 3
> parameters and the interaction has 8 parameters. Do I understand it
> correctly that when the interaction effect is included in the model, R
just
> fits a parameter for every possible combination of the predictor
categories
> minus 1 for the reference group (so 3*4-1 parameters in this case)? And
that
> therefore the attribution of parameters to main or interaction predictors
is
> kind of arbitrary?
>
> I am aware that there are several ways of encoding the contrasts of
> categorical predictors, e.g. contr.sum(4); contr.treatment(4), however I
do
> not understand how the interaction effects are then coded (just products
for
> every combination?), and how this influences the SS type III values, when
it
> does in fact result in the same model.

If by type-III SS, you mean dropping particular columns of the model matrix and noting the increase in the residual sum of squares, this can indeed be affected by contrast coding, because when you drop columns corresponding to a term that's marginal to an included term, the resulting model is not variant with respect to changes in contrast coding. More fundamentally, the hypotheses that you're testing can change with the contrasts, and if you're using a set of contrasts that are not orthogonal in the row basis of the model matrix, the tests for "main effects" will not test what are reasonably interpreted as main effects.

>
> And does this same issue then also arise for all other anova models? For
> example, if I fit a normal model in glm (I have done the same thing
numerous
> times in SPSS glm univariate), am I then allowed to interpret the main
> effects SS3 values as genuine main effects? For example:
>
> options(contrasts = c("contr.sum", "contr.poly"));
> mymodel <- glm(Freq~Infl*Type,data=housing);
> drop1(mymodel,attributes(mymodel$terms)$term.labels,test="Chisq");

Why a chi-square test? Why use glm() for a linear model? contr.sum and contr.poly do indeed create contrasts that are orthogonal in the row basis of the model matrix, so you could interpret tests of main effects, though these imply averaging over the levels of the other factor. In the presence of an interaction, that might not be of much interest. What is the attraction of "type-III" tests as opposed to "type-II" tests?

These questions are sufficiently complicated that they're hard to address in an email. They are addressed in Sec. 8.2 and 10.4 of my Applied Regression text (either edition).

Regards,
 John

>
> Thanks again, all comments are very much appreciated!
>
> Jeroen.
>
>
>
>
>
>
>
> John Fox-6 wrote:
> >
> > Dear Jeroen,
> >
> > drop1() respects relationships of marginality among the terms in a model
> > and
> > won't drop a lower-order term (e.g., a main effect) when a higher-order
> > relative (e.g., an interaction to which the main effect is marginal) is
in
> > the model. If you want a complete "type-III" ANOVA table and are careful
> > with contrast coding (which you were not, since contr.treatment produces
> > contrasts that are not orthogonal in the row basis of the model matrix

--

> > you could use, e.g., contr.sum), then you can use the Anova() function
in
> > the car package, specifying type="III".
> >
> > Note: I'm responding to this message on the r-devel list where it was
> > posted, but this is really a question for r-help.
> >
> > I hope this helps,
> > John
> >
> > ------------------------------
> > John Fox, Professor
> > Department of Sociology
> > McMaster University
> > Hamilton, Ontario, Canada
> > web: socserv.mcmaster.ca/jfox
> >
> >
> >> -----Original Message-----
> >> From: r-devel-bounces_at_r-project.org
> >> [mailto:r-devel-bounces_at_r-project.org]
> > On
> >> Behalf Of Jeroen Ooms
> >> Sent: September-30-08 9:28 AM
> >> To: r-devel_at_r-project.org
> >> Subject: [Rd] weird behavior of drop1() for polr models (MASS)
> >>
> >>
> >> I would like to do a SS type III analysis on a proportional odds
logistic
> >> regression model. I use drop1(), but dropterm() shows the same
behaviour.
> > It
> >> works as expected for regular main effects models, however when the
model
> >> includes an interaction effect it seems to have problems with matching
> >> the
> >> parameters to the predictor terms. An example:
> >>
> >> library("MASS");
> >> options(contrasts = c("contr.treatment", "contr.poly"));
> >>
> >> house.plr1 <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data =
> >> housing);
> >>
drop1(house.plr1,attributes(house.plr1$terms)$term.labels,test="Chisq");
> >>
> >> house.plr2 <- polr(Sat ~ Infl * Type + Cont, weights = Freq, data =
> >> housing);
> >>
drop1(house.plr2,attributes(house.plr2$terms)$term.labels,test="Chisq");
> >>
> >> Notice that model 2 has a * instead of a + between predictors Infl and
> > Type.
> >> In model 1, estimated parameters are nicely attributed to the right
term,
> >> however in house.plr2, only 2 of the 4 terms are evaluated.
> >>
> >> I am using R version 2.7.2 (2008-08-25) i386-pc-mingw32, and
MASS_7.2-44.
> >>
> >> --
> >> View this message in context: http://www.nabble.com/weird-behavior-of-
> >> drop1%28%29-for-polr-models-%28MASS%29-tp19742120p19742120.html
> >> Sent from the R devel mailing list archive at Nabble.com.
> >>
> >> ______________________________________________
> >> R-devel_at_r-project.org mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-devel
> >
> > ______________________________________________
> > R-devel_at_r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-devel
> >
> >
>
> --
> View this message in context: http://www.nabble.com/weird-behavior-of-
> drop1%28%29-for-polr-models-%28MASS%29-tp19742120p19748496.html
> Sent from the R devel mailing list archive at Nabble.com.
>
> ______________________________________________
> R-devel_at_r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
______________________________________________ R-devel_at_r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Received on Tue 30 Sep 2008 - 21:17:12 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 01 Oct 2008 - 03:30:17 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.

list of date sections of archive