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:
- I'm confused. Why is this happening?
- Which test-statistic should I be reporting when using these comparisons?
- 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
