# Inconsistent Z-scores from multiple comparisons of factors in glmmadmb models

First time poster, long time desperate reader. This place has been so helpful that I can usually find that someone has had the same problem as myself and had it answered here. However, I haven't found an answer to my problem below in a way that I understand.

I am having trouble finding consistent z-scores with glmmadmb when used in conjunction with glht. The z-scores attached to the reported p-values of a 6-level categorical variable match up with the p-values pinpointed on a z-table using the z-scores from the output of a glmmadmb model. But after a multiple comparisons, they do not. When testing a simpler 2-level factor the z-scores seem to be accurate for the original model and a multiple comparisons. Examples below. My query is three-fold:

1. I'm confused. Why is this happening?
2. Which test-statistic should I be reporting when using these comparisons?
3. How do I find the values of those test statistics?

I have been using glmmadmb (glmmADMB) to model bee visit duration (minutes) of bees to flowers with a categorical fixed factor (hour of morning: 6AM,7AM,8AM,9AM,10AM,11AM) and random factor (site). I chose to use a gamma distribution because the response variable is right-skewed, non-zero, and continuous. Code and output below:

sb<- glmmadmb(min~hour + (1|site), family = "gamma", data=sbsub)
summary(sb)

Call:
glmmadmb(formula = min ~ hour + (1 | site), data = sbsub, family = "gamma")

AIC: 134

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)    0.117      0.291    0.40  0.68708
hour7         -0.794      0.322   -2.47  0.01363 *
hour8         -1.112      0.315   -3.53  0.00042 ***
hour9         -0.642      0.323   -1.99  0.04660 *
hour10        -1.178      0.311   -3.79  0.00015 ***
hour11         0.350      0.329    1.07  0.28646
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Number of observations: total=93, site=11
Random effect variance(s):
Group=site
Variance StdDev
(Intercept)   0.3774 0.6143

Gamma shape parameter: 1.4184 (std. err.: 0.19914)

Log-likelihood: -58.9969


The p-values in the output can be found in a z-table after using the listed z-score from the outputs and multiplying the z-table value by 2...But, I did not just want to see how bee visit duration at each hour compares to the base level of 6AM. I wanted to know how visit duration compares across all combinations of hours. So I used glht and cld from the multcomp package and got the output below. Looks sexy, and informative. I could plop those letters right on a bar chart.

sb2<-glht(sb, linfct = mcp(hour = "Tukey"))
summary(sb2)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: glmmadmb(formula = min ~ hour + (1 | site), data = sbsub, family = "gamma")

Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
7 - 6 == 0   -0.79358    0.32168  -2.467   0.1206
8 - 6 == 0   -1.11241    0.31528  -3.528    <0.01 **
9 - 6 == 0   -0.64218    0.32272  -1.990   0.3215
10 - 6 == 0  -1.17755    0.31065  -3.791    <0.01 **
11 - 6 == 0   0.35035    0.32868   1.066   0.8793
8 - 7 == 0   -0.31883    0.44769  -0.712   0.9767
9 - 7 == 0    0.15140    0.45244   0.335   0.9993
10 - 7 == 0  -0.38398    0.44428  -0.864   0.9469
11 - 7 == 0   1.14393    0.46486   2.461   0.1232
9 - 8 == 0    0.47023    0.45019   1.044   0.8881
10 - 8 == 0  -0.06515    0.44124  -0.148   1.0000
11 - 8 == 0   1.46276    0.45547   3.212   0.0151 *
10 - 9 == 0  -0.53537    0.43020  -1.244   0.7919
11 - 9 == 0   0.99253    0.46182   2.149   0.2399
11 - 10 == 0  1.52790    0.46197   3.307   0.0110 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

cld(sb2)
6    7    8    9   10   11
"b" "ab"  "a" "ab"  "a"  "b"


But, these z-scores don't seem to mean anything. Looking them up in the z-table doesn't reveal the same p-value/2 in the output. So, reporting them seems moot. However, I feel I should report some kind of test-statistic to accompany the p-values? This also seems to happen when using glmer with a binomial distribution modeling predation and parasitism. In that case, the z-scores from a multiple comparisons between a 3-level categorical variable and a 2-level categorical variable do not have accurate z-scores.

I tried it with another smaller pollination model testing visit duration of bees in male and female flowers.

sb<- glmmadmb(min~sex + (1|site), family = "gamma", data=sbsub)
summary(sb)

Call:
glmmadmb(formula = min ~ sex + (1 | site), data = sbsub, family = "gamma")

AIC: 154.4

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   -0.154      0.257   -0.60    0.550
sex1          -0.367      0.215   -1.71    0.087 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Number of observations: total=93, site=11
Random effect variance(s):
Group=site
Variance StdDev
(Intercept)   0.4517 0.6721

Gamma shape parameter: 1.102 (std. err.: 0.15293)

Log-likelihood: -73.1965

sb2<-glht(sb, linfct = mcp(sex = "Tukey"))
summary(sb2)

Fit: glmmadmb(formula = min ~ sex + (1 | site), data = sbsub, family = "gamma")

Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
1 - 0 == 0  -0.3669     0.2147  -1.709   0.0875 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

cld(sb2)
0    1
"a"  "a"


The p-values/2 in the output of both the original model and the multiple comparison can be found after using the listed z-score from the outputs to find a value in the z-table.

Not sure what's happening here. Any thoughts?

Much appreciated, Nava

• The $p$-values in the output of glht are meaningful and you can report them: The reported $p$-values are adjusted for multiple comparison (the output states it explicitly). In the last output, you only compare two groups and thus, the $p$-value has not to be adjusted. – COOLSerdash Jul 16 '13 at 8:11
• Ok. The relationships reported to be different do look quite different, so that's good. How about the test-statistic, z? – Nava Jul 16 '13 at 22:34
• What do you mean? The $z$-values are the estimates divided by the corresponding standard errors. – COOLSerdash Jul 16 '13 at 22:42
• Ah-ha! I wasn't aware of that. I was double-checking the outputs with a z-table and they were not corresponding with the reported p-values. I assumed reviewers would do the same thing. In the cases where I am to report adjusted p-values from a multiple comparison, is it prudent to include the z-score as a test-statistic? Thanks for your help – Nava Jul 16 '13 at 23:57
• You may want to have a look at that. Yes, I would include the $z$ or $t$ values as well as the adjusted $p$-values in a paper. The adjustment in multiple comparisons is a standard and reviewers should be aware of that. I would also state the adjustment procedure, which was "single-step" in this case. For technical details, see here on page 2 (the `pmv function). Please ask if you have further questions. – COOLSerdash Jul 17 '13 at 6:55