I think I've managed to reverse-engineer your original data frame from the cross-tabs:
df
#> Zone Year total_count empty_count
#> 1 Crocodile 2016 1 0
#> 2 Rankin 2016 14 3
#> 3 West 2016 29 7
#> 4 Whipray 2016 9 4
#> 5 Crocodile 2017 18 8
#> 6 Rankin 2017 46 17
#> 7 West 2017 67 31
#> 8 Whipray 2017 7 4
#> 9 Crocodile 2018 7 2
#> 10 Rankin 2018 69 8
#> 11 West 2018 58 17
#> 12 Whipray 2018 15 0
It seems to me that rather than attempting multiple pairwise comparisons you could carry out a single logistic regression to find out where the significant differences are. Just ensure you have a "non-empty" count column, and that your years are factors:
df$non_empty_count <- df$total_count - df$empty_count
df$Year <- as.factor(df$Year)
And your logistic regression looks like this:
model <- glm(cbind(empty_count, non_empty_count) ~ Zone + Year,
data = df, family = binomial)
summary(model)
#>
#> Call:
#> glm(formula = cbind(empty_count, non_empty_count) ~ Zone + Year,
#> family = binomial, data = df)
#>
#> Deviance Residuals:
#> Min 1Q Median 3Q Max
#> -2.47563 -0.59548 0.04512 0.56564 1.26509
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -0.9588 0.5360 -1.789 0.0737 .
#> ZoneRankin -0.4864 0.4740 -1.026 0.3048
#> ZoneWest 0.1281 0.4552 0.281 0.7784
#> ZoneWhipray -0.1395 0.6038 -0.231 0.8173
#> Year2017 0.7975 0.3659 2.180 0.0293 *
#> Year2018 -0.3861 0.3831 -1.008 0.3136
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> Null deviance: 39.255 on 11 degrees of freedom
#> Residual deviance: 11.800 on 6 degrees of freedom
#> AIC: 57.905
#>
#> Number of Fisher Scoring iterations: 5
Which you can interpret as saying that while the proportion of empty stomachs didn't vary significantly between sites, there was a significantly higher proportion of fish with empty stomachs in 2017 compared to other years across all sites.
Reproducible data
df <- structure(list(Zone = structure(c(1L, 2L, 3L, 4L, 1L, 2L, 3L,
4L, 1L, 2L, 3L, 4L), .Label = c("Crocodile", "Rankin", "West",
"Whipray"), class = "factor"), Year = c(2016, 2016, 2016, 2016,
2017, 2017, 2017, 2017, 2018, 2018, 2018, 2018), total_count = c(1L,
14L, 29L, 9L, 18L, 46L, 67L, 7L, 7L, 69L, 58L, 15L), empty_count = c(0L,
3L, 7L, 4L, 8L, 17L, 31L, 4L, 2L, 8L, 17L, 0L)), row.names = c(NA,
-12L), class = "data.frame")
Created on 2022-02-09 by the reprex package (v2.0.1)