[R] Post-hoc tests in MASS using glm.nb

[Bryony Tolhurst]

Dear Bill

Many thanks. I will try this.

One question: why is the attach()function problematic? I have always done it that way (well in my very limited R-using capacity!) as dictated by 'Statistics, an Introduction using R' by Michael Crawley. My dataset is called Side ('Side.txt') so how do I import this data without using attach(data)? I have tried:

side<-read.table('Side.txt',T)
attach(side)

instead of:

data<-read.table('Side.txt',T) 
attach(data)

But obviously I am still using the attach function, if not with 'data'!!

Thanks again

Bryony Tolhurst

-----Original Message-----
From: Bill.Venables at csiro.au [mailto:Bill.Venables at csiro.au] 
Sent: 17 May 2011 03:21
To: Bryony Tolhurst; r-help at r-project.org
Subject: RE: [R] Post-hoc tests in MASS using glm.nb

?relevel

Also, you might want to fit the models as follows

Model1 <- glm.nb(Cells ~ Cryogel*Day, data = myData)

myData2 <- within(myData, Cryogel <- relevel(Cryogel, ref = "2"))
Model2 <- update(Model1, data = myData1) 

&c

You should always spedify the data set when you fit a model if at all possible.  I would recommend you NEVER use attach() to put it on the search path, (under all but the most exceptional circumstances).

You could fit your model as 

Model0 <- glm.nv(Cells ~ interaction(Cryogel, Day) - 1, data = myData)

This will give you the subclass means as the regression coefficients.  You can then use vcov(Model0) to get the variance matrix and compare any two you like using directly calculated t-statistics.  This is pretty straightforward as well.

Bill Venables.


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of bryony
Sent: Tuesday, 17 May 2011 3:46 AM
To: r-help at r-project.org
Subject: [R] Post-hoc tests in MASS using glm.nb

I am struggling to generate p values for comparisons of levels (post-hoc
tests) in a glm with a negative binomial distribution

I am trying to compare cell counts on different days as grown on different media (e.g. types of cryogel) so I have 2 explanatory variables (Day and Cryogel), which are both factors, and an over-dispersed count variable (number of cells) as the response. I know that both variables are significant, and that there is a significant interaction between them.
However, I seem unable to generate multiple comparisons between the days and cryogels. 

So my model is 

Model1<-glm.nb(Cells~Cryogel+Day+Day:Cryogel)

The output gives me comparisons between levels of the factors relative to a reference level (Day 0 and Cryogel 1) as follows:

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)      1.2040     0.2743   4.389 1.14e-05 ***
Day14            3.3226     0.3440   9.658  < 2e-16 ***
Day28            3.3546     0.3440   9.752  < 2e-16 ***
Day7             3.3638     0.3440   9.779  < 2e-16 ***
Cryogel2         0.7097     0.3655   1.942  0.05215 .  
Cryogel3         0.7259     0.3651   1.988  0.04677 *  
Cryogel4         1.4191     0.3539   4.010 6.07e-05 ***
Day14:Cryogel2  -0.7910     0.4689  -1.687  0.09162 .  
Day28:Cryogel2  -0.5272     0.4685  -1.125  0.26053    
Day7:Cryogel2   -1.1794     0.4694  -2.512  0.01199 *  
Day14:Cryogel3  -1.0833     0.4691  -2.309  0.02092 *  
Day28:Cryogel3   0.1735     0.4733   0.367  0.71395    
Day7:Cryogel3   -1.0907     0.4690  -2.326  0.02003 *  
Day14:Cryogel4  -1.2834     0.4655  -2.757  0.00583 ** 
Day28:Cryogel4  -0.6300     0.4591  -1.372  0.16997    
Day7:Cryogel4   -1.3436     0.4596  -2.923  0.00347 ** 


HOWEVER I want ALL the comparisons e.g. Cryogel 2 versus 4, 3 versus 2 etc on each of the days. I realise that such multiple comparsions need to be approached with care to avoid Type 1 error, however it is easy to do this in other programmes (e.g. SPSS, Genstat) and I'm frustrated that it appears to be difficult in R. I have tried the glht (multcomp) function but it gives me the same results. I assume that there is some way of entering the data differently so as to tell R to use a different reference level each time and re-run the analysis for each level, but don't know what this is.
Please help!

Many thanks for your input

Bryony

--
View this message in context: http://r.789695.n4.nabble.com/Post-hoc-tests-in-MASS-using-glm-nb-tp3526934p3526934.html
Sent from the R help mailing list archive at Nabble.com.

______________________________________________
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.