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

From: Jeroen Ooms <j.c.l.ooms_at_uu.nl>
Date: Tue, 30 Sep 2008 12:10:16 -0700 (PDT)

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.

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");

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
Received on Tue 30 Sep 2008 - 19:18:35 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 Tue 30 Sep 2008 - 23:30:16 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