## Search this blog

### Multiple comparisons for generalized linear model in R

#### An example

##### Defining y and x
> y <- c (rpois(10,1), rpois(10,2), rpois(10,3), rpois(10,5))
> y
[1] 0 1 0 2 1 2 0 1 1 0 3 2 2 4 2 3 3 1 1 1 3 0 4 4 2 3 1 4 3 6 3 8 6 3 2 6 5 5 2 8
> x <- gl(4,10)
> x
[1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4
Levels: 1 2 3 4
##### Fitting a Poisson regression
> f <- glm(y~x, family=poisson)
> summary(f)

Call:
glm(formula = y ~ x, family = poisson)

Deviance Residuals:
Min        1Q    Median        3Q       Max
-2.44949  -0.90724   0.04533   0.52699   1.52242

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.2231     0.3536  -0.631 0.527945
x2            1.0116     0.4129   2.450 0.014277 *
x3            1.3218     0.3979   3.322 0.000895 ***
x4            1.7918     0.3819   4.692 2.71e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 67.967  on 39  degrees of freedom
Residual deviance: 34.884  on 36  degrees of freedom
AIC: 142.6

Number of Fisher Scoring iterations: 5
##### Multiple comparisons
> # install.packages("multcomp")
> require(multcomp)
> f.mc <- glht(f, linfct = mcp(x = "Tukey"))
> summary(f.mc)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: glm(formula = y ~ x, family = poisson)

Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
2 - 1 == 0   1.0116     0.4129   2.450  0.06433 .
3 - 1 == 0   1.3218     0.3979   3.322  0.00481 **
4 - 1 == 0   1.7918     0.3819   4.692  < 0.001 ***
3 - 2 == 0   0.3102     0.2807   1.105  0.67729
4 - 2 == 0   0.7802     0.2575   3.030  0.01198 *
4 - 3 == 0   0.4700     0.2327   2.019  0.17303
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

> confint(f.mc)

Simultaneous Confidence Intervals

Multiple Comparisons of Means: Tukey Contrasts

Fit: glm(formula = y ~ x, family = poisson)

Quantile = 2.5471
95% family-wise confidence level

Linear Hypotheses:
Estimate lwr      upr
2 - 1 == 0  1.01160 -0.04002  2.06322
3 - 1 == 0  1.32176  0.30822  2.33529
4 - 1 == 0  1.79176  0.81905  2.76447
3 - 2 == 0  0.31015 -0.40481  1.02512
4 - 2 == 0  0.78016  0.12436  1.43596
4 - 3 == 0  0.47000 -0.12281  1.06282

> coef(f.mc)
2 - 1     3 - 1     4 - 1     3 - 2     4 - 2     4 - 3
1.0116009 1.3217558 1.7917595 0.3101549 0.7801586 0.4700036
> vcov(f.mc)
2 - 1         3 - 1         4 - 1         3 - 2         4 - 2         4 - 3
2 - 1  1.704542e-01  1.249996e-01  1.249996e-01 -4.545455e-02 -4.545455e-02  9.714451e-17
3 - 1  1.249996e-01  1.583330e-01  1.249996e-01  3.333333e-02  6.938894e-17 -3.333333e-02
4 - 1  1.249996e-01  1.249996e-01  1.458330e-01 -2.775558e-17  2.083333e-02  2.083333e-02
3 - 2 -4.545455e-02  3.333333e-02 -2.775558e-17  7.878788e-02  4.545455e-02 -3.333333e-02
4 - 2 -4.545455e-02  6.938894e-17  2.083333e-02  4.545455e-02  6.628788e-02  2.083333e-02
4 - 3  9.714451e-17 -3.333333e-02  2.083333e-02 -3.333333e-02  2.083333e-02  5.416667e-02