From j.c.l.ooms at uu.nl  Tue Sep 30 21:10:16 2008
From: j.c.l.ooms at uu.nl (Jeroen Ooms)
Date: Tue, 30 Sep 2008 12:10:16 -0700 (PDT)
Subject: [Rd] weird behavior of drop1() for polr models (MASS)
In-Reply-To: <003801c92304$7ea05d50$7be117f0$@ca>
References: <19742120.post@talk.nabble.com>
	<003801c92304$7ea05d50$7be117f0$@ca>
Message-ID: <19748496.post@talk.nabble.com>


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.



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