HemiMarsh
Nick Masto
3/20/2021
1. Background
1.1 The hemi-marsh concept
The idea of the hemi-marsh has revolutionized waterfowl management. A theory spear-headed in the late 1960s-early 1970s, predicted that ducks and other waterfowl preferentially selected marshes with approximately equal vegetation-to-water interspersion (Weller and Spatcher 1965, Weller and Fredrickson 1974). However, it wasn’t until Kaminski and Prince (1981) tested this hypothesis experimentally at Delta Marsh Research Station in Manitoba, that the hemi-marsh concept was universally accepted by the waterfowl community as “best-practice”.
1.2 Analytical Approach in Kaminski and Prince (1981)
Kaminski and Prince (1981) analyzed their experiment using traditional means testing Analysis of Variance (ANOVA). However, the experimental design in which avian abundance and diversity was collected in Kaminski and Prince (1981), called for a repeated measures ANOVA where dabbling duck indicated breeding pairs were counted each 3 times per day, once per week, across 5-6 weeks per year. In the mid-1970s when the data were being analyzed, repeated measures ANOVA was not programatically or computational available. We wish to reanalyze dabbling duck density and diversity accounting for the experimental design; we expect, reanalysis may refine aspects of the hemi-marsh theory.
2. Data and Design
The data were collected in 1977 and 1978 at the Delta Marsh Research Station. Experimental 1-ha impoundments were established and randomly assigned treatments of 30:70, 50:50, 70:30 water-to-vegetation interspersion (water:veg) by either mowing or tilling (treatment). The authors counted indicated breeding pairs (ibp) of blue-winged teal, gadwall, mallards, northern pintail, northern shoveler in experimental plots. IBPs were counted in the morning, afternoon, and evening (period) across 6 weeks and 5 weeks in 1977 and 1978, respectively (week).
dat <- read.csv("hemimarsh.csv")
str(dat)## 'data.frame': 990 obs. of 7 variables:
## $ ip : int 8 3 4 2 4 3 5 4 6 1 ...
## $ species : chr "MALL" "MALL" "MALL" "MALL" ...
## $ treatment: chr "mow" "mow" "mow" "mow" ...
## $ veg.water: chr "a" "a" "a" "a" ...
## $ year : int 1977 1977 1977 1977 1977 1977 1977 1977 1977 1977 ...
## $ week : int 1 2 3 4 5 6 1 2 3 4 ...
## $ period : chr "morning" "morning" "morning" "morning" ...#look at and change variables to factors
names(dat) ## [1] "ip" "species" "treatment" "veg.water" "year" "week"
## [7] "period"factors <- c(2:7)
dat[,factors] <- lapply(dat[,factors], factor)
str(dat)## 'data.frame': 990 obs. of 7 variables:
## $ ip : int 8 3 4 2 4 3 5 4 6 1 ...
## $ species : Factor w/ 5 levels "BWTE","GADW",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ treatment: Factor w/ 2 levels "mow","till": 1 1 1 1 1 1 1 1 1 1 ...
## $ veg.water: Factor w/ 3 levels "a","b","c": 1 1 1 1 1 1 1 1 1 1 ...
## $ year : Factor w/ 2 levels "1977","1978": 1 1 1 1 1 1 1 1 1 1 ...
## $ week : Factor w/ 6 levels "1","2","3","4",..: 1 2 3 4 5 6 1 2 3 4 ...
## $ period : Factor w/ 3 levels "afternoon","evening",..: 3 3 3 3 3 3 1 1 1 1 ...2.1 Transforming IBPs
Kaminski and Prince (1981) transformed indicated breeding pairs to convert to IBPs per hectare by adding species count by 0.5 and dividing by 3. We followed their methodology exactly.
dat <- dat %>%
mutate(ip.ha = (ip + 0.5) / 3)
boxplot(ip.ha ~ species,
data = dat,
main = "IBPs by species",
xlab = "Dabbling Duck Species",
ylab = "Indicated Breeding Pairs",
col = "orange",
border = "purple"
)
#remove that MALL outlier
#dat <- dat %>%
# filter(ip.ha < 10)3. Means Testing with IBPs/ha as response
Because Kaminski and Prince (1981) transformed all IBPs by square-rooting to fit assumptions of homogenous residuals, we will also square-root IBPs. We will examine residuals for each species. For this exercise, we also will replicate Kaminski and Prince (1981) by veg:water, treatment, period, and their interactions as fixed variables. We also will analyze data for each year separately. Thus, any differences in our results and the original authors should be a direct consequence of including week as a random variable in all models.
3.1 Blue-winged Teal
Below and from the plots, we see that no evidence is found for interactive effects. However, blue-winged teal there is evidence for veg:water ratio and treatment in 1977 and an effect found for all three fixed variables in 1978. Point estimates across all years was greater for 50:50 veg to water interspersion and greater BWTE densities were generally found in mowed versus tilled plots.
#1977----------
bwte.m1977 <- lmer(sqrt(ip.ha) ~ veg.water * treatment * period + (1|week), data = bwte1977)
plot(bwte.m1977,
main = "BWTE 1977 Fit",
xlab = "Residuals",
ylab = "Fitted")
anova(bwte.m1977)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## veg.water 3.8853 1.94265 2 85 14.1644 4.907e-06 ***
## treatment 0.6553 0.65530 1 85 4.7780 0.03158 *
## period 0.7771 0.38855 2 85 2.8330 0.06440 .
## veg.water:treatment 0.1250 0.06250 2 85 0.4557 0.63552
## veg.water:period 0.4499 0.11247 4 85 0.8200 0.51595
## treatment:period 0.3034 0.15168 2 85 1.1059 0.33562
## veg.water:treatment:period 0.1066 0.02664 4 85 0.1942 0.94079
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(bwte.m1977, pairwise~veg.water, type = "response")## NOTE: Results may be misleading due to involvement in interactions## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## veg.water response SE df lower.CL upper.CL
## a 1.59 0.460 5.86 0.664 2.93
## b 2.45 0.570 5.86 1.251 4.06
## c 1.23 0.404 5.86 0.438 2.43
##
## Results are averaged over the levels of: treatment, period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## a - b -0.303 0.0873 85 -3.476 0.0023
## a - c 0.153 0.0873 85 1.753 0.1918
## b - c 0.456 0.0873 85 5.229 <.0001
##
## Results are averaged over the levels of: treatment, period
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmeans(bwte.m1977, pairwise~veg.water|treatment, type = "response")## NOTE: Results may be misleading due to involvement in interactions
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## treatment = mow:
## veg.water response SE df lower.CL upper.CL
## a 1.32 0.442 7.27 0.487 2.56
## b 2.35 0.590 7.27 1.173 3.94
## c 1.04 0.393 7.27 0.326 2.17
##
## treatment = till:
## veg.water response SE df lower.CL upper.CL
## a 1.90 0.529 7.27 0.857 3.34
## b 2.56 0.614 7.27 1.317 4.20
## c 1.43 0.460 7.27 0.558 2.72
##
## Results are averaged over the levels of: period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## treatment = mow:
## contrast estimate SE df t.ratio p.value
## a - b -0.385 0.123 85 -3.119 0.0069
## a - c 0.127 0.123 85 1.029 0.5608
## b - c 0.512 0.123 85 4.147 0.0002
##
## treatment = till:
## contrast estimate SE df t.ratio p.value
## a - b -0.222 0.123 85 -1.797 0.1768
## a - c 0.179 0.123 85 1.450 0.3201
## b - c 0.401 0.123 85 3.247 0.0047
##
## Results are averaged over the levels of: period
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmip(bwte.m1977, ~veg.water | treatment, CIs = TRUE, type = "response",
linearg = list(linetype = NA),
CIarg = list(lwd = 1, alpha = 1))## NOTE: Results may be misleading due to involvement in interactions
#1978----------
bwte.m1978 <- lmer(sqrt(ip.ha) ~ veg.water * treatment * period + (1|week), data = bwte1978)
plot(bwte.m1978,
main = "BWTE 1978 Fit",
xlab = "Residuals",
ylab = "Fitted")
anova(bwte.m1978)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## veg.water 2.09274 1.04637 2 68 11.8883 3.736e-05 ***
## treatment 0.50288 0.50288 1 68 5.7135 0.0196113 *
## period 1.58625 0.79312 2 68 9.0110 0.0003377 ***
## veg.water:treatment 0.24068 0.12034 2 68 1.3672 0.2617263
## veg.water:period 0.09555 0.02389 4 68 0.2714 0.8954241
## treatment:period 0.45125 0.22563 2 68 2.5634 0.0844661 .
## veg.water:treatment:period 0.26031 0.06508 4 68 0.7394 0.5684135
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(bwte.m1978, pairwise~veg.water, type = "response")## NOTE: Results may be misleading due to involvement in interactions
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## veg.water response SE df lower.CL upper.CL
## a 0.866 0.132 10.4 0.598 1.18
## b 1.602 0.180 10.4 1.228 2.03
## c 0.913 0.136 10.4 0.637 1.24
##
## Results are averaged over the levels of: treatment, period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## a - b -0.3352 0.0766 68 -4.376 0.0001
## a - c -0.0248 0.0766 68 -0.324 0.9437
## b - c 0.3103 0.0766 68 4.051 0.0004
##
## Results are averaged over the levels of: treatment, period
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmeans(bwte.m1978, pairwise~treatment, type = "response")## NOTE: Results may be misleading due to involvement in interactions
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## treatment response SE df lower.CL upper.CL
## mow 1.266 0.144 6.88 0.948 1.63
## till 0.952 0.125 6.88 0.680 1.27
##
## Results are averaged over the levels of: veg.water, period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## mow - till 0.15 0.0625 68 2.390 0.0196
##
## Results are averaged over the levels of: veg.water, period
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-rogeremmeans(bwte.m1978, pairwise~period, type = "response")## NOTE: Results may be misleading due to involvement in interactions
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## period response SE df lower.CL upper.CL
## afternoon 1.390 0.168 10.4 1.043 1.79
## evening 1.221 0.157 10.4 0.898 1.59
## morning 0.753 0.123 10.4 0.504 1.05
##
## Results are averaged over the levels of: veg.water, treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## afternoon - evening 0.0741 0.0766 68 0.968 0.5997
## afternoon - morning 0.3113 0.0766 68 4.064 0.0004
## evening - morning 0.2372 0.0766 68 3.096 0.0079
##
## Results are averaged over the levels of: veg.water, treatment
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmeans(bwte.m1978, pairwise~veg.water|treatment|period, type = "response")## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## treatment = mow, period = afternoon:
## veg.water response SE df lower.CL upper.CL
## a 1.443 0.337 60.2 0.847 2.196
## b 2.304 0.426 60.2 1.530 3.236
## c 1.503 0.344 60.2 0.893 2.271
##
## treatment = till, period = afternoon:
## veg.water response SE df lower.CL upper.CL
## a 0.772 0.247 60.2 0.358 1.345
## b 1.922 0.389 60.2 1.222 2.780
## c 0.747 0.243 60.2 0.340 1.311
##
## treatment = mow, period = evening:
## veg.water response SE df lower.CL upper.CL
## a 1.430 0.336 60.2 0.837 2.180
## b 1.711 0.367 60.2 1.055 2.525
## c 1.321 0.323 60.2 0.754 2.046
##
## treatment = till, period = evening:
## veg.water response SE df lower.CL upper.CL
## a 0.647 0.226 60.2 0.274 1.178
## b 1.713 0.368 60.2 1.057 2.528
## c 0.745 0.242 60.2 0.339 1.309
##
## treatment = mow, period = morning:
## veg.water response SE df lower.CL upper.CL
## a 0.395 0.177 60.2 0.121 0.827
## b 1.002 0.281 60.2 0.519 1.643
## c 0.810 0.253 60.2 0.384 1.395
##
## treatment = till, period = morning:
## veg.water response SE df lower.CL upper.CL
## a 0.765 0.246 60.2 0.353 1.336
## b 1.150 0.301 60.2 0.626 1.831
## c 0.532 0.205 60.2 0.201 1.021
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## treatment = mow, period = afternoon:
## contrast estimate SE df t.ratio p.value
## a - b -0.3168 0.188 68 -1.688 0.2170
## a - c -0.0248 0.188 68 -0.132 0.9904
## b - c 0.2919 0.188 68 1.556 0.2718
##
## treatment = till, period = afternoon:
## contrast estimate SE df t.ratio p.value
## a - b -0.5076 0.188 68 -2.705 0.0232
## a - c 0.0148 0.188 68 0.079 0.9966
## b - c 0.5224 0.188 68 2.784 0.0188
##
## treatment = mow, period = evening:
## contrast estimate SE df t.ratio p.value
## a - b -0.1124 0.188 68 -0.599 0.8211
## a - c 0.0462 0.188 68 0.246 0.9671
## b - c 0.1586 0.188 68 0.845 0.6763
##
## treatment = till, period = evening:
## contrast estimate SE df t.ratio p.value
## a - b -0.5047 0.188 68 -2.690 0.0241
## a - c -0.0591 0.188 68 -0.315 0.9469
## b - c 0.4456 0.188 68 2.375 0.0525
##
## treatment = mow, period = morning:
## contrast estimate SE df t.ratio p.value
## a - b -0.3723 0.188 68 -1.984 0.1239
## a - c -0.2715 0.188 68 -1.447 0.3230
## b - c 0.1008 0.188 68 0.537 0.8532
##
## treatment = till, period = morning:
## contrast estimate SE df t.ratio p.value
## a - b -0.1974 0.188 68 -1.052 0.5471
## a - c 0.1453 0.188 68 0.774 0.7200
## b - c 0.3426 0.188 68 1.826 0.1688
##
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmip(bwte.m1978, ~veg.water | treatment | period, CIs = TRUE, type = "response",
linearg = list(linetype = NA),
CIarg = list(lwd = 1, alpha = 1))
3.2 Gadwall
In 1977, gadwall indicated breeding pairs were unaffected by veg-to-water interspersion, treatment, or survey timing. However, gadwall were observed in greater densities in 50:50 veg-to-water plots than others.
#1977----------
gadw.m1977 <- lmer(sqrt(ip.ha) ~ veg.water * treatment * period + (1|week), data = gadw1977)
#that fit isn't great but it'll probably do
#ANOVA is fairly robust against departures of normal residuals
plot(gadw.m1977,
main = "GADW 1977 Fit",
xlab = "Residuals",
ylab = "Fitted")
anova(gadw.m1977)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## veg.water 0.063357 0.031679 2 85 1.2558 0.29006
## treatment 0.046473 0.046473 1 85 1.8423 0.17827
## period 0.040594 0.020297 2 85 0.8046 0.45062
## veg.water:treatment 0.028208 0.014104 2 85 0.5591 0.57380
## veg.water:period 0.018393 0.004598 4 85 0.1823 0.94700
## treatment:period 0.128825 0.064413 2 85 2.5535 0.08376 .
## veg.water:treatment:period 0.089758 0.022439 4 85 0.8896 0.47386
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(gadw.m1977, pairwise~veg.water, type = "response")## NOTE: Results may be misleading due to involvement in interactions## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## veg.water response SE df lower.CL upper.CL
## a 0.230 0.0452 7.99 0.138 0.346
## b 0.291 0.0508 7.99 0.185 0.419
## c 0.263 0.0483 7.99 0.164 0.386
##
## Results are averaged over the levels of: treatment, period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## a - b -0.0592 0.0374 85 -1.581 0.2594
## a - c -0.0332 0.0374 85 -0.887 0.6499
## b - c 0.0260 0.0374 85 0.694 0.7677
##
## Results are averaged over the levels of: treatment, period
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmeans(gadw.m1977, pairwise~veg.water|treatment, type = "response")## NOTE: Results may be misleading due to involvement in interactions
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## treatment = mow:
## veg.water response SE df lower.CL upper.CL
## a 0.220 0.0507 13.6 0.125 0.343
## b 0.282 0.0573 13.6 0.172 0.419
## c 0.220 0.0507 13.6 0.125 0.343
##
## treatment = till:
## veg.water response SE df lower.CL upper.CL
## a 0.240 0.0530 13.6 0.140 0.368
## b 0.300 0.0591 13.6 0.186 0.440
## c 0.310 0.0601 13.6 0.194 0.453
##
## Results are averaged over the levels of: period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## treatment = mow:
## contrast estimate SE df t.ratio p.value
## a - b -0.06124 0.0529 85 -1.157 0.4822
## a - c 0.00000 0.0529 85 0.000 1.0000
## b - c 0.06124 0.0529 85 1.157 0.4822
##
## treatment = till:
## contrast estimate SE df t.ratio p.value
## a - b -0.05712 0.0529 85 -1.079 0.5297
## a - c -0.06641 0.0529 85 -1.254 0.4248
## b - c -0.00929 0.0529 85 -0.176 0.9832
##
## Results are averaged over the levels of: period
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmip(gadw.m1977, ~veg.water, CIs = TRUE, type = "response",
linearg = list(linetype = NA),
CIarg = list(lwd = 1, alpha = 1))## NOTE: Results may be misleading due to involvement in interactions
#1978----------
gadw.m1978 <- lmer(sqrt(ip.ha) ~ veg.water * treatment * period + (1|week), data = gadw1978)## boundary (singular) fit: see ?isSingularplot(gadw.m1978,
main = "GADW 1978 Fit",
xlab = "Residuals",
ylab = "Fitted")
anova(gadw.m1978)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## veg.water 1.31869 0.65935 2 72 12.7929 1.762e-05 ***
## treatment 0.06291 0.06291 1 72 1.2207 0.2729
## period 0.12216 0.06108 2 72 1.1851 0.3116
## veg.water:treatment 0.18602 0.09301 2 72 1.8046 0.1719
## veg.water:period 0.06257 0.01564 4 72 0.3035 0.8747
## treatment:period 0.19586 0.09793 2 72 1.9001 0.1570
## veg.water:treatment:period 0.01943 0.00486 4 72 0.0942 0.9840
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(gadw.m1978, pairwise~veg.water, type = "response")## NOTE: Results may be misleading due to involvement in interactions
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## veg.water response SE df lower.CL upper.CL
## a 0.358 0.0496 29.1 0.263 0.466
## b 0.606 0.0646 29.1 0.482 0.746
## c 0.235 0.0402 29.1 0.160 0.324
##
## Results are averaged over the levels of: treatment, period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## a - b -0.181 0.0586 68 -3.084 0.0082
## a - c 0.113 0.0586 68 1.930 0.1380
## b - c 0.294 0.0586 68 5.014 <.0001
##
## Results are averaged over the levels of: treatment, period
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmip(gadw.m1978, ~veg.water, CIs = TRUE, type = "response",
linearg = list(linetype = NA),
CIarg = list(lwd = 1, alpha = 1))## NOTE: Results may be misleading due to involvement in interactions
3.3 Mallards
In 1977, we found mallard breeding pairs to be affected by veg:water but it depended on whether the experimental plots were mowed or tilled. We also found that the effect of treatment on mallard IBPs depended on survey time. However, in 1978 mallard IBPs were clearly positively affected by 50:50 veg-to-water interspersion.
#1977----------
mall.m1977 <- lmer(sqrt(ip.ha) ~ veg.water * treatment * period + (1|week), data = mall1977)
#that fit isn't great but it'll probably do
plot(mall.m1977,
main = "MALL 1977 Fit",
xlab = "Residuals",
ylab = "Fitted")
anova(mall.m1977)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## veg.water 4.4461 2.22303 2 85 17.4456 4.486e-07 ***
## treatment 0.0800 0.08000 1 85 0.6278 0.430367
## period 0.2210 0.11052 2 85 0.8673 0.423754
## veg.water:treatment 1.2568 0.62838 2 85 4.9313 0.009414 **
## veg.water:period 0.6035 0.15086 4 85 1.1839 0.323753
## treatment:period 0.9742 0.48710 2 85 3.8226 0.025723 *
## veg.water:treatment:period 0.1998 0.04996 4 85 0.3920 0.813806
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(mall.m1977, pairwise~veg.water, type = "response")## NOTE: Results may be misleading due to involvement in interactions## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## veg.water response SE df lower.CL upper.CL
## a 1.17 0.140 23.9 0.896 1.47
## b 2.33 0.197 23.9 1.940 2.76
## c 1.24 0.144 23.9 0.961 1.56
##
## Results are averaged over the levels of: treatment, period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## a - b -0.4462 0.0841 85 -5.303 <.0001
## a - c -0.0336 0.0841 85 -0.399 0.9161
## b - c 0.4126 0.0841 85 4.904 <.0001
##
## Results are averaged over the levels of: treatment, period
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmeans(mall.m1977, pairwise~veg.water:treatment, type = "response")## NOTE: Results may be misleading due to involvement in interactions
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## veg.water treatment response SE df lower.CL upper.CL
## a mow 1.320 0.202 55.5 0.947 1.76
## b mow 2.413 0.273 55.5 1.897 2.99
## c mow 0.876 0.165 55.5 0.577 1.24
## a till 1.022 0.178 55.5 0.697 1.41
## b till 2.247 0.264 55.5 1.750 2.81
## c till 1.668 0.227 55.5 1.244 2.15
##
## Results are averaged over the levels of: period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## a mow - b mow -0.4043 0.119 85 -3.398 0.0128
## a mow - c mow 0.2133 0.119 85 1.793 0.4760
## a mow - a till 0.1381 0.119 85 1.160 0.8542
## a mow - b till -0.3500 0.119 85 -2.942 0.0467
## a mow - c till -0.1424 0.119 85 -1.197 0.8373
## b mow - c mow 0.6177 0.119 85 5.191 <.0001
## b mow - a till 0.5424 0.119 85 4.558 0.0002
## b mow - b till 0.0543 0.119 85 0.457 0.9974
## b mow - c till 0.2620 0.119 85 2.202 0.2480
## c mow - a till -0.0752 0.119 85 -0.632 0.9883
## c mow - b till -0.5633 0.119 85 -4.734 0.0001
## c mow - c till -0.3557 0.119 85 -2.989 0.0412
## a till - b till -0.4881 0.119 85 -4.102 0.0013
## a till - c till -0.2805 0.119 85 -2.357 0.1835
## b till - c till 0.2076 0.119 85 1.745 0.5065
##
## Results are averaged over the levels of: period
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 6 estimatesemmeans(mall.m1977, pairwise~treatment:period, type = "response")## NOTE: Results may be misleading due to involvement in interactions
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## treatment period response SE df lower.CL upper.CL
## mow afternoon 1.44 0.211 55.5 1.047 1.89
## till afternoon 1.61 0.223 55.5 1.195 2.09
## mow evening 1.34 0.204 55.5 0.965 1.78
## till evening 2.07 0.253 55.5 1.591 2.60
## mow morning 1.64 0.225 55.5 1.220 2.12
## till morning 1.20 0.192 55.5 0.844 1.62
##
## Results are averaged over the levels of: veg.water
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## mow afternoon - till afternoon -0.0698 0.119 85 -0.586 0.9917
## mow afternoon - mow evening 0.0407 0.119 85 0.342 0.9994
## mow afternoon - till evening -0.2383 0.119 85 -2.003 0.3495
## mow afternoon - mow morning -0.0812 0.119 85 -0.682 0.9835
## mow afternoon - till morning 0.1043 0.119 85 0.877 0.9511
## till afternoon - mow evening 0.1105 0.119 85 0.929 0.9381
## till afternoon - till evening -0.1686 0.119 85 -1.417 0.7170
## till afternoon - mow morning -0.0114 0.119 85 -0.096 1.0000
## till afternoon - till morning 0.1741 0.119 85 1.463 0.6885
## mow evening - till evening -0.2790 0.119 85 -2.345 0.1879
## mow evening - mow morning -0.1219 0.119 85 -1.025 0.9085
## mow evening - till morning 0.0636 0.119 85 0.534 0.9946
## till evening - mow morning 0.1571 0.119 85 1.321 0.7729
## till evening - till morning 0.3426 0.119 85 2.879 0.0549
## mow morning - till morning 0.1855 0.119 85 1.559 0.6275
##
## Results are averaged over the levels of: veg.water
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 6 estimatesemmip(mall.m1977, veg.water~treatment, CIs = TRUE, type = "response",
linearg = list(linetype = "dashed"),
CIarg = list(lwd = 1, alpha = 1))## NOTE: Results may be misleading due to involvement in interactions
emmip(mall.m1977, period~treatment, CIs = TRUE,
linearg = list(linetype = "dashed"),
CIarg = list(lwd = 1, alpha = 1))## NOTE: Results may be misleading due to involvement in interactions
#1978----------
mall.m1978 <- lmer(sqrt(ip.ha) ~ veg.water * treatment * period + (1|week), data = mall1978)
plot(mall.m1978,
main = "MALL 1978 Fit",
xlab = "Residuals",
ylab = "Fitted")
anova(mall.m1978)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## veg.water 1.09989 0.54994 2 68 12.3887 2.584e-05 ***
## treatment 0.01073 0.01073 1 68 0.2418 0.6245
## period 0.04048 0.02024 2 68 0.4560 0.6358
## veg.water:treatment 0.14377 0.07188 2 68 1.6194 0.2056
## veg.water:period 0.09960 0.02490 4 68 0.5609 0.6918
## treatment:period 0.06457 0.03229 2 68 0.7273 0.4869
## veg.water:treatment:period 0.02367 0.00592 4 68 0.1333 0.9696
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(mall.m1978, pairwise~veg.water, type = "response")## NOTE: Results may be misleading due to involvement in interactions
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## veg.water response SE df lower.CL upper.CL
## a 0.501 0.0784 8.58 0.338 0.696
## b 0.862 0.1028 8.58 0.644 1.112
## c 0.465 0.0755 8.58 0.309 0.653
##
## Results are averaged over the levels of: treatment, period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## a - b -0.2207 0.0544 68 -4.056 0.0004
## a - c 0.0256 0.0544 68 0.470 0.8854
## b - c 0.2463 0.0544 68 4.527 0.0001
##
## Results are averaged over the levels of: treatment, period
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmip(mall.m1978, ~veg.water, CIs = TRUE, type = "response",
linearg = list(linetype = NA),
CIarg = list(lwd = 1, alpha = 1))## NOTE: Results may be misleading due to involvement in interactions
3.4 Northern Pintail
In 1977 we detected no evidence of any factors affecting northern pintail densities. However, the model resulted in singular fit which suggests possible overfitting where fixed and random variables explained nearly all of the variation in NOPI IBPs. In particular, we did not detect any affect with our repeated measures ANOVA similar to Kaminski and Prince (1981); thus, the addition of “week” as the random variable may have overcomplicated the model and is now fitting the noise in the data. Therefore, our model may be suspect. Nevertheless, overfitting is generally a problem when the investigator worries about Type I errors. But we reached the same conclusion as Kaminski and Prince (1981), so we chose not to fit the model again using regularization terms in a Bayesian context.
In 1978 we detected an affect of treatment, in which northern pintail densities were greatest in mowed plots.
#1977----------
nopi.m1977 <- lmer(sqrt(ip.ha) ~ veg.water * treatment * period + (1|week), data = nopi1977)## boundary (singular) fit: see ?isSingularnopi.m1977.1 <- lmer(sqrt(ip.ha) ~ veg.water * treatment + (1|week), data = nopi1977)## boundary (singular) fit: see ?isSingularnopi.m1977.2 <- lm(sqrt(ip.ha) ~ veg.water * treatment * period , data = nopi1977)
#boundary is singlular meaning the model overfit (not generalizing well)
plot(nopi.m1977,
main = "NOPI 1977 Fit",
xlab = "Residuals",
ylab = "Fitted")
#1978----------
nopi.m1978 <- lmer(sqrt(ip.ha) ~ veg.water * treatment * period + (1|week), data = nopi1978)
plot(nopi.m1978,
main = "NOPI 1978 Fit",
xlab = "Residuals",
ylab = "Fitted")
anova(nopi.m1978)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## veg.water 0.14817 0.07409 2 68 1.8723 0.161617
## treatment 0.33147 0.33147 1 68 8.3767 0.005104 **
## period 0.11665 0.05832 2 68 1.4739 0.236246
## veg.water:treatment 0.09605 0.04802 2 68 1.2136 0.303475
## veg.water:period 0.07497 0.01874 4 68 0.4737 0.754897
## treatment:period 0.02026 0.01013 2 68 0.2561 0.774837
## veg.water:treatment:period 0.15141 0.03785 4 68 0.9566 0.437062
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(nopi.m1978, pairwise~treatment, type = "response")## NOTE: Results may be misleading due to involvement in interactions## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## treatment response SE df lower.CL upper.CL
## mow 0.371 0.0425 9.61 0.282 0.473
## till 0.238 0.0341 9.61 0.168 0.321
##
## Results are averaged over the levels of: veg.water, period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## mow - till 0.121 0.0419 68 2.894 0.0051
##
## Results are averaged over the levels of: veg.water, period
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-rogeremmip(nopi.m1978, ~treatment, CIs = TRUE, type = "response",
linearg = list(linetype = NA),
CIarg = list(lwd = 1, alpha = 1))## NOTE: Results may be misleading due to involvement in interactions
## 3.5 Northern Shoveler
We detected effects of veg.water ratio, treatment, and their interaction in 1977 for northern shovelers. Interestingly, tilled plots resulted in greatest IBP densities of NOSH likely due to their feeding ecology and preference for aquatic invertebrates. We did not detect any affect on northern shoveler density in 1978
#1977----------
nosh.m1977 <- lmer(sqrt(ip.ha) ~ veg.water * treatment * period + (1|week), data = nosh1977)
#that fit isn't great but it'll probably do
plot(nosh.m1977,
main = "NOSH 1977 Fit",
xlab = "Residuals",
ylab = "Fitted")
anova(nosh.m1977)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## veg.water 0.85729 0.42864 2 85 8.2870 0.0005153 ***
## treatment 1.54843 1.54843 1 85 29.9358 4.424e-07 ***
## period 0.17740 0.08870 2 85 1.7148 0.1861653
## veg.water:treatment 0.43479 0.21740 2 85 4.2029 0.0181719 *
## veg.water:period 0.08854 0.02214 4 85 0.4279 0.7880758
## treatment:period 0.02066 0.01033 2 85 0.1997 0.8193658
## veg.water:treatment:period 0.13302 0.03325 4 85 0.6429 0.6333727
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(nosh.m1977, pairwise~veg.water, type = "response")## NOTE: Results may be misleading due to involvement in interactions## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## veg.water response SE df lower.CL upper.CL
## a 0.473 0.0855 8.78 0.298 0.687
## b 0.630 0.0987 8.78 0.426 0.874
## c 0.331 0.0716 8.78 0.189 0.514
##
## Results are averaged over the levels of: treatment, period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## a - b -0.106 0.0536 85 -1.986 0.1220
## a - c 0.112 0.0536 85 2.085 0.0990
## b - c 0.218 0.0536 85 4.071 0.0003
##
## Results are averaged over the levels of: treatment, period
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmeans(nosh.m1977, pairwise~treatment, type = "response")## NOTE: Results may be misleading due to involvement in interactions
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## treatment response SE df lower.CL upper.CL
## mow 0.320 0.0658 6.77 0.183 0.496
## till 0.649 0.0937 6.77 0.445 0.891
##
## Results are averaged over the levels of: veg.water, period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## mow - till -0.239 0.0438 85 -5.471 <.0001
##
## Results are averaged over the levels of: veg.water, period
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-rogeremmeans(nosh.m1977, pairwise~veg.water|treatment, type = "response")## NOTE: Results may be misleading due to involvement in interactions
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## treatment = mow:
## veg.water response SE df lower.CL upper.CL
## a 0.248 0.0725 16 0.1183 0.426
## b 0.575 0.1104 16 0.3645 0.832
## c 0.195 0.0643 16 0.0824 0.355
##
## treatment = till:
## veg.water response SE df lower.CL upper.CL
## a 0.768 0.1276 16 0.5217 1.063
## b 0.688 0.1208 16 0.4561 0.968
## c 0.504 0.1034 16 0.3086 0.747
##
## Results are averaged over the levels of: period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## treatment = mow:
## contrast estimate SE df t.ratio p.value
## a - b -0.2598 0.0758 85 -3.427 0.0027
## a - c 0.0568 0.0758 85 0.749 0.7350
## b - c 0.3166 0.0758 85 4.176 0.0002
##
## treatment = till:
## contrast estimate SE df t.ratio p.value
## a - b 0.0469 0.0758 85 0.619 0.8100
## a - c 0.1668 0.0758 85 2.200 0.0770
## b - c 0.1198 0.0758 85 1.580 0.2596
##
## Results are averaged over the levels of: period
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmip(nosh.m1977, veg.water~treatment, CIs = TRUE, type = "response",
linearg = list(linetype = "dashed"),
CIarg = list(lwd = 1, alpha = 1))## NOTE: Results may be misleading due to involvement in interactions
#1978----------
nosh.m1978 <- lmer(sqrt(ip.ha) ~ veg.water * treatment * period + (1|week), data = nosh1978)
plot(nosh.m1978,
main = "NOSH 1978 Fit",
xlab = "Residuals",
ylab = "Fitted")
anova(nosh.m1978)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## veg.water 0.34697 0.173487 2 68 2.8893 0.06247 .
## treatment 0.08027 0.080270 1 68 1.3368 0.25164
## period 0.10182 0.050912 2 68 0.8479 0.43280
## veg.water:treatment 0.22116 0.110580 2 68 1.8416 0.16638
## veg.water:period 0.03219 0.008046 4 68 0.1340 0.96932
## treatment:period 0.10696 0.053479 2 68 0.8906 0.41513
## veg.water:treatment:period 0.23624 0.059060 4 68 0.9836 0.42244
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(nosh.m1978, pairwise~veg.water, type = "response")## NOTE: Results may be misleading due to involvement in interactions
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## veg.water response SE df lower.CL upper.CL
## a 0.491 0.0664 21.6 0.363 0.639
## b 0.543 0.0698 21.6 0.407 0.697
## c 0.349 0.0560 21.6 0.242 0.475
##
## Results are averaged over the levels of: treatment, period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## a - b -0.0359 0.0633 68 -0.568 0.8378
## a - c 0.1100 0.0633 68 1.739 0.1982
## b - c 0.1459 0.0633 68 2.307 0.0616
##
## Results are averaged over the levels of: treatment, period
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmip(nosh.m1978, ~veg.water, CIs = TRUE, type = "response",
linearg = list(linetype = NA),
CIarg = list(lwd = 1, alpha = 1))## NOTE: Results may be misleading due to involvement in interactions
4. Means Testing with diversity as response variable
We also analyze dabbling duck diversity in response to experimentally manipulated marshes to stay consistent with Kaminski and Prince (1981). We maintain fixed and random variables as in previous analyses, but now we take the Broullion’s Diversity Index for dabbling duck species (which are weighted by there densities).
## [1] "treatment" "veg.water" "year" "week" "period" "bwte"
## [7] "gadw" "mall" "nopi" "nosh"4.1 Calculating Broullion’s Diversity Index
Kaminksi and Prince (1981) calculated dabbling duck diversity using Brillouin’s Diversity Index because samples were convenient where MERP cells were located and thus, non-random. Therefore, we also calculated Brillouin’s Index similarly as the natural logarithm as factorial of the frequency of duck observations minus the summed total and divided by n.
#write a function that will calculate Broullion's Index
Brillouin_Index <- function(x){
N <- sum(x)
(log10(factorial(N)) - sum(log10(factorial(x)))) / N
}
dat2$bindex <- apply(dat2[6:10], 1, Brillouin_Index)
rcompanion::plotNormalHistogram(dat2$bindex,
xlab = "Brillouin's Diversity Index")
dat21977 <- dat2 %>%
filter(year == "1977")
dat21978 <- dat2 %>%
filter(year == "1978")4.2 Diversity Relative to Experimental Marsh Manipulation
We found evidence that all veg:water ratio, treatment, and their interaction influenced dabbling duck diversity in 1977. Further, observed dabbler diversity by treatment plots depended on the timing of the survey in 1977. In 1978, we found support for effects of veg:water interspersion, treatment, and survey timing but not for any interaction.
#1977----------
bindex1977 <- lmer(bindex ~ veg.water * treatment * period + (1|week), data = dat21977)
#that fit isn't great but it'll probably do
plot(bindex1977,
main = "Diversity 1977 Fit",
xlab = "Residuals",
ylab = "Fitted")
anova(bindex1977)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## veg.water 0.0280085 0.0140043 2 85 5.6881 0.004804 **
## treatment 0.0217664 0.0217664 1 85 8.8408 0.003832 **
## period 0.0084128 0.0042064 2 85 1.7085 0.187300
## veg.water:treatment 0.0254725 0.0127363 2 85 5.1730 0.007585 **
## veg.water:period 0.0023198 0.0005799 4 85 0.2356 0.917579
## treatment:period 0.0193539 0.0096770 2 85 3.9305 0.023302 *
## veg.water:treatment:period 0.0025189 0.0006297 4 85 0.2558 0.905387
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(bindex1977, pairwise~veg.water, type = "response")## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## veg.water emmean SE df lower.CL upper.CL
## a 0.278 0.0188 6.59 0.233 0.323
## b 0.308 0.0188 6.59 0.263 0.353
## c 0.271 0.0188 6.59 0.226 0.316
##
## Results are averaged over the levels of: treatment, period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## a - b -0.0305 0.0117 85 -2.609 0.0287
## a - c 0.0064 0.0117 85 0.547 0.8482
## b - c 0.0369 0.0117 85 3.156 0.0062
##
## Results are averaged over the levels of: treatment, period
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmeans(bindex1977, pairwise~treatment, type = "response")## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## treatment emmean SE df lower.CL upper.CL
## mow 0.272 0.0182 5.77 0.227 0.317
## till 0.300 0.0182 5.77 0.255 0.345
##
## Results are averaged over the levels of: veg.water, period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## mow - till -0.0284 0.00955 85 -2.973 0.0038
##
## Results are averaged over the levels of: veg.water, period
## Degrees-of-freedom method: kenward-rogeremmeans(bindex1977, pairwise~veg.water | treatment, type = "response")## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## treatment = mow:
## veg.water emmean SE df lower.CL upper.CL
## a 0.257 0.0205 9.32 0.211 0.303
## b 0.315 0.0205 9.32 0.269 0.362
## c 0.243 0.0205 9.32 0.197 0.289
##
## treatment = till:
## veg.water emmean SE df lower.CL upper.CL
## a 0.299 0.0205 9.32 0.253 0.345
## b 0.301 0.0205 9.32 0.255 0.347
## c 0.300 0.0205 9.32 0.254 0.346
##
## Results are averaged over the levels of: period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## treatment = mow:
## contrast estimate SE df t.ratio p.value
## a - b -0.05845 0.0165 85 -3.534 0.0019
## a - c 0.01424 0.0165 85 0.861 0.6662
## b - c 0.07269 0.0165 85 4.395 0.0001
##
## treatment = till:
## contrast estimate SE df t.ratio p.value
## a - b -0.00257 0.0165 85 -0.155 0.9868
## a - c -0.00144 0.0165 85 -0.087 0.9958
## b - c 0.00113 0.0165 85 0.068 0.9975
##
## Results are averaged over the levels of: period
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmeans(bindex1977, pairwise~period | treatment, type = "response")## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## treatment = mow:
## period emmean SE df lower.CL upper.CL
## afternoon 0.265 0.0205 9.32 0.219 0.312
## evening 0.260 0.0205 9.32 0.214 0.306
## morning 0.290 0.0205 9.32 0.244 0.336
##
## treatment = till:
## period emmean SE df lower.CL upper.CL
## afternoon 0.281 0.0205 9.32 0.235 0.328
## evening 0.325 0.0205 9.32 0.279 0.371
## morning 0.293 0.0205 9.32 0.247 0.340
##
## Results are averaged over the levels of: veg.water
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## treatment = mow:
## contrast estimate SE df t.ratio p.value
## afternoon - evening 0.00568 0.0165 85 0.344 0.9370
## afternoon - morning -0.02442 0.0165 85 -1.476 0.3074
## evening - morning -0.03010 0.0165 85 -1.820 0.1693
##
## treatment = till:
## contrast estimate SE df t.ratio p.value
## afternoon - evening -0.04396 0.0165 85 -2.658 0.0252
## afternoon - morning -0.01213 0.0165 85 -0.734 0.7443
## evening - morning 0.03183 0.0165 85 1.925 0.1380
##
## Results are averaged over the levels of: veg.water
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmip(bindex1977, veg.water~treatment, CIs = TRUE,
linearg = list(linetype = "dashed"),
CIarg = list(lwd = 1, alpha = 1))## NOTE: Results may be misleading due to involvement in interactions
emmip(bindex1977, period~treatment, CIs = TRUE,
linearg = list(linetype = "dashed"),
CIarg = list(lwd = 1, alpha = 1))## NOTE: Results may be misleading due to involvement in interactions
#1978----------
bindex1978 <- lmer(bindex ~ veg.water * treatment * period + (1|week), data = dat21978)
plot(bindex1978,
main = "Diversity 1978 Fit",
xlab = "Residuals",
ylab = "Fitted")
anova(bindex1978)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## veg.water 0.093656 0.046828 2 68 18.2908 4.404e-07 ***
## treatment 0.012150 0.012150 1 68 4.7457 0.03284 *
## period 0.015696 0.007848 2 68 3.0654 0.05313 .
## veg.water:treatment 0.008249 0.004124 2 68 1.6110 0.20722
## veg.water:period 0.003779 0.000945 4 68 0.3690 0.82991
## treatment:period 0.008549 0.004274 2 68 1.6695 0.19596
## veg.water:treatment:period 0.006790 0.001698 4 68 0.6630 0.61982
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(bindex1978, pairwise~veg.water, type = "response")## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## veg.water emmean SE df lower.CL upper.CL
## a 0.280 0.0108 14.6 0.257 0.303
## b 0.326 0.0108 14.6 0.303 0.349
## c 0.248 0.0108 14.6 0.225 0.271
##
## Results are averaged over the levels of: treatment, period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## a - b -0.0462 0.0131 68 -3.535 0.0021
## a - c 0.0324 0.0131 68 2.483 0.0406
## b - c 0.0786 0.0131 68 6.018 <.0001
##
## Results are averaged over the levels of: treatment, period
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmip(bindex1978, ~veg.water | treatment | period, CIs = TRUE,
linearg = list(linetype = NA),
CIarg = list(lwd = 1, alpha = 1))
emmip(bindex1978, ~veg.water, CIs = TRUE,
linearg = list(linetype = NA),
CIarg = list(lwd = 1, alpha = 1))## NOTE: Results may be misleading due to involvement in interactions
5. Re-creating graphs from Kaminski and Prince (1981)
vegwater <- read.csv("vegwater.csv")
vegwater1977 <- vegwater %>%
filter(year == "1977")
vegwater1978 <- vegwater %>%
filter(year == "1978")
vw1 <- ggplot(vegwater1977, aes(x = species, y = ibp, fill = veg.water)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(x="Dabbling Duck Species", y ="Indicated Breeding Pairs per Hectare") +
geom_errorbar(aes(ymin = lower.cl, ymax = upper.cl), width = 0.2,
position = position_dodge(0.9))+
theme_bw() +
theme(legend.position = "none") + ggtitle("Summer 1977")
vw2 <- ggplot(vegwater1978, aes(x = species, y = ibp, fill = veg.water)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(x="Dabbling Duck Species", y ="") +
geom_errorbar(aes(ymin = lower.cl, ymax = upper.cl), width = 0.2,
position = position_dodge(0.9))+ ylim(0, 4) +
theme_bw()+
ggtitle("Summer 1978")
ggpubr::ggarrange(vw1, vw2, ncol=2, nrow=1, align="v", legend = "bottom") 