John Sorkin
2009-Sep-08 01:41 UTC
[R] Omnibus test for main effects in the face ofaninteraction containing the main effects.
Daniel, When Group is entered as a factor, and the factor has two levels, the ANOVA table gives a p value for each level of the factor. What I am looking for is the omnibus p value for the factor, i.e. the test that the factor (with all its levels) improves the prediction of the outcome. You are correct that normally one could rely on the fact that the model Post<-Time+as.factor(Group)+as.factor(Group)*Time contains the model Post<-Time+as.factor(Group) and compare the two models using anova(model1,model2). However, my model is has a random effect, the comparison is not so easy. The REML comparions of nested random effects models is not valid when the fixed effects are not the same in the models, which is the essence of the problem in my case. In addition to the REML problem if one wants to perform an omnibus test for Group, one would want to compare nested models, one containing Group, and the other not containing group. This would suggest comparing Post<-Time+ as.factor(Group)*Time to Post<-Time+Group+as.factor(Group)*Time The quandry here is whether one should or not "allow" the first model as it is poorly specified - one term of the interaction, as.factor(Group)*Time, as.factor(Group) does not appear as a main effect - a no-no in model building. John John Sorkin M.D., Ph.D. Chief, Biostatistics and Informatics Baltimore VA Medical Center GRECC, University of Maryland School of Medicine Claude D. Pepper OAIC, University of Maryland Clinical Nutrition Research Unit, and Baltimore VA Center Stroke of Excellence University of Maryland School of Medicine Division of Gerontology Baltimore VA Medical Center 10 North Greene Street GRECC (BT/18/GR) Baltimore, MD 21201-1524 (Phone) 410-605-7119 (Fax) 410-605-7913 (Please call phone number above prior to faxing) jsorkin at grecc.umaryland.edu>>> "Daniel Malter" <daniel at umd.edu> 09/07/09 9:23 PM >>>John, your question is confusing. After reading it twice, I still cannot figure out what exactly you want to compare. Your model "a" is the unrestricted model, and model "b" is a restricted version of model "a" (i.e., b is a hiearchically reduced version of a, or put differently, all coefficients of b are in a with a having additional coefficients). Thus, it is appropriate to compare the models (also called nested models). Comparing c with a and d with a is also appropriate for the same reason. However, note that depedent on discipline, it may be highly unconventional to fit an interaction without all direct effects of the interacted variables (the reason for this being that you may get biased estimates). What you might consider is: 1. Run an intercept only model 2. Run a model with group and time 3. Run a model with group, time, and the interaction Then compare 2 to 1, and 3 to 2. This tells you whether including more variables (hierarchically) makes your model better. HTH, Daniel On a different note, if lme fits with "restricted maximum likelihood," I think I remember that you cannot compare them. You have to fit them with "maximum likelihood." I am pointing this out because lmer with restricted maximum likelihood by standard, so lme might too. ------------------------- cuncta stricte discussurus ------------------------- -----Urspr?ngliche Nachricht----- Von: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] Im Auftrag von John Sorkin Gesendet: Monday, September 07, 2009 4:00 PM An: r-help at r-project.org Betreff: [R] Omnibus test for main effects in the face of aninteraction containing the main effects. R 2.9.1 Windows XP UPDATE, Even my first suggestion anova(fita,fitb) is probably not appropriate as the fixed effects are different in the two model, so I don't even know how to perform the ombnibus test for the interaction! I am fitting a random effects ANOVA with two factors Group which has two levels and Time which has three levels: fita<-lme(Post~Time+factor(Group)+factor(Group)*Time, random=~1|SS,data=blah$alldata) I wantinteraction. I believe I can get the omnibus test for the interaction by running the model: fitb<-lme(Post~Time+factor(Group), random=~1|SS,data=blah$alldata) followed by anova(fita,fitb). How do I get the omnibus test for the main effects i.e. for Time and factor(Group)? I could drop each from the model, i.e. fitc<-lme(Post~ factor(Group)+factor(Group)*Time, random=~1|SS,data=blah$alldata) fitd<-lme(Post~Time+ factor(Group)*Time, random=~1|SS,data=blah$alldata) and then run anova(fita,fitc) anova(fita,fitd) but I don't like this option as it will have in interaction that contains a factor that is not included in the model as a main effect. How then do I get the omnibus test for Time and factor(Group)? Thanks John John David Sorkin M.D., Ph.D. Chief, Biostatistics and Informatics University of Maryland School of Medicine Division of Gerontology Baltimore VA Medical Center 10 North Greene Street GRECC (BT/18/GR) Baltimore, MD 21201-1524 (Phone) 410-605-7119 (Fax) 410-605-7913 (Please call phone number above prior to faxing) Confidentiality Statement: This email message, including any attachments, is for th...{{dropped:16}}
Ben Bolker
2009-Sep-08 02:14 UTC
[R] Omnibus test for main effects in the face ofaninteraction containing the main effects.
John Sorkin wrote:> > Daniel, > When Group is entered as a factor, and the factor has two levels, the > ANOVA table gives a p value for each level of the factor. What I am > looking for is the omnibus p value for the factor, i.e. the test that > the factor (with all its levels) improves the prediction of the outcome. > > You are correct that normally one could rely on the fact that the model > Post<-Time+as.factor(Group)+as.factor(Group)*Time > contains the model > Post<-Time+as.factor(Group) > and compare the two models using anova(model1,model2). However, my model > is has a random effect, the comparison is not so easy. The REML > compari[s]ons of nested random effects models is not valid when the fixed > effects are not the same in the models, which is the essence of the > problem in my case. >So why not use ML comparisons if this is what you want to do ... ?> In addition to the REML problem if one wants to perform an omnibus test > for Group, one would want to compare nested models, one containing > Group, and the other not containing group. This would suggest comparing > Post<-Time+ as.factor(Group)*Time to > Post<-Time+Group+as.factor(Group)*Time > The quand[a]ry here is whether one should or not "allow" the first model > as > it is poorly specified - one term of the interaction, > as.factor(Group)*Time, as.factor(Group) does not appear as a main effect > - a no-no in model building. > John >Agreed. I don't think the question makes sense -- "the overall effect of group" includes both the main effect & the interaction. (re) reading Venables, W. N. ?Exegeses on Linear Models.? 1998 International S-PLUS User Conference. Washington, DC, 1998. http://www.stats.ox.ac.uk/pub/MASS3/Exegeses.pdf. might be useful ... -- View this message in context: http://www.nabble.com/Re%3A-Omnibus-test-for-main-effects-in-the-face%09ofaninteraction%09containing-the-main-effects.-tp25338728p25338952.html Sent from the R help mailing list archive at Nabble.com.
Daniel Malter
2009-Sep-08 04:36 UTC
[R] Omnibus test for main effects in the faceofaninteraction containing the main effects.
John, as I wrote in the post sciptum, an anova on ML (but not REML) fitted models seems permissible (Faraway 2006, "Extending the linear model with R", p. 158). I am certainly not an expert on this and there are better "sources" of information on why and when (e.g., Deepayan Sarkar, Julian Faraway, Doug Bates, and certainly many others). As for your second problem, I think my statement remains valid, too. Because this may yield biased estimates, you should not allow for such a model. set.seed(1750) x1=rbinom(120,2,c(1/3,1/2)) x2=rbinom(120,1,1/2) e=rnorm(120,0,0.5) x1=factor(x1) x2=factor(x2) y=1*(x1==1)+(-1)*(x1==2)+1.5*(x2==1)+(x1==1)*(x2==1)*(-2)+(x1==2)*(x2==1)*(3 )+e reg=lm(y~x1+I((x1==1)*(x2==1))+I((x1==2)*(x2==1))) summary(reg) ##Gives biased estimates reg1=lm(y~x1+x2+x1:x2) summary(reg1) ##Gives unbiased estimates Therefore, I would not use the model with interaction effects if not all direct effects are included. And in fact, I never see this done in analyses in my field (business/economics). Daniel ------------------------- cuncta stricte discussurus ------------------------- -----Urspr?ngliche Nachricht----- Von: John Sorkin [mailto:Jsorkin at grecc.umaryland.edu] Gesendet: Monday, September 07, 2009 9:42 PM An: r-help at r-project.org; daniel at umd.edu Betreff: Re: [R] Omnibus test for main effects in the faceofaninteraction containing the main effects. Daniel, When Group is entered as a factor, and the factor has two levels, the ANOVA table gives a p value for each level of the factor. What I am looking for is the omnibus p value for the factor, i.e. the test that the factor (with all its levels) improves the prediction of the outcome. You are correct that normally one could rely on the fact that the model Post<-Time+as.factor(Group)+as.factor(Group)*Time contains the model Post<-Time+as.factor(Group) and compare the two models using anova(model1,model2). However, my model is has a random effect, the comparison is not so easy. The REML comparions of nested random effects models is not valid when the fixed effects are not the same in the models, which is the essence of the problem in my case. In addition to the REML problem if one wants to perform an omnibus test for Group, one would want to compare nested models, one containing Group, and the other not containing group. This would suggest comparing Post<-Time+ as.factor(Group)*Time to Post<-Time+Group+as.factor(Group)*Time The quandry here is whether one should or not "allow" the first model as it is poorly specified - one term of the interaction, as.factor(Group)*Time, as.factor(Group) does not appear as a main effect - a no-no in model building. John John Sorkin M.D., Ph.D. Chief, Biostatistics and Informatics Baltimore VA Medical Center GRECC, University of Maryland School of Medicine Claude D. Pepper OAIC, University of Maryland Clinical Nutrition Research Unit, and Baltimore VA Center Stroke of Excellence University of Maryland School of Medicine Division of Gerontology Baltimore VA Medical Center 10 North Greene Street GRECC (BT/18/GR) Baltimore, MD 21201-1524 (Phone) 410-605-7119 (Fax) 410-605-7913 (Please call phone number above prior to faxing) jsorkin at grecc.umaryland.edu>>> "Daniel Malter" <daniel at umd.edu> 09/07/09 9:23 PM >>>John, your question is confusing. After reading it twice, I still cannot figure out what exactly you want to compare. Your model "a" is the unrestricted model, and model "b" is a restricted version of model "a" (i.e., b is a hiearchically reduced version of a, or put differently, all coefficients of b are in a with a having additional coefficients). Thus, it is appropriate to compare the models (also called nested models). Comparing c with a and d with a is also appropriate for the same reason. However, note that depedent on discipline, it may be highly unconventional to fit an interaction without all direct effects of the interacted variables (the reason for this being that you may get biased estimates). What you might consider is: 1. Run an intercept only model 2. Run a model with group and time 3. Run a model with group, time, and the interaction Then compare 2 to 1, and 3 to 2. This tells you whether including more variables (hierarchically) makes your model better. HTH, Daniel On a different note, if lme fits with "restricted maximum likelihood," I think I remember that you cannot compare them. You have to fit them with "maximum likelihood." I am pointing this out because lmer with restricted maximum likelihood by standard, so lme might too. ------------------------- cuncta stricte discussurus ------------------------- -----Urspr?ngliche Nachricht----- Von: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] Im Auftrag von John Sorkin Gesendet: Monday, September 07, 2009 4:00 PM An: r-help at r-project.org Betreff: [R] Omnibus test for main effects in the face of aninteraction containing the main effects. R 2.9.1 Windows XP UPDATE, Even my first suggestion anova(fita,fitb) is probably not appropriate as the fixed effects are different in the two model, so I don't even know how to perform the ombnibus test for the interaction! I am fitting a random effects ANOVA with two factors Group which has two levels and Time which has three levels: fita<-lme(Post~Time+factor(Group)+factor(Group)*Time, random=~1|SS,data=blah$alldata) I wantinteraction. I believe I can get the omnibus test for the interaction by running the model: fitb<-lme(Post~Time+factor(Group), random=~1|SS,data=blah$alldata) followed by anova(fita,fitb). How do I get the omnibus test for the main effects i.e. for Time and factor(Group)? I could drop each from the model, i.e. fitc<-lme(Post~ factor(Group)+factor(Group)*Time, random=~1|SS,data=blah$alldata) fitd<-lme(Post~Time+ factor(Group)*Time, random=~1|SS,data=blah$alldata) and then run anova(fita,fitc) anova(fita,fitd) but I don't like this option as it will have in interaction that contains a factor that is not included in the model as a main effect. How then do I get the omnibus test for Time and factor(Group)? Thanks John John David Sorkin M.D., Ph.D. Chief, Biostatistics and Informatics University of Maryland School of Medicine Division of Gerontology Baltimore VA Medical Center 10 North Greene Street GRECC (BT/18/GR) Baltimore, MD 21201-1524 (Phone) 410-605-7119 (Fax) 410-605-7913 (Please call phone number above prior to faxing) Confidentiality Statement: This email message, including any attachments, is for th...{{dropped:17}}
Mark Difford
2009-Sep-08 07:09 UTC
[R] Omnibus test for main effects in the face ofaninteraction containing the main effects.
Hi John,>> When Group is entered as a factor, and the factor has two levels, the >> ANOVA table gives a p value for each level of the factor.This does not (normally) happen so you are doing something strange. ## From your first posting on this subject fita<-lme(Post~Time+factor(Group)+factor(Group)*Time, random=~1|SS,data=blah$alldata) To begin with, what is blah$alldata? lme() asks for a data frame if the formula interface is used. This looks like part of a list, or the column of a data frame. Have a look at the output below, from a syntactically correct model that has essentially the same structure as yours. The data set should be loaded with nlme so you can run it directly to see the result. Sex, with two levels, is not split in the anova table. ## str(Orthodont) anova( lme(distance ~ age * Sex, data = Orthodont, random = ~ 1) ) Surely this is essentially what you are looking for? If Sex is not already a factor, and it really is better to make it one when you build your data set, then you can use as.factor as you have done, with the same result. (Note: age * Sex expands to age + Sex + age:Sex, which equals the messy and unnecessary age + Sex + age * Sex.) Regards, Mark. John Sorkin wrote:> > Daniel, > When Group is entered as a factor, and the factor has two levels, the > ANOVA table gives a p value for each level of the factor. What I am > looking for is the omnibus p value for the factor, i.e. the test that > the factor (with all its levels) improves the prediction of the outcome. > > You are correct that normally one could rely on the fact that the model > Post<-Time+as.factor(Group)+as.factor(Group)*Time > contains the model > Post<-Time+as.factor(Group) > and compare the two models using anova(model1,model2). However, my model > is has a random effect, the comparison is not so easy. The REML > comparions of nested random effects models is not valid when the fixed > effects are not the same in the models, which is the essence of the > problem in my case. > > In addition to the REML problem if one wants to perform an omnibus test > for Group, one would want to compare nested models, one containing > Group, and the other not containing group. This would suggest comparing > Post<-Time+ as.factor(Group)*Time to > Post<-Time+Group+as.factor(Group)*Time > The quandry here is whether one should or not "allow" the first model as > it is poorly specified - one term of the interaction, > as.factor(Group)*Time, as.factor(Group) does not appear as a main effect > - a no-no in model building. > John > > > John Sorkin M.D., Ph.D. > Chief, Biostatistics and Informatics > Baltimore VA Medical Center GRECC, > University of Maryland School of Medicine Claude D. Pepper OAIC, > University of Maryland Clinical Nutrition Research Unit, and > Baltimore VA Center Stroke of Excellence > > University of Maryland School of Medicine > Division of Gerontology > Baltimore VA Medical Center > 10 North Greene Street > GRECC (BT/18/GR) > Baltimore, MD 21201-1524 > > (Phone) 410-605-7119 > (Fax) 410-605-7913 (Please call phone number above prior to faxing) > jsorkin at grecc.umaryland.edu >>>> "Daniel Malter" <daniel at umd.edu> 09/07/09 9:23 PM >>> > John, your question is confusing. After reading it twice, I still cannot > figure out what exactly you want to compare. > > Your model "a" is the unrestricted model, and model "b" is a restricted > version of model "a" (i.e., b is a hiearchically reduced version of a, > or > put differently, all coefficients of b are in a with a having additional > coefficients). Thus, it is appropriate to compare the models (also > called > nested models). > > Comparing c with a and d with a is also appropriate for the same reason. > However, note that depedent on discipline, it may be highly > unconventional > to fit an interaction without all direct effects of the interacted > variables > (the reason for this being that you may get biased estimates). > > What you might consider is: > 1. Run an intercept only model > 2. Run a model with group and time > 3. Run a model with group, time, and the interaction > > Then compare 2 to 1, and 3 to 2. This tells you whether including more > variables (hierarchically) makes your model better. > > HTH, > Daniel > > On a different note, if lme fits with "restricted maximum likelihood," I > think I remember that you cannot compare them. You have to fit them with > "maximum likelihood." I am pointing this out because lmer with > restricted > maximum likelihood by standard, so lme might too. > > ------------------------- > cuncta stricte discussurus > ------------------------- > > -----Urspr?ngliche Nachricht----- > Von: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] > Im > Auftrag von John Sorkin > Gesendet: Monday, September 07, 2009 4:00 PM > An: r-help at r-project.org > Betreff: [R] Omnibus test for main effects in the face of aninteraction > containing the main effects. > > R 2.9.1 > Windows XP > > UPDATE, > Even my first suggestion > anova(fita,fitb) is probably not appropriate as the fixed effects are > different in the two model, so I don't even know how to perform the > ombnibus > test for the interaction! > > > > I am fitting a random effects ANOVA with two factors Group which has two > levels and Time which has three levels: > fita<-lme(Post~Time+factor(Group)+factor(Group)*Time, > random=~1|SS,data=blah$alldata) > > I wantinteraction. I believe I can get the omnibus test for the > interaction by > running the model: > > fitb<-lme(Post~Time+factor(Group), random=~1|SS,data=blah$alldata) > followed > by anova(fita,fitb). > > How do I get the omnibus test for the main effects i.e. for Time and > factor(Group)? I could drop each from the model, i.e. > fitc<-lme(Post~ factor(Group)+factor(Group)*Time, > random=~1|SS,data=blah$alldata) > fitd<-lme(Post~Time+ factor(Group)*Time, > random=~1|SS,data=blah$alldata) > > and then run > anova(fita,fitc) > anova(fita,fitd) > but I don't like this option as it will have in interaction that > contains a > factor that is not included in the model as a main effect. How then do I > get > the omnibus test for Time and factor(Group)? > > Thanks > John > > > > > John David Sorkin M.D., Ph.D. > Chief, Biostatistics and Informatics > University of Maryland School of Medicine Division of Gerontology > Baltimore VA Medical Center > 10 North Greene Street > GRECC (BT/18/GR) > Baltimore, MD 21201-1524 > (Phone) 410-605-7119 > (Fax) 410-605-7913 (Please call phone number above prior to faxing) > > Confidentiality Statement: > This email message, including any attachments, is for th...{{dropped:16}} > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > >-- View this message in context: http://www.nabble.com/Re%3A-Omnibus-test-for-main-effects-in-the-face%09ofaninteraction%09containing-the-main-effects.-tp25338728p25340897.html Sent from the R help mailing list archive at Nabble.com.