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

Tue Sep 30 23:14:09 CEST 2008

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