5
$\begingroup$

I am running a GLM, using the function glm.nb (pscl package) trying to figure how what could influence a particular trait in several locations and years. The output as follow (with slight modification and removing the things beyond this question)

          Estimate Std. Error z value Pr(>|z|)    
          state2        -0.163190   0.807307  -0.202  0.00081    
          state3        -0.530588   0.783023  -0.678  0.00302 
          state4         0.942953   0.737328   1.279  0.00094    
          year2         -1.759102   0.733214  -2.399  0.01643 
          year3         -0.633870   0.662903  -0.956  0.33897

I understand the output as such: compared to year1, year 2 is statistically different, but year3 is not different from year1. but what about year2 and year3?

For the state factor, state2, state3 and state4 are all statistically different than state1. so I need a post-hoc test between state2 state3 and state4, am I right?

thanks

Improve this question
$\endgroup$
3
  • $\begingroup$ Could you please be more specific about your data? Explain all your variables and hypothesis, and if possible, add raw data (or at least FULL output from R). $\endgroup$
    – Ladislav Naďo
    Commented Jun 10, 2014 at 21:18
  • $\begingroup$ Thanks for the reply. basically I am trying to see what factors could explain variations in a given disease. I am using States, Years, Sex as well as body weight and some other body measurements. For each of the factors that had several levels which were significantly different from each other, should I run a post-hoc test between them? in other words, the output suggests that state2, 3 and 4 are all different from state1 (is this interpretation correct? if yes, is a post-hoc possible to see if they are different from each others? $\endgroup$
    – user36491
    Commented Jun 10, 2014 at 21:59
  • $\begingroup$ There's something very strange here: your p-values for the states don't make sense given their z-scores. However, for years, they do line up. Maybe a cut and paste error? $\endgroup$
    – Cliff AB
    Commented Jul 3, 2015 at 17:56

2 Answers 2

Reset to default
2
$\begingroup$

The summary function is not the best method to get post-hoc results. It is better to use something made for the task, like the emmeans package.

The following is a toy example. It uses the glm.nb function from the MASS package. MASS::glm.nb is supported by emmeans. I don't know if pscl::glm.nb would work as well.

The model in this example throws some errors. I'm ignoring them for this example.

To get an anova table you can use the anova function. Instead, I'm using car::Anova here. glm.nb isn't explicitly supported by car::Anova, but it appears to work okay.

Note that with emmeans you can compare treatments for a main effect or an interaction effect from the model. You can get estimates and p-values for individual contrasts (pairs) or have the results displayed as a compact letter display (cld).

if(!require(MASS)){install.packages("MASS")}
if(!require(emmeans)){install.packages("emmeans")}
if(!require(car)){install.packages("car")}

State = c(rep("state1",10), rep("state2",10), rep("state3",10), rep("state4",10))
Year  = rep(c("year1","year2"),20)
Value = c(1,3,3,5,9,11,2,7,4,5,6,7,8,9,9,10,7,8,8,10,
          3,2,3,2,4,2, 5,3,5,3,8,6,8,4,9,4, 8,5,8, 5)

Data = data.frame(State, Year, Value)

library(MASS)

model = glm.nb(Value ~ Year + State + Year:State, data = Data)

# anova(model)

library(car)
Anova(model)

library(emmeans)
marginal = emmeans(model, ~ State)

pairs(marginal)

### contrast          estimate        SE  df z.ratio p.value
### state1 - state2 -0.5216748 0.1829872 Inf  -2.851  0.0226
### state1 - state3  0.4488936 0.2335741 Inf   1.922  0.2188
### state1 - state4 -0.2565999 0.1942621 Inf  -1.321  0.5495
### state2 - state3  0.9705684 0.2135304 Inf   4.545  <.0001
### state2 - state4  0.2650749 0.1696353 Inf   1.563  0.4002
### state3 - state4 -0.7054935 0.2232682 Inf  -3.160  0.0086
###
### Results are averaged over the levels of: Year 
### Results are given on the log (not the response) scale. 
### P value adjustment: tukey method for comparing a family of 4 estimates 

cld(marginal, Letters=letters)

### State    emmean        SE  df asymp.LCL asymp.UCL .group
### state3 1.130882 0.1825757 Inf 0.7730397  1.488723  a    
### state1 1.579775 0.1456811 Inf 1.2942455  1.865305  ab   
### state4 1.836375 0.1285099 Inf 1.5845002  2.088250   bc  
### state2 2.101450 0.1107309 Inf 1.8844214  2.318479    c  

### Results are averaged over the levels of: Year 
### Results are given on the log (not the response) scale. 
### Confidence level used: 0.95 
### P value adjustment: tukey method for comparing a family of 4 estimates 
### significance level used: alpha = 0.05  

marginal = emmeans(model, ~ State:Year)
# pairs(marginal)
cld(marginal, Letters=letters)
Improve this answer
$\endgroup$
1
  • $\begingroup$ There have been some changes to the emmeans package since I wrote this. Particularly with the cld() function. You can still call a cld with emmeans, but you may need to call the multcomp package first. $\endgroup$
    – Sal Mangiafico
    Commented Nov 10, 2023 at 15:10
1
$\begingroup$

From what I see you are trying to use glm.nb which is a modification of the system function glm(). Such function includes estimation of the additional parameter, theta, for a Negative Binomial generalized linear model.

Are you sure that it is necessary? I do not know the exact details, but please consider the using of "classical" glm.

I think you should:

  1. construct a full model including interactions (one dependent variable and many explanatory variables)

  2. try to simplify the model into MAM (minimal adequate model) by removing non-significant interactions (or main effects). This link might help you in this procedure: How to get only desirable comparisons from post-hoc

When you have MAM next step is to take a look inside the model to see which group(s) are different from each other. You may either make all possible comparisons (this should work well if you have orthogonal design in you data), or see only comparisons in which you are interested by setting contrasts (see again the link given above).

R output in your question suggests that state2, 3 and 4 are all different from state1. So your interpretation is correct. This is (most probably) due to default setting of "treatment" contrast in your model which compares only first group with each other (how to change such contrast see again the given link above).

If you need any additional help, please do not hesitate to ask.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Ive used glm.nb as the error distribution follows negative binomial. I didnt simplify the model as I am only interested in statistically significant factors (State and Year). and needed to a post-hoc for each one of these alone. So not a full comprasions with all the factors, this is why I didnt think I need to simplify it. thanks. $\endgroup$
    – user36491
    Commented Jun 12, 2014 at 19:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.