Q1Model_RMarkdown

####Packages

library(MASS)
library(lme4)

Loading required package: Matrix

library(car)

Loading required package: carData

library(multcomp)

Loading required package: mvtnorm

Loading required package: survival

Loading required package: TH.data

Attaching package: 'TH.data'

The following object is masked from 'package:MASS':

    geyser

library(emmeans)

Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'

library(broom)
library(effects)

lattice theme set by effectsTheme()
See ?effectsTheme for details.

library(broom.mixed)
library(glmmTMB)
library(DHARMa)

This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')

subsite_table <- readRDS("/Users/emmabacon/Documents/EB_UrbanTrees/output/diversitymetrics_table_rare.rds")

#TREE DENSITY MODEL OPTIONS ##OPTION 1: Glmer, with negative binomial(theta = 1). Chose this orignally because this is what Kayleigh had done, but there are differences in random effect when compared to glmer.nb below

#OPTION 1
#runnning a linear mixed model to test the effect of green space type on tree abundance, including logged sample area as a covariate because abudance is not rarefied 
modelgst_den0 <- glmer(tree_density_plant_ha ~ GreenSpace + (1|Plot), data= subsite_table, family = "negative.binomial"(theta = 1))

boundary (singular) fit: see help('isSingular')

print(modelgst_den0) #Random effect is zero

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Negative Binomial(1)  ( log )
Formula: tree_density_plant_ha ~ GreenSpace + (1 | Plot)
   Data: subsite_table
      AIC       BIC    logLik -2*log(L)  df.resid 
1302.5191 1323.4401 -643.2596 1286.5191        93 
Random effects:
 Groups Name        Std.Dev.
 Plot   (Intercept) 0       
Number of obs: 101, groups:  Plot, 22
Fixed Effects:
                  (Intercept)        GreenSpaceInstitutional  
                      5.11818                       -0.20067  
               GreenSpacePark  GreenSpacePublic Right-Of-Way  
                     -0.05607                        0.98157  
        GreenSpaceResidential           GreenSpaceVacant Lot  
                      0.24508                        0.24802  
optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 

#boundary (singular) fit: see help('isSingular') - the random effect variance is near-zero
VarCorr(modelgst_den0) #near-zero, so random effect is not contributing meaningfully

 Groups Name        Std.Dev.
 Plot   (Intercept) 0       

lme4::isSingular(modelgst_den0, tol = 1e-4) 

[1] TRUE

#Random effect is singular, so remove

#Diagnostics
plot(modelgst_den0, which = 1)

shapiro.test(resid(modelgst_den0))

    Shapiro-Wilk normality test

data:  resid(modelgst_den0)
W = 0.98574, p-value = 0.3519

hist(residuals(modelgst_den0))

qqnorm(residuals(modelgst_den0))
qqline(residuals(modelgst_den0))

#Just ran this to compare with below
#check effect of green space type on abundance using model outputs
Anova(modelgst_den0)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: tree_density_plant_ha
            Chisq Df Pr(>Chisq)   
GreenSpace 17.614  5   0.003472 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Green space statistically significant portion of the variation
tbl.den0 <- broom::tidy(modelgst_den0)
tbl.den0

# A tibble: 7 × 7
  effect   group term                     estimate std.error statistic   p.value
  <chr>    <chr> <chr>                       <dbl>     <dbl>     <dbl>     <dbl>
1 fixed    <NA>  (Intercept)                5.12       0.259    19.8    6.06e-87
2 fixed    <NA>  GreenSpaceInstitutional   -0.201      0.351    -0.572  5.67e- 1
3 fixed    <NA>  GreenSpacePark            -0.0561     0.373    -0.150  8.80e- 1
4 fixed    <NA>  GreenSpacePublic Right-…   0.982      0.336     2.93   3.44e- 3
5 fixed    <NA>  GreenSpaceResidential      0.245      0.336     0.730  4.65e- 1
6 fixed    <NA>  GreenSpaceVacant Lot       0.248      0.409     0.606  5.45e- 1
7 ran_pars Plot  sd__(Intercept)            0         NA        NA     NA       

#tukey test to check pairwise comparisons of abundances of each green space type
Pairwise.den0 <- emmeans(modelgst_den0, pairwise ~ GreenSpace, type = "response")
#summarizing pairwise differences across green space types
summary(Pairwise.den0) 

$emmeans
 GreenSpace          response   SE  df asymp.LCL asymp.UCL
 Commercial               167 43.3 Inf     100.5       277
 Institutional            137 32.3 Inf      86.0       217
 Park                     158 42.3 Inf      93.4       267
 Public Right-Of-Way      446 95.1 Inf     293.4       677
 Residential              213 45.6 Inf     140.4       324
 Vacant Lot               214 67.8 Inf     115.0       398

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                              ratio     SE  df null z.ratio p.value
 Commercial / Institutional            1.222 0.4290 Inf    1   0.572  0.9928
 Commercial / Park                     1.058 0.3940 Inf    1   0.150  1.0000
 Commercial / (Public Right-Of-Way)    0.375 0.1260 Inf    1  -2.925  0.0403
 Commercial / Residential              0.783 0.2630 Inf    1  -0.730  0.9783
 Commercial / Vacant Lot               0.780 0.3190 Inf    1  -0.606  0.9906
 Institutional / Park                  0.865 0.3090 Inf    1  -0.404  0.9986
 Institutional / (Public Right-Of-Way) 0.307 0.0977 Inf    1  -3.711  0.0028
 Institutional / Residential           0.640 0.2040 Inf    1  -1.398  0.7282
 Institutional / Vacant Lot            0.638 0.2530 Inf    1  -1.134  0.8671
 Park / (Public Right-Of-Way)          0.354 0.1210 Inf    1  -3.028  0.0297
 Park / Residential                    0.740 0.2540 Inf    1  -0.878  0.9518
 Park / Vacant Lot                     0.738 0.3060 Inf    1  -0.732  0.9780
 (Public Right-Of-Way) / Residential   2.089 0.6310 Inf    1   2.438  0.1431
 (Public Right-Of-Way) / Vacant Lot    2.082 0.7960 Inf    1   1.920  0.3900
 Residential / Vacant Lot              0.997 0.3810 Inf    1  -0.008  1.0000

P value adjustment: tukey method for comparing a family of 6 estimates 
Tests are performed on the log scale 

##OPTION 2: Glmer.nb, where theta is set by package. I do not get the singularity warnings with this model.

modelgst_den <- glmer.nb(tree_density_plant_ha ~ GreenSpace + (1|Plot), data= subsite_table)
summary(modelgst_den)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Negative Binomial(2.323)  ( log )
Formula: tree_density_plant_ha ~ GreenSpace + (1 | Plot)
   Data: subsite_table

      AIC       BIC    logLik -2*log(L)  df.resid 
   1272.5    1293.4    -628.2    1256.5        93 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4101 -0.7865 -0.0450  0.3388  4.7168 

Random effects:
 Groups Name        Variance Std.Dev.
 Plot   (Intercept) 0.01896  0.1377  
Number of obs: 101, groups:  Plot, 22

Fixed effects:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    5.10781    0.17624  28.982  < 2e-16 ***
GreenSpaceInstitutional       -0.20234    0.23337  -0.867    0.386    
GreenSpacePark                -0.03856    0.25057  -0.154    0.878    
GreenSpacePublic Right-Of-Way  0.96269    0.22855   4.212 2.53e-05 ***
GreenSpaceResidential          0.24938    0.22378   1.114    0.265    
GreenSpaceVacant Lot           0.25663    0.27365   0.938    0.348    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) GrnSpI GrnSpP GSPR-O GrnSpR
GrnSpcInstt -0.725                            
GreenSpcPrk -0.689  0.508                     
GrnSpPR-O-W -0.743  0.566  0.510              
GrnSpcRsdnt -0.764  0.571  0.537  0.592       
GrnSpcVcntL -0.623  0.466  0.441  0.474  0.487

#Diagnostics
plot(modelgst_den, which = 1)

shapiro.test(resid(modelgst_den))

    Shapiro-Wilk normality test

data:  resid(modelgst_den)
W = 0.9882, p-value = 0.5159

hist(residuals(modelgst_den))

qqnorm(residuals(modelgst_den))
qqline(residuals(modelgst_den))

#check effect of green space type on tree density using model outputs
Anova(modelgst_den)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: tree_density_plant_ha
            Chisq Df Pr(>Chisq)    
GreenSpace 35.703  5  1.089e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Green space statistically significant portion of the variation
tbl.den <- broom::tidy(modelgst_den)
tbl.den

# A tibble: 7 × 7
  effect   group term                    estimate std.error statistic    p.value
  <chr>    <chr> <chr>                      <dbl>     <dbl>     <dbl>      <dbl>
1 fixed    <NA>  (Intercept)               5.11       0.176    29.0    1.11e-184
2 fixed    <NA>  GreenSpaceInstitutional  -0.202      0.233    -0.867  3.86e-  1
3 fixed    <NA>  GreenSpacePark           -0.0386     0.251    -0.154  8.78e-  1
4 fixed    <NA>  GreenSpacePublic Right…   0.963      0.229     4.21   2.53e-  5
5 fixed    <NA>  GreenSpaceResidential     0.249      0.224     1.11   2.65e-  1
6 fixed    <NA>  GreenSpaceVacant Lot      0.257      0.274     0.938  3.48e-  1
7 ran_pars Plot  sd__(Intercept)           0.138     NA        NA     NA        

#tukey test to check pairwise comparisons of abundances of each green space type
Pairwise.den <- emmeans(modelgst_den, pairwise ~ GreenSpace, type = "response")
#summarizing pairwise differences across green space types
summary(Pairwise.den) 

$emmeans
 GreenSpace          response   SE  df asymp.LCL asymp.UCL
 Commercial               165 29.1 Inf     117.0       234
 Institutional            135 21.7 Inf      98.5       185
 Park                     159 28.9 Inf     111.4       227
 Public Right-Of-Way      433 66.3 Inf     320.7       584
 Residential              212 30.6 Inf     159.8       282
 Vacant Lot               214 45.8 Inf     140.4       325

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                              ratio     SE  df null z.ratio p.value
 Commercial / Institutional            1.224 0.2860 Inf    1   0.867  0.9543
 Commercial / Park                     1.039 0.2600 Inf    1   0.154  1.0000
 Commercial / (Public Right-Of-Way)    0.382 0.0873 Inf    1  -4.212  0.0004
 Commercial / Residential              0.779 0.1740 Inf    1  -1.114  0.8756
 Commercial / Vacant Lot               0.774 0.2120 Inf    1  -0.938  0.9368
 Institutional / Park                  0.849 0.2040 Inf    1  -0.681  0.9841
 Institutional / (Public Right-Of-Way) 0.312 0.0671 Inf    1  -5.414  <.0001
 Institutional / Residential           0.637 0.1350 Inf    1  -2.132  0.2707
 Institutional / Vacant Lot            0.632 0.1670 Inf    1  -1.737  0.5073
 Park / (Public Right-Of-Way)          0.367 0.0874 Inf    1  -4.210  0.0004
 Park / Residential                    0.750 0.1720 Inf    1  -1.255  0.8096
 Park / Vacant Lot                     0.744 0.2070 Inf    1  -1.062  0.8963
 (Public Right-Of-Way) / Residential   2.041 0.4170 Inf    1   3.492  0.0064
 (Public Right-Of-Way) / Vacant Lot    2.026 0.5280 Inf    1   2.710  0.0733
 Residential / Vacant Lot              0.993 0.2540 Inf    1  -0.028  1.0000

P value adjustment: tukey method for comparing a family of 6 estimates 
Tests are performed on the log scale 

##OPTION 3: glm.nb - Finally, here is the model without random effecs, which I was using before, but now that glmer.nb does show random effects, I am not sure. I also get the warning because this is negative binomial on non-count data.

modelgst_den1 <- glm.nb(tree_density_plant_ha ~ GreenSpace, data = subsite_table)
summary(modelgst_den1)

Call:
glm.nb(formula = tree_density_plant_ha ~ GreenSpace, data = subsite_table, 
    init.theta = 2.237817151, link = log)

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    5.11818    0.17375  29.457  < 2e-16 ***
GreenSpaceInstitutional       -0.20067    0.23542  -0.852    0.394    
GreenSpacePark                -0.05607    0.25012  -0.224    0.823    
GreenSpacePublic Right-Of-Way  0.98157    0.22495   4.363 1.28e-05 ***
GreenSpaceResidential          0.24508    0.22520   1.088    0.276    
GreenSpaceVacant Lot           0.24802    0.27449   0.904    0.366    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(2.2378) family taken to be 1)

    Null deviance: 151.48  on 100  degrees of freedom
Residual deviance: 108.17  on  95  degrees of freedom
AIC: 1270.7

Number of Fisher Scoring iterations: 1

              Theta:  2.238 
          Std. Err.:  0.299 

 2 x log-likelihood:  -1256.717 

# MODEL DIAGNOSTICS
plot(modelgst_den1, which = 1)

shapiro.test(resid(modelgst_den1))

    Shapiro-Wilk normality test

data:  resid(modelgst_den1)
W = 0.98557, p-value = 0.3421

hist(residuals(modelgst_den1))

qqnorm(residuals(modelgst_den1))
qqline(residuals(modelgst_den1))

#When comparing with the glmer.nb, this looks a bit better? But similar, so let me know what you think

#check effect of green space type on tree density using model outputs
Anova(modelgst_den1)

Analysis of Deviance Table (Type II tests)

Response: tree_density_plant_ha
           LR Chisq Df Pr(>Chisq)    
GreenSpace   43.311  5  3.196e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Green space statistically significant portion of the variation
tbl.den1 <- broom::tidy(modelgst_den1)
tbl.den1

# A tibble: 6 × 5
  term                          estimate std.error statistic   p.value
  <chr>                            <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)                     5.12       0.174    29.5   1.03e-190
2 GreenSpaceInstitutional        -0.201      0.235    -0.852 3.94e-  1
3 GreenSpacePark                 -0.0561     0.250    -0.224 8.23e-  1
4 GreenSpacePublic Right-Of-Way   0.982      0.225     4.36  1.28e-  5
5 GreenSpaceResidential           0.245      0.225     1.09  2.76e-  1
6 GreenSpaceVacant Lot            0.248      0.274     0.904 3.66e-  1

#tukey test to check pairwise comparisons of abundances of each green space type
Pairwise.den1 <- emmeans(modelgst_den1, pairwise ~ GreenSpace, type = "response")
#summarizing pairwise differences across green space types
summary(Pairwise.den1) 

$emmeans
 GreenSpace          response   SE  df asymp.LCL asymp.UCL
 Commercial               167 29.0 Inf       119       235
 Institutional            137 21.7 Inf       100       187
 Park                     158 28.4 Inf       111       225
 Public Right-Of-Way      446 63.7 Inf       337       590
 Residential              213 30.6 Inf       161       283
 Vacant Lot               214 45.5 Inf       141       325

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                              ratio     SE  df null z.ratio p.value
 Commercial / Institutional            1.222 0.2880 Inf    1   0.852  0.9575
 Commercial / Park                     1.058 0.2650 Inf    1   0.224  0.9999
 Commercial / (Public Right-Of-Way)    0.375 0.0843 Inf    1  -4.363  0.0002
 Commercial / Residential              0.783 0.1760 Inf    1  -1.088  0.8863
 Commercial / Vacant Lot               0.780 0.2140 Inf    1  -0.904  0.9457
 Institutional / Park                  0.865 0.2080 Inf    1  -0.602  0.9909
 Institutional / (Public Right-Of-Way) 0.307 0.0655 Inf    1  -5.534  <.0001
 Institutional / Residential           0.640 0.1370 Inf    1  -2.084  0.2957
 Institutional / Vacant Lot            0.638 0.1690 Inf    1  -1.691  0.5376
 Park / (Public Right-Of-Way)          0.354 0.0814 Inf    1  -4.516  0.0001
 Park / Residential                    0.740 0.1700 Inf    1  -1.309  0.7800
 Park / Vacant Lot                     0.738 0.2050 Inf    1  -1.092  0.8847
 (Public Right-Of-Way) / Residential   2.089 0.4230 Inf    1   3.640  0.0037
 (Public Right-Of-Way) / Vacant Lot    2.082 0.5330 Inf    1   2.865  0.0479
 Residential / Vacant Lot              0.997 0.2560 Inf    1  -0.011  1.0000

P value adjustment: tukey method for comparing a family of 6 estimates 
Tests are performed on the log scale 

#In conclusion, for tree density, I am leaning towards glmer.nb. Anova/pairwise results are similar for all, although glm.nb is the most significant

#TREE ABUNDANCE MODEL OPTIONS ##OPTION 1: Glmer, with negative binomial(theta = 1). Chose this orignally because this is what Kayleigh had done, but there are differences in random effect when compared to glmer.nb below

modelgst_abu0 <- glmer(tree_count ~ GreenSpace + log(area_ha) + (1|Plot), data= subsite_table, family = "negative.binomial"(theta = 1))

boundary (singular) fit: see help('isSingular')

print(modelgst_abu0)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Negative Binomial(1)  ( log )
Formula: tree_count ~ GreenSpace + log(area_ha) + (1 | Plot)
   Data: subsite_table
      AIC       BIC    logLik -2*log(L)  df.resid 
1186.7657 1210.3018 -584.3829 1168.7657        92 
Random effects:
 Groups Name        Std.Dev.
 Plot   (Intercept) 1.43e-05
Number of obs: 101, groups:  Plot, 22
Fixed Effects:
                  (Intercept)        GreenSpaceInstitutional  
                       3.7386                         0.4249  
               GreenSpacePark  GreenSpacePublic Right-Of-Way  
                       0.8049                         1.0076  
        GreenSpaceResidential           GreenSpaceVacant Lot  
                       1.4172                         0.8651  
                 log(area_ha)  
                       0.8226  
optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 

#boundary (singular) fit: see help('isSingular') - the random effect variance is near-zero
VarCorr(modelgst_abu0) #near-zero, so random effect is not contributing meaningfully

 Groups Name        Std.Dev.  
 Plot   (Intercept) 1.4298e-05

lme4::isSingular(modelgst_abu0, tol = 1e-4) #Random effect is singular

[1] TRUE

#Diagnostics
plot(modelgst_abu0, which = 1)

shapiro.test(resid(modelgst_abu0))

    Shapiro-Wilk normality test

data:  resid(modelgst_abu0)
W = 0.98318, p-value = 0.2275

hist(residuals(modelgst_abu0))

qqnorm(residuals(modelgst_abu0))
qqline(residuals(modelgst_abu0))

#check effect of green space type on abundance using model outputs
Anova(modelgst_abu0)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: tree_count
              Chisq Df Pr(>Chisq)    
GreenSpace   13.152  5    0.02199 *  
log(area_ha) 45.393  1  1.612e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Green space statistically significant portion of the variation
tbl.abu0 <- broom::tidy(modelgst_abu0)
tbl.abu0

# A tibble: 8 × 7
  effect   group term                     estimate std.error statistic   p.value
  <chr>    <chr> <chr>                       <dbl>     <dbl>     <dbl>     <dbl>
1 fixed    <NA>  (Intercept)               3.74e+0     0.302     12.4   3.73e-35
2 fixed    <NA>  GreenSpaceInstitutional   4.25e-1     0.394      1.08  2.81e- 1
3 fixed    <NA>  GreenSpacePark            8.05e-1     0.391      2.06  3.95e- 2
4 fixed    <NA>  GreenSpacePublic Right-…  1.01e+0     0.454      2.22  2.63e- 2
5 fixed    <NA>  GreenSpaceResidential     1.42e+0     0.496      2.86  4.26e- 3
6 fixed    <NA>  GreenSpaceVacant Lot      8.65e-1     0.428      2.02  4.33e- 2
7 fixed    <NA>  log(area_ha)              8.23e-1     0.122      6.74  1.61e-11
8 ran_pars Plot  sd__(Intercept)           1.43e-5    NA         NA    NA       

#tukey test to check pairwise comparisons of abundances of each green space type
Pairwise.abu0 <- emmeans(modelgst_abu0, pairwise ~ GreenSpace, type = "response")
#summarizing pairwise differences across green space types
summary(Pairwise.abu0) 

$emmeans
 GreenSpace          response    SE  df asymp.LCL asymp.UCL
 Commercial              94.6  35.3 Inf      45.5       196
 Institutional          144.6  37.3 Inf      87.3       240
 Park                   211.5  67.7 Inf     112.9       396
 Public Right-Of-Way    259.0  56.1 Inf     169.3       396
 Residential            390.1  91.7 Inf     246.0       619
 Vacant Lot             224.6 105.0 Inf      90.1       560

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                              ratio    SE  df null z.ratio p.value
 Commercial / Institutional            0.654 0.258 Inf    1  -1.077  0.8905
 Commercial / Park                     0.447 0.175 Inf    1  -2.058  0.3094
 Commercial / (Public Right-Of-Way)    0.365 0.166 Inf    1  -2.222  0.2276
 Commercial / Residential              0.242 0.120 Inf    1  -2.858  0.0489
 Commercial / Vacant Lot               0.421 0.180 Inf    1  -2.021  0.3301
 Institutional / Park                  0.684 0.252 Inf    1  -1.031  0.9076
 Institutional / (Public Right-Of-Way) 0.558 0.194 Inf    1  -1.677  0.5471
 Institutional / Residential           0.371 0.139 Inf    1  -2.643  0.0872
 Institutional / Vacant Lot            0.644 0.301 Inf    1  -0.940  0.9360
 Park / (Public Right-Of-Way)          0.817 0.329 Inf    1  -0.503  0.9961
 Park / Residential                    0.542 0.237 Inf    1  -1.398  0.7282
 Park / Vacant Lot                     0.942 0.427 Inf    1  -0.133  1.0000
 (Public Right-Of-Way) / Residential   0.664 0.205 Inf    1  -1.330  0.7686
 (Public Right-Of-Way) / Vacant Lot    1.153 0.620 Inf    1   0.265  0.9998
 Residential / Vacant Lot              1.737 1.010 Inf    1   0.949  0.9338

P value adjustment: tukey method for comparing a family of 6 estimates 
Tests are performed on the log scale 

##OPTION 2: Glmer.nb, where theta is set by package. I do not get the singularity warnings with this model

modelgst_abu <- glmer.nb(tree_count ~ GreenSpace + (1|Plot), data = subsite_table)
print(modelgst_abu)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Negative Binomial(1.3138)  ( log )
Formula: tree_count ~ GreenSpace + (1 | Plot)
   Data: subsite_table
      AIC       BIC    logLik -2*log(L)  df.resid 
 1219.348  1240.269  -601.674  1203.348        93 
Random effects:
 Groups Name        Std.Dev.
 Plot   (Intercept) 0.1591  
Number of obs: 101, groups:  Plot, 22
Fixed Effects:
                  (Intercept)        GreenSpaceInstitutional  
                        3.051                          1.566  
               GreenSpacePark  GreenSpacePublic Right-Of-Way  
                        1.550                          2.767  
        GreenSpaceResidential           GreenSpaceVacant Lot  
                        3.585                          0.203  

summary(modelgst_abu)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Negative Binomial(1.3138)  ( log )
Formula: tree_count ~ GreenSpace + (1 | Plot)
   Data: subsite_table

      AIC       BIC    logLik -2*log(L)  df.resid 
   1219.3    1240.3    -601.7    1203.3        93 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.1402 -0.7221 -0.1814  0.3948  3.3922 

Random effects:
 Groups Name        Variance Std.Dev.
 Plot   (Intercept) 0.02531  0.1591  
Number of obs: 101, groups:  Plot, 22

Fixed effects:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     3.0508     0.2448  12.461  < 2e-16 ***
GreenSpaceInstitutional         1.5660     0.3150   4.971 6.65e-07 ***
GreenSpacePark                  1.5502     0.3438   4.509 6.50e-06 ***
GreenSpacePublic Right-Of-Way   2.7670     0.3055   9.057  < 2e-16 ***
GreenSpaceResidential           3.5846     0.3020  11.871  < 2e-16 ***
GreenSpaceVacant Lot            0.2030     0.3716   0.546    0.585    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) GrnSpI GrnSpP GSPR-O GrnSpR
GrnSpcInstt -0.695                            
GreenSpcPrk -0.710  0.487                     
GrnSpPR-O-W -0.783  0.556  0.563              
GrnSpcRsdnt -0.768  0.571  0.548  0.617       
GrnSpcVcntL -0.589  0.463  0.418  0.473  0.485

#Error - potentially due to adding log(area_ha)? When I remove area, the model runs without warnings.
#I want to stay consistent with other models, so what do you suggest?

# Abundance MODEL DIAGNOSTICS
plot(modelgst_abu, which = 1)

shapiro.test(resid(modelgst_abu))

    Shapiro-Wilk normality test

data:  resid(modelgst_abu)
W = 0.99151, p-value = 0.7794

hist(residuals(modelgst_abu))

qqnorm(residuals(modelgst_abu))
qqline(residuals(modelgst_abu))

#check effect of green space type on abundance using model outputs
Anova(modelgst_abu)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: tree_count
            Chisq Df Pr(>Chisq)    
GreenSpace 199.16  5  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Green space statistically significant portion of the variation
tbl.abu <- broom::tidy(modelgst_abu)
tbl.abu

# A tibble: 7 × 7
  effect   group term                     estimate std.error statistic   p.value
  <chr>    <chr> <chr>                       <dbl>     <dbl>     <dbl>     <dbl>
1 fixed    <NA>  (Intercept)                 3.05      0.245    12.5    1.22e-35
2 fixed    <NA>  GreenSpaceInstitutional     1.57      0.315     4.97   6.65e- 7
3 fixed    <NA>  GreenSpacePark              1.55      0.344     4.51   6.50e- 6
4 fixed    <NA>  GreenSpacePublic Right-…    2.77      0.305     9.06   1.34e-19
5 fixed    <NA>  GreenSpaceResidential       3.58      0.302    11.9    1.68e-32
6 fixed    <NA>  GreenSpaceVacant Lot        0.203     0.372     0.546  5.85e- 1
7 ran_pars Plot  sd__(Intercept)             0.159    NA        NA     NA       

#tukey test to check pairwise comparisons of abundances of each green space type
Pairwise.abu <- emmeans(modelgst_abu, pairwise ~ GreenSpace, type = "response")
#summarizing pairwise differences across green space types
summary(Pairwise.abu) 

$emmeans
 GreenSpace          response     SE  df asymp.LCL asymp.UCL
 Commercial              21.1   5.17 Inf      13.1      34.1
 Institutional          101.2  23.10 Inf      64.7     158.2
 Park                    99.6  24.10 Inf      62.0     160.0
 Public Right-Of-Way    336.2  63.90 Inf     231.7     488.0
 Residential            761.6 148.00 Inf     520.9    1113.5
 Vacant Lot              25.9   7.81 Inf      14.3      46.7

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                                ratio       SE  df null z.ratio
 Commercial / Institutional             0.2089  0.06580 Inf    1  -4.971
 Commercial / Park                      0.2122  0.07300 Inf    1  -4.509
 Commercial / (Public Right-Of-Way)     0.0629  0.01920 Inf    1  -9.057
 Commercial / Residential               0.0277  0.00838 Inf    1 -11.871
 Commercial / Vacant Lot                0.8163  0.30300 Inf    1  -0.546
 Institutional / Park                   1.0159  0.34000 Inf    1   0.047
 Institutional / (Public Right-Of-Way)  0.3009  0.08800 Inf    1  -4.108
 Institutional / Residential            0.1328  0.03800 Inf    1  -7.056
 Institutional / Vacant Lot             3.9079  1.40000 Inf    1   3.797
 Park / (Public Right-Of-Way)           0.2962  0.09040 Inf    1  -3.985
 Park / Residential                     0.1308  0.04040 Inf    1  -6.580
 Park / Vacant Lot                      3.8467  1.49000 Inf    1   3.485
 (Public Right-Of-Way) / Residential    0.4415  0.11700 Inf    1  -3.074
 (Public Right-Of-Way) / Vacant Lot    12.9874  4.58000 Inf    1   7.279
 Residential / Vacant Lot              29.4188 10.20000 Inf    1   9.742
 p.value
  <.0001
  0.0001
  <.0001
  <.0001
  0.9942
  1.0000
  0.0006
  <.0001
  0.0020
  0.0010
  <.0001
  0.0065
  0.0258
  <.0001
  <.0001

P value adjustment: tukey method for comparing a family of 6 estimates 
Tests are performed on the log scale 

##OPTION 3: glm.nb - Finally, here is the model without random effecs, which I was using before, but now that glmer.nb does show random effects, I am not sure

modelgst_abu1 <- glm.nb(tree_count ~ GreenSpace + log(area_ha), data = subsite_table)
summary(modelgst_abu1)

Call:
glm.nb(formula = tree_count ~ GreenSpace + log(area_ha), data = subsite_table, 
    init.theta = 2.211662881, link = log)

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    3.71208    0.18812  19.732  < 2e-16 ***
GreenSpaceInstitutional        0.46022    0.24957   1.844 0.065179 .  
GreenSpacePark                 0.82848    0.26084   3.176 0.001492 ** 
GreenSpacePublic Right-Of-Way  1.03564    0.27231   3.803 0.000143 ***
GreenSpaceResidential          1.44621    0.28783   5.025 5.05e-07 ***
GreenSpaceVacant Lot           0.89869    0.31329   2.869 0.004123 ** 
log(area_ha)                   0.82140    0.07792  10.541  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(2.2117) family taken to be 1)

    Null deviance: 483.90  on 100  degrees of freedom
Residual deviance: 108.21  on  94  degrees of freedom
AIC: 1157.6

Number of Fisher Scoring iterations: 1

              Theta:  2.212 
          Std. Err.:  0.311 

 2 x log-likelihood:  -1141.618 

# Abundance MODEL DIAGNOSTICS
plot(modelgst_abu1, which = 1)

shapiro.test(resid(modelgst_abu1))

    Shapiro-Wilk normality test

data:  resid(modelgst_abu1)
W = 0.98413, p-value = 0.2681

hist(residuals(modelgst_abu1))

qqnorm(residuals(modelgst_abu1))
qqline(residuals(modelgst_abu1))

#check effect of green space type on abundance using model outputs
Anova(modelgst_abu1)

Analysis of Deviance Table (Type II tests)

Response: tree_count
             LR Chisq Df Pr(>Chisq)    
GreenSpace     26.999  5  5.707e-05 ***
log(area_ha)   81.762  1  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Green space statistically significant portion of the variation
tbl.abu <- broom::tidy(modelgst_abu1)
tbl.abu

# A tibble: 7 × 5
  term                          estimate std.error statistic  p.value
  <chr>                            <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)                      3.71     0.188      19.7  1.14e-86
2 GreenSpaceInstitutional          0.460    0.250       1.84 6.52e- 2
3 GreenSpacePark                   0.828    0.261       3.18 1.49e- 3
4 GreenSpacePublic Right-Of-Way    1.04     0.272       3.80 1.43e- 4
5 GreenSpaceResidential            1.45     0.288       5.02 5.05e- 7
6 GreenSpaceVacant Lot             0.899    0.313       2.87 4.12e- 3
7 log(area_ha)                     0.821    0.0779     10.5  5.59e-26

#tukey test to check pairwise comparisons of abundances of each green space type
Pairwise.abu1 <- emmeans(modelgst_abu1, pairwise ~ GreenSpace, type = "response")
#summarizing pairwise differences across green space types
summary(Pairwise.abu1) 

$emmeans
 GreenSpace          response   SE  df asymp.LCL asymp.UCL
 Commercial                92 20.0 Inf      60.1       141
 Institutional            146 26.5 Inf     102.1       208
 Park                     211 43.8 Inf     140.1       317
 Public Right-Of-Way      259 37.8 Inf     194.7       345
 Residential              391 59.7 Inf     289.5       527
 Vacant Lot               226 72.9 Inf     120.0       425

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                              ratio     SE  df null z.ratio p.value
 Commercial / Institutional            0.631 0.1580 Inf    1  -1.844  0.4373
 Commercial / Park                     0.437 0.1140 Inf    1  -3.176  0.0187
 Commercial / (Public Right-Of-Way)    0.355 0.0967 Inf    1  -3.803  0.0020
 Commercial / Residential              0.235 0.0678 Inf    1  -5.025  <.0001
 Commercial / Vacant Lot               0.407 0.1280 Inf    1  -2.869  0.0474
 Institutional / Park                  0.692 0.1710 Inf    1  -1.489  0.6713
 Institutional / (Public Right-Of-Way) 0.562 0.1350 Inf    1  -2.389  0.1599
 Institutional / Residential           0.373 0.0948 Inf    1  -3.882  0.0015
 Institutional / Vacant Lot            0.645 0.2050 Inf    1  -1.377  0.7412
 Park / (Public Right-Of-Way)          0.813 0.2140 Inf    1  -0.788  0.9697
 Park / Residential                    0.539 0.1490 Inf    1  -2.229  0.2243
 Park / Vacant Lot                     0.932 0.2990 Inf    1  -0.219  0.9999
 (Public Right-Of-Way) / Residential   0.663 0.1360 Inf    1  -2.000  0.3419
 (Public Right-Of-Way) / Vacant Lot    1.147 0.4240 Inf    1   0.371  0.9991
 Residential / Vacant Lot              1.729 0.6730 Inf    1   1.406  0.7235

P value adjustment: tukey method for comparing a family of 6 estimates 
Tests are performed on the log scale 

#NOTE: Pairwise comparisons are very different for glmer vs. glmer.nb vs glm.nb, I am very confused

#SPECIES RICHNESS MODEL OPTIONS ##Since species richness is rarefied, it is not count data. Am I still able to do negative binomial? ##OPTION 1: with negative binomial(theta = 1) ##do not need to include area as a covariate because we already calculated rarefied diversity metrics based on sample completeness

modelgst_sr0 <- glmer(Cov_rich ~ GreenSpace + (1|Plot), data= subsite_table, family = "negative.binomial"(theta = 1))

boundary (singular) fit: see help('isSingular')

print(modelgst_sr0)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Negative Binomial(1)  ( log )
Formula: Cov_rich ~ GreenSpace + (1 | Plot)
   Data: subsite_table
      AIC       BIC    logLik -2*log(L)  df.resid 
 875.7234  896.6444 -429.8617  859.7234        93 
Random effects:
 Groups Name        Std.Dev.
 Plot   (Intercept) 0       
Number of obs: 101, groups:  Plot, 22
Fixed Effects:
                  (Intercept)        GreenSpaceInstitutional  
                       2.0492                         1.0454  
               GreenSpacePark  GreenSpacePublic Right-Of-Way  
                       0.9426                         1.9757  
        GreenSpaceResidential           GreenSpaceVacant Lot  
                       2.0732                        -0.1994  
optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 

#boundary (singular) fit: see help('isSingular') - the random effect variance is near-zero
VarCorr(modelgst_sr0) #near-zero, so random effect is not contributing meaningfully

 Groups Name        Std.Dev.
 Plot   (Intercept) 0       

lme4::isSingular(modelgst_sr0, tol = 1e-4) #Random effect is singular

[1] TRUE

plot(modelgst_sr0, which = 1)

shapiro.test(resid(modelgst_sr0))

    Shapiro-Wilk normality test

data:  resid(modelgst_sr0)
W = 0.97559, p-value = 0.05772

hist(residuals(modelgst_sr0))

qqnorm(residuals(modelgst_sr0))
qqline(residuals(modelgst_sr0))

Anova(modelgst_sr0)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: Cov_rich
            Chisq Df Pr(>Chisq)    
GreenSpace 66.838  5  4.656e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise.sr0 <- emmeans(modelgst_sr0, pairwise ~ GreenSpace, type = "response")
summary(Pairwise.sr0) 

$emmeans
 GreenSpace          response    SE  df asymp.LCL asymp.UCL
 Commercial              7.76  2.13 Inf      4.53      13.3
 Institutional          22.08  5.32 Inf     13.77      35.4
 Park                   19.92  5.46 Inf     11.65      34.1
 Public Right-Of-Way    55.98 12.00 Inf     36.72      85.3
 Residential            61.71 13.30 Inf     40.49      94.0
 Vacant Lot              6.36  2.16 Inf      3.26      12.4

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                              ratio     SE  df null z.ratio p.value
 Commercial / Institutional            0.352 0.1280 Inf    1  -2.863  0.0481
 Commercial / Park                     0.390 0.1510 Inf    1  -2.432  0.1453
 Commercial / (Public Right-Of-Way)    0.139 0.0483 Inf    1  -5.668  <.0001
 Commercial / Residential              0.126 0.0438 Inf    1  -5.949  <.0001
 Commercial / Vacant Lot               1.221 0.5330 Inf    1   0.456  0.9975
 Institutional / Park                  1.108 0.4040 Inf    1   0.282  0.9998
 Institutional / (Public Right-Of-Way) 0.394 0.1270 Inf    1  -2.880  0.0459
 Institutional / Residential           0.358 0.1160 Inf    1  -3.183  0.0183
 Institutional / Vacant Lot            3.472 1.4500 Inf    1   2.986  0.0337
 Park / (Public Right-Of-Way)          0.356 0.1240 Inf    1  -2.967  0.0357
 Park / Residential                    0.323 0.1120 Inf    1  -3.247  0.0148
 Park / Vacant Lot                     3.133 1.3700 Inf    1   2.615  0.0936
 (Public Right-Of-Way) / Residential   0.907 0.2760 Inf    1  -0.321  0.9996
 (Public Right-Of-Way) / Vacant Lot    8.803 3.5400 Inf    1   5.404  <.0001
 Residential / Vacant Lot              9.704 3.9000 Inf    1   5.648  <.0001

P value adjustment: tukey method for comparing a family of 6 estimates 
Tests are performed on the log scale 

##OPTION 2: Glmer.nb, where theta is set by package. I do not get the singularity warnings with this model

modelgst_sr <- glmer.nb(Cov_rich ~ GreenSpace + (1|Plot), data= subsite_table)
print(modelgst_sr)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Negative Binomial(3.4047)  ( log )
Formula: Cov_rich ~ GreenSpace + (1 | Plot)
   Data: subsite_table
      AIC       BIC    logLik -2*log(L)  df.resid 
 833.6377  854.5587 -408.8189  817.6377        93 
Random effects:
 Groups Name        Std.Dev.
 Plot   (Intercept) 0.1809  
Number of obs: 101, groups:  Plot, 22
Fixed Effects:
                  (Intercept)        GreenSpaceInstitutional  
                       2.0401                         1.0498  
               GreenSpacePark  GreenSpacePublic Right-Of-Way  
                       0.9156                         1.9639  
        GreenSpaceResidential           GreenSpaceVacant Lot  
                       2.0573                        -0.1605  

#Warnings due to non count data

# SR MODEL DIAGNOSTICS
plot(modelgst_sr, which = 1)

shapiro.test(resid(modelgst_sr))

    Shapiro-Wilk normality test

data:  resid(modelgst_sr)
W = 0.97567, p-value = 0.05857

hist(residuals(modelgst_sr))

qqnorm(residuals(modelgst_sr))
qqline(residuals(modelgst_sr))

Anova(modelgst_sr)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: Cov_rich
            Chisq Df Pr(>Chisq)    
GreenSpace 178.13  5  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise.sr <- emmeans(modelgst_sr, pairwise ~ GreenSpace, type = "response")
summary(Pairwise.sr) 

$emmeans
 GreenSpace          response   SE  df asymp.LCL asymp.UCL
 Commercial              7.69 1.35 Inf      5.46      10.8
 Institutional          21.98 3.18 Inf     16.54      29.2
 Park                   19.22 3.20 Inf     13.87      26.6
 Public Right-Of-Way    54.82 6.92 Inf     42.80      70.2
 Residential            60.18 7.63 Inf     46.95      77.1
 Vacant Lot              6.55 1.45 Inf      4.25      10.1

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                              ratio     SE  df null z.ratio p.value
 Commercial / Institutional            0.350 0.0771 Inf    1  -4.768  <.0001
 Commercial / Park                     0.400 0.0939 Inf    1  -3.902  0.0013
 Commercial / (Public Right-Of-Way)    0.140 0.0293 Inf    1  -9.389  <.0001
 Commercial / Residential              0.128 0.0267 Inf    1  -9.845  <.0001
 Commercial / Vacant Lot               1.174 0.3230 Inf    1   0.583  0.9922
 Institutional / Park                  1.144 0.2460 Inf    1   0.624  0.9893
 Institutional / (Public Right-Of-Way) 0.401 0.0740 Inf    1  -4.951  <.0001
 Institutional / Residential           0.365 0.0677 Inf    1  -5.434  <.0001
 Institutional / Vacant Lot            3.355 0.8710 Inf    1   4.659  <.0001
 Park / (Public Right-Of-Way)          0.351 0.0702 Inf    1  -5.232  <.0001
 Park / Residential                    0.319 0.0638 Inf    1  -5.711  <.0001
 Park / Vacant Lot                     2.933 0.8050 Inf    1   3.922  0.0012
 (Public Right-Of-Way) / Residential   0.911 0.1540 Inf    1  -0.553  0.9939
 (Public Right-Of-Way) / Vacant Lot    8.368 2.0900 Inf    1   8.493  <.0001
 Residential / Vacant Lot              9.187 2.3100 Inf    1   8.837  <.0001

P value adjustment: tukey method for comparing a family of 6 estimates 
Tests are performed on the log scale 

##OPTION 3: glm.nb - Finally, here is the model without random effecs, which I was using before

modelgst_sr1 <- glm.nb(Cov_rich ~ GreenSpace, data = subsite_table)
print(modelgst_sr1)

Call:  glm.nb(formula = Cov_rich ~ GreenSpace, data = subsite_table, 
    init.theta = 3.069647171, link = log)

Coefficients:
                  (Intercept)        GreenSpaceInstitutional  
                       2.0492                         1.0454  
               GreenSpacePark  GreenSpacePublic Right-Of-Way  
                       0.9426                         1.9757  
        GreenSpaceResidential           GreenSpaceVacant Lot  
                       2.0732                        -0.1994  

Degrees of Freedom: 100 Total (i.e. Null);  95 Residual
Null Deviance:      274.5 
Residual Deviance: 110.2    AIC: 832.7

summary(modelgst_sr1)

Call:
glm.nb(formula = Cov_rich ~ GreenSpace, data = subsite_table, 
    init.theta = 3.069647171, link = log)

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     2.0492     0.1741  11.771  < 2e-16 ***
GreenSpaceInstitutional         1.0454     0.2257   4.633 3.61e-06 ***
GreenSpacePark                  0.9426     0.2391   3.943 8.06e-05 ***
GreenSpacePublic Right-Of-Way   1.9757     0.2143   9.219  < 2e-16 ***
GreenSpaceResidential           2.0732     0.2141   9.682  < 2e-16 ***
GreenSpaceVacant Lot           -0.1994     0.2804  -0.711    0.477    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(3.0696) family taken to be 1)

    Null deviance: 274.48  on 100  degrees of freedom
Residual deviance: 110.24  on  95  degrees of freedom
AIC: 832.7

Number of Fisher Scoring iterations: 1

              Theta:  3.070 
          Std. Err.:  0.510 

 2 x log-likelihood:  -818.699 

confint(modelgst_sr1, level = 0.95)

                                   2.5 %    97.5 %
(Intercept)                    1.7155475 2.3997402
GreenSpaceInstitutional        0.6012506 1.4876548
GreenSpacePark                 0.4743896 1.4137443
GreenSpacePublic Right-Of-Way  1.5522864 2.3941008
GreenSpaceResidential          1.6500930 2.4912218
GreenSpaceVacant Lot          -0.7468153 0.3549816

# SR MODEL DIAGNOSTICS
plot(modelgst_sr1, which = 1)

shapiro.test(resid(modelgst_sr1))

    Shapiro-Wilk normality test

data:  resid(modelgst_sr1)
W = 0.97941, p-value = 0.1157

hist(residuals(modelgst_sr1))

qqnorm(residuals(modelgst_sr1))
qqline(residuals(modelgst_sr1))

Anova(modelgst_sr1)

Analysis of Deviance Table (Type II tests)

Response: Cov_rich
           LR Chisq Df Pr(>Chisq)    
GreenSpace   164.24  5  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise.sr1 <- emmeans(modelgst_sr1, pairwise ~ GreenSpace, type = "response")
summary(Pairwise.sr1) 

$emmeans
 GreenSpace          response   SE  df asymp.LCL asymp.UCL
 Commercial              7.76 1.35 Inf      5.52     10.92
 Institutional          22.08 3.17 Inf     16.66     29.25
 Park                   19.92 3.26 Inf     14.45     27.47
 Public Right-Of-Way    55.98 7.00 Inf     43.81     71.51
 Residential            61.71 7.69 Inf     48.33     78.79
 Vacant Lot              6.36 1.40 Inf      4.13      9.78

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                              ratio     SE  df null z.ratio p.value
 Commercial / Institutional            0.352 0.0793 Inf    1  -4.633  0.0001
 Commercial / Park                     0.390 0.0932 Inf    1  -3.943  0.0011
 Commercial / (Public Right-Of-Way)    0.139 0.0297 Inf    1  -9.219  <.0001
 Commercial / Residential              0.126 0.0269 Inf    1  -9.682  <.0001
 Commercial / Vacant Lot               1.221 0.3420 Inf    1   0.711  0.9807
 Institutional / Park                  1.108 0.2410 Inf    1   0.472  0.9971
 Institutional / (Public Right-Of-Way) 0.394 0.0751 Inf    1  -4.887  <.0001
 Institutional / Residential           0.358 0.0680 Inf    1  -5.405  <.0001
 Institutional / Vacant Lot            3.472 0.9120 Inf    1   4.742  <.0001
 Park / (Public Right-Of-Way)          0.356 0.0733 Inf    1  -5.013  <.0001
 Park / Residential                    0.323 0.0665 Inf    1  -5.491  <.0001
 Park / Vacant Lot                     3.133 0.8590 Inf    1   4.165  0.0004
 (Public Right-Of-Way) / Residential   0.907 0.1600 Inf    1  -0.552  0.9939
 (Public Right-Of-Way) / Vacant Lot    8.803 2.2300 Inf    1   8.603  <.0001
 Residential / Vacant Lot              9.704 2.4500 Inf    1   8.994  <.0001

P value adjustment: tukey method for comparing a family of 6 estimates 
Tests are performed on the log scale 

#OPTION 4: Gamma GLMM with log link, random effects are back and not zero

modelgst_sr_gamma <- glmer(Cov_rich ~ GreenSpace + (1 | Plot), data = subsite_table, family = Gamma(link = "log"))
summary(modelgst_sr_gamma) 

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Gamma  ( log )
Formula: Cov_rich ~ GreenSpace + (1 | Plot)
   Data: subsite_table

      AIC       BIC    logLik -2*log(L)  df.resid 
    836.5     857.4    -410.2     820.5        93 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6862 -0.5642 -0.1629  0.4564  3.3659 

Random effects:
 Groups   Name        Variance Std.Dev.
 Plot     (Intercept) 0.04961  0.2227  
 Residual             0.33535  0.5791  
Number of obs: 101, groups:  Plot, 22

Fixed effects:
                              Estimate Std. Error t value Pr(>|z|)    
(Intercept)                     2.0038     0.1784  11.230  < 2e-16 ***
GreenSpaceInstitutional         1.0927     0.2229   4.901 9.52e-07 ***
GreenSpacePark                  0.9239     0.2363   3.911 9.21e-05 ***
GreenSpacePublic Right-Of-Way   1.9896     0.2123   9.372  < 2e-16 ***
GreenSpaceResidential           2.0745     0.2115   9.810  < 2e-16 ***
GreenSpaceVacant Lot           -0.1268     0.2593  -0.489    0.625    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) GrnSpI GrnSpP GSPR-O GrnSpR
GrnSpcInstt -0.684                            
GreenSpcPrk -0.639  0.497                     
GrnSpPR-O-W -0.727  0.576  0.542              
GrnSpcRsdnt -0.722  0.565  0.543  0.617       
GrnSpcVcntL -0.569  0.443  0.412  0.482  0.475

#Residuals
resid_gamma <- residuals(modelgst_sr_gamma, type = "pearson")
qqnorm(resid_gamma); qqline(resid_gamma)

hist(resid_gamma)

plot(resid_gamma)

Anova(modelgst_sr_gamma)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: Cov_rich
            Chisq Df Pr(>Chisq)    
GreenSpace 172.64  5  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise.srgamma <- emmeans(modelgst_sr_gamma, pairwise ~ GreenSpace, type = "response")
summary(Pairwise.srgamma) 

$emmeans
 GreenSpace          response   SE  df asymp.LCL asymp.UCL
 Commercial              7.42 1.32 Inf      5.23     10.52
 Institutional          22.12 3.64 Inf     16.02     30.55
 Park                   18.68 3.43 Inf     13.03     26.78
 Public Right-Of-Way    54.24 8.01 Inf     40.61     72.44
 Residential            59.05 8.78 Inf     44.12     79.02
 Vacant Lot              6.53 1.41 Inf      4.28      9.97

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                              ratio     SE  df null z.ratio p.value
 Commercial / Institutional            0.335 0.0748 Inf    1  -4.901  <.0001
 Commercial / Park                     0.397 0.0938 Inf    1  -3.911  0.0013
 Commercial / (Public Right-Of-Way)    0.137 0.0290 Inf    1  -9.372  <.0001
 Commercial / Residential              0.126 0.0266 Inf    1  -9.810  <.0001
 Commercial / Vacant Lot               1.135 0.2940 Inf    1   0.489  0.9966
 Institutional / Park                  1.184 0.2730 Inf    1   0.732  0.9780
 Institutional / (Public Right-Of-Way) 0.408 0.0819 Inf    1  -4.468  0.0001
 Institutional / Residential           0.375 0.0760 Inf    1  -4.841  <.0001
 Institutional / Vacant Lot            3.385 0.8680 Inf    1   4.757  <.0001
 Park / (Public Right-Of-Way)          0.344 0.0743 Inf    1  -4.942  <.0001
 Park / Residential                    0.316 0.0681 Inf    1  -5.348  <.0001
 Park / Vacant Lot                     2.860 0.7700 Inf    1   3.901  0.0013
 (Public Right-Of-Way) / Residential   0.919 0.1700 Inf    1  -0.458  0.9975
 (Public Right-Of-Way) / Vacant Lot    8.301 2.0200 Inf    1   8.696  <.0001
 Residential / Vacant Lot              9.037 2.2100 Inf    1   8.997  <.0001

P value adjustment: tukey method for comparing a family of 6 estimates 
Tests are performed on the log scale 

#EFFECTIVE NUMBER OF SPECIES MODEL OPTIONS ##OPTION 1: lmer

modelgst_ens <- lmer(Cov_q1 ~ GreenSpace + (1|Plot), data= subsite_table)
summary(modelgst_ens) #Random effect explains some variance

Linear mixed model fit by REML ['lmerMod']
Formula: Cov_q1 ~ GreenSpace + (1 | Plot)
   Data: subsite_table

REML criterion at convergence: 676.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.8216 -0.6527 -0.1132  0.5126  2.9882 

Random effects:
 Groups   Name        Variance Std.Dev.
 Plot     (Intercept)  6.28    2.506   
 Residual             55.38    7.442   
Number of obs: 101, groups:  Plot, 22

Fixed effects:
                              Estimate Std. Error t value
(Intercept)                     5.1888     2.0183   2.571
GreenSpaceInstitutional         9.4367     2.6174   3.605
GreenSpacePark                  6.0149     2.8002   2.148
GreenSpacePublic Right-Of-Way  12.7440     2.5111   5.075
GreenSpaceResidential           9.7798     2.5111   3.895
GreenSpaceVacant Lot            0.2622     3.0703   0.085

Correlation of Fixed Effects:
            (Intr) GrnSpI GrnSpP GSPR-O GrnSpR
GrnSpcInstt -0.711                            
GreenSpcPrk -0.669  0.508                     
GrnSpPR-O-W -0.747  0.572  0.538              
GrnSpcRsdnt -0.747  0.572  0.538  0.601       
GrnSpcVcntL -0.600  0.462  0.430  0.483  0.483

# ENS MODEL DIAGNOSTICS
qqnorm(residuals(modelgst_ens))
qqline(residuals(modelgst_ens))

hist(residuals(modelgst_ens))

plot(residuals(modelgst_ens))

shapiro.test(resid(modelgst_ens))

    Shapiro-Wilk normality test

data:  resid(modelgst_ens)
W = 0.95901, p-value = 0.003218

#Diagnostics do not look great, did not continue with this
Anova(modelgst_ens)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: Cov_q1
            Chisq Df Pr(>Chisq)    
GreenSpace 38.421  5  3.105e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise.ens <- emmeans(modelgst_ens, pairwise ~ GreenSpace, type = "response")
summary(Pairwise.ens) 

$emmeans
 GreenSpace          emmean   SE   df lower.CL upper.CL
 Commercial            5.19 2.03 94.2    1.162     9.22
 Institutional        14.63 1.85 93.2   10.949    18.30
 Park                 11.20 2.10 94.6    7.037    15.37
 Public Right-Of-Way  17.93 1.67 91.4   14.607    21.26
 Residential          14.97 1.67 91.4   11.643    18.29
 Vacant Lot            5.45 2.48 95.0    0.522    10.38

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
 contrast                              estimate   SE   df t.ratio p.value
 Commercial - Institutional              -9.437 2.62 78.2  -3.596  0.0072
 Commercial - Park                       -6.015 2.81 82.0  -2.137  0.2790
 Commercial - (Public Right-Of-Way)     -12.744 2.52 79.2  -5.059  <.0001
 Commercial - Residential                -9.780 2.52 79.2  -3.882  0.0028
 Commercial - Vacant Lot                 -0.262 3.09 80.5  -0.085  1.0000
 Institutional - Park                     3.422 2.71 83.5   1.263  0.8039
 Institutional - (Public Right-Of-Way)   -3.307 2.38 77.1  -1.390  0.7327
 Institutional - Residential             -0.343 2.38 77.1  -0.144  1.0000
 Institutional - Vacant Lot               9.174 2.99 82.1   3.065  0.0338
 Park - (Public Right-Of-Way)            -6.729 2.58 79.6  -2.612  0.1062
 Park - Residential                      -3.765 2.58 79.6  -1.462  0.6894
 Park - Vacant Lot                        5.753 3.17 85.3   1.816  0.4610
 (Public Right-Of-Way) - Residential      2.964 2.24 74.5   1.321  0.7725
 (Public Right-Of-Way) - Vacant Lot      12.482 2.90 82.7   4.308  0.0006
 Residential - Vacant Lot                 9.518 2.90 82.7   3.285  0.0181

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 6 estimates 

##OPTION 2: Glmer.nb, where theta is set by package. I do not get the singularity warnings with this model, but do get the warnings for non-count data

modelgst_ens_glmernb <- glmer.nb(Cov_q1 ~ GreenSpace + (1|Plot), data = subsite_table)
#Again, warnings due to non count data
summary(modelgst_ens_glmernb) #random effects are not zero

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Negative Binomial(3.9374)  ( log )
Formula: Cov_q1 ~ GreenSpace + (1 | Plot)
   Data: subsite_table

      AIC       BIC    logLik -2*log(L)  df.resid 
    673.5     694.4    -328.8     657.5        93 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6036 -0.6572 -0.1372  0.6647  2.5985 

Random effects:
 Groups Name        Variance Std.Dev.
 Plot   (Intercept) 0.04965  0.2228  
Number of obs: 101, groups:  Plot, 22

Fixed effects:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    1.56506    0.18446   8.485  < 2e-16 ***
GreenSpaceInstitutional        1.08677    0.22384   4.855 1.20e-06 ***
GreenSpacePark                 0.80937    0.24067   3.363 0.000771 ***
GreenSpacePublic Right-Of-Way  1.29374    0.21467   6.027 1.67e-09 ***
GreenSpaceResidential          1.11093    0.21580   5.148 2.63e-07 ***
GreenSpaceVacant Lot           0.01055    0.28434   0.037 0.970401    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) GrnSpI GrnSpP GSPR-O GrnSpR
GrnSpcInstt -0.761                            
GreenSpcPrk -0.712  0.582                     
GrnSpPR-O-W -0.802  0.659  0.617              
GrnSpcRsdnt -0.795  0.649  0.615  0.686       
GrnSpcVcntL -0.601  0.482  0.446  0.512  0.507

qqnorm(residuals(modelgst_ens_glmernb))
qqline(residuals(modelgst_ens_glmernb))

hist(residuals(modelgst_ens_glmernb))

plot(residuals(modelgst_ens_glmernb))

shapiro.test(resid(modelgst_ens_glmernb))

    Shapiro-Wilk normality test

data:  resid(modelgst_ens_glmernb)
W = 0.99086, p-value = 0.7282

#Residuals look good enough?

Anova(modelgst_ens_glmernb)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: Cov_q1
           Chisq Df Pr(>Chisq)    
GreenSpace 56.78  5  5.614e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise.ensnb <- emmeans(modelgst_ens_glmernb, pairwise ~ GreenSpace, type = "response")
summary(Pairwise.ensnb) 

$emmeans
 GreenSpace          response    SE  df asymp.LCL asymp.UCL
 Commercial              4.78 0.882 Inf      3.33      6.87
 Institutional          14.18 2.070 Inf     10.66     18.87
 Park                   10.74 1.820 Inf      7.71     14.98
 Public Right-Of-Way    17.44 2.250 Inf     13.55     22.45
 Residential            14.53 1.910 Inf     11.22     18.80
 Vacant Lot              4.83 1.100 Inf      3.09      7.55

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                              ratio     SE  df null z.ratio p.value
 Commercial / Institutional            0.337 0.0755 Inf    1  -4.855  <.0001
 Commercial / Park                     0.445 0.1070 Inf    1  -3.363  0.0100
 Commercial / (Public Right-Of-Way)    0.274 0.0589 Inf    1  -6.027  <.0001
 Commercial / Residential              0.329 0.0711 Inf    1  -5.148  <.0001
 Commercial / Vacant Lot               0.990 0.2810 Inf    1  -0.037  1.0000
 Institutional / Park                  1.320 0.2810 Inf    1   1.302  0.7839
 Institutional / (Public Right-Of-Way) 0.813 0.1470 Inf    1  -1.141  0.8641
 Institutional / Residential           0.976 0.1800 Inf    1  -0.131  1.0000
 Institutional / Vacant Lot            2.934 0.7740 Inf    1   4.079  0.0006
 Park / (Public Right-Of-Way)          0.616 0.1240 Inf    1  -2.413  0.1516
 Park / Residential                    0.740 0.1490 Inf    1  -1.496  0.6666
 Park / Vacant Lot                     2.223 0.6200 Inf    1   2.865  0.0479
 (Public Right-Of-Way) / Residential   1.201 0.2050 Inf    1   1.072  0.8925
 (Public Right-Of-Way) / Vacant Lot    3.608 0.9160 Inf    1   5.054  <.0001
 Residential / Vacant Lot              3.005 0.7670 Inf    1   4.311  0.0002

P value adjustment: tukey method for comparing a family of 6 estimates 
Tests are performed on the log scale 

##OPTION 3: try Gamma GLMM with log link

modelgst_ens_gamma <- glmer(Cov_q1 ~ GreenSpace + (1 | Plot), data = subsite_table, family = Gamma(link = "log"))
summary(modelgst_ens_gamma) #Random effect is small

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Gamma  ( log )
Formula: Cov_q1 ~ GreenSpace + (1 | Plot)
   Data: subsite_table

      AIC       BIC    logLik -2*log(L)  df.resid 
    668.4     689.3    -326.2     652.4        93 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.5899 -0.6317 -0.1469  0.6302  2.5609 

Random effects:
 Groups   Name        Variance Std.Dev.
 Plot     (Intercept) 0.07615  0.2759  
 Residual             0.33569  0.5794  
Number of obs: 101, groups:  Plot, 22

Fixed effects:
                              Estimate Std. Error t value Pr(>|z|)    
(Intercept)                    1.50418    0.17866   8.419  < 2e-16 ***
GreenSpaceInstitutional        1.13054    0.21241   5.322 1.02e-07 ***
GreenSpacePark                 0.82831    0.22831   3.628 0.000286 ***
GreenSpacePublic Right-Of-Way  1.34044    0.20232   6.625 3.47e-11 ***
GreenSpaceResidential          1.14824    0.20208   5.682 1.33e-08 ***
GreenSpaceVacant Lot           0.06252    0.25014   0.250 0.802620    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) GrnSpI GrnSpP GSPR-O GrnSpR
GrnSpcInstt -0.648                            
GreenSpcPrk -0.612  0.492                     
GrnSpPR-O-W -0.695  0.578  0.540              
GrnSpcRsdnt -0.688  0.557  0.547  0.613       
GrnSpcVcntL -0.543  0.436  0.405  0.484  0.476

#Residuals
resid_gamma <- residuals(modelgst_ens_gamma, type = "pearson")
qqnorm(resid_gamma); qqline(resid_gamma)

hist(resid_gamma)

plot(resid_gamma)

Anova(modelgst_ens_gamma)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: Cov_q1
            Chisq Df Pr(>Chisq)    
GreenSpace 68.215  5   2.41e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise.ensgamma <- emmeans(modelgst_ens_gamma, pairwise ~ GreenSpace, type = "response")
summary(Pairwise.ensgamma) 

$emmeans
 GreenSpace          response    SE  df asymp.LCL asymp.UCL
 Commercial              4.50 0.804 Inf      3.17      6.39
 Institutional          13.94 2.330 Inf     10.05     19.33
 Park                   10.30 1.900 Inf      7.17     14.80
 Public Right-Of-Way    17.20 2.580 Inf     12.81     23.08
 Residential            14.19 2.150 Inf     10.54     19.11
 Vacant Lot              4.79 1.030 Inf      3.15      7.29

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast                              ratio     SE  df null z.ratio p.value
 Commercial / Institutional            0.323 0.0686 Inf    1  -5.322  <.0001
 Commercial / Park                     0.437 0.0997 Inf    1  -3.628  0.0039
 Commercial / (Public Right-Of-Way)    0.262 0.0530 Inf    1  -6.625  <.0001
 Commercial / Residential              0.317 0.0641 Inf    1  -5.682  <.0001
 Commercial / Vacant Lot               0.939 0.2350 Inf    1  -0.250  0.9999
 Institutional / Park                  1.353 0.3010 Inf    1   1.359  0.7518
 Institutional / (Public Right-Of-Way) 0.811 0.1550 Inf    1  -1.101  0.8812
 Institutional / Residential           0.982 0.1920 Inf    1  -0.091  1.0000
 Institutional / Vacant Lot            2.910 0.7210 Inf    1   4.312  0.0002
 Park / (Public Right-Of-Way)          0.599 0.1240 Inf    1  -2.466  0.1343
 Park / Residential                    0.726 0.1500 Inf    1  -1.551  0.6308
 Park / Vacant Lot                     2.151 0.5630 Inf    1   2.928  0.0400
 (Public Right-Of-Way) / Residential   1.212 0.2160 Inf    1   1.081  0.8893
 (Public Right-Of-Way) / Vacant Lot    3.589 0.8380 Inf    1   5.471  <.0001
 Residential / Vacant Lot              2.962 0.6960 Inf    1   4.619  0.0001

P value adjustment: tukey method for comparing a family of 6 estimates 
Tests are performed on the log scale 

#EVENNESS MODEL OPTIONS - STILL NOT SURE ABOUT WHAT TO DO HERE

##OPTION 1: lmer

modelgst_eve <- lmer(Cov_pie ~ GreenSpace + (1|Plot), data= subsite_table)
summary(modelgst_eve)

Linear mixed model fit by REML ['lmerMod']
Formula: Cov_pie ~ GreenSpace + (1 | Plot)
   Data: subsite_table

REML criterion at convergence: 3.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.5827 -0.2477  0.2054  0.5927  1.4268 

Random effects:
 Groups   Name        Variance Std.Dev.
 Plot     (Intercept) 0.005031 0.07093 
 Residual             0.046892 0.21654 
Number of obs: 101, groups:  Plot, 22

Fixed effects:
                              Estimate Std. Error t value
(Intercept)                    0.62166    0.05859  10.610
GreenSpaceInstitutional        0.16802    0.07614   2.207
GreenSpacePark                 0.18011    0.08144   2.211
GreenSpacePublic Right-Of-Way  0.24213    0.07305   3.315
GreenSpaceResidential          0.12605    0.07305   1.726
GreenSpaceVacant Lot          -0.02886    0.08931  -0.323

Correlation of Fixed Effects:
            (Intr) GrnSpI GrnSpP GSPR-O GrnSpR
GrnSpcInstt -0.713                            
GreenSpcPrk -0.670  0.508                     
GrnSpPR-O-W -0.749  0.572  0.538              
GrnSpcRsdnt -0.749  0.572  0.538  0.601       
GrnSpcVcntL -0.602  0.462  0.431  0.483  0.483

# Evenness MODEL DIAGNOSTICS
qqnorm(residuals(modelgst_eve))
qqline(residuals(modelgst_eve))

hist(residuals(modelgst_eve))

plot(residuals(modelgst_eve))

shapiro.test(resid(modelgst_eve)) #Residuals are not normal, did not continue

    Shapiro-Wilk normality test

data:  resid(modelgst_eve)
W = 0.83919, p-value = 4.246e-09

Anova(modelgst_eve)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: Cov_pie
            Chisq Df Pr(>Chisq)   
GreenSpace 17.508  5   0.003631 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise.eve <- emmeans(modelgst_eve, pairwise ~ GreenSpace, type = "response")
summary(Pairwise.eve) 

$emmeans
 GreenSpace          emmean     SE   df lower.CL upper.CL
 Commercial           0.622 0.0589 94.3    0.505    0.739
 Institutional        0.790 0.0537 93.4    0.683    0.896
 Park                 0.802 0.0609 94.6    0.681    0.923
 Public Right-Of-Way  0.864 0.0486 91.8    0.767    0.960
 Residential          0.748 0.0486 91.8    0.651    0.844
 Vacant Lot           0.593 0.0721 95.0    0.450    0.736

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
 contrast                              estimate     SE   df t.ratio p.value
 Commercial - Institutional             -0.1680 0.0763 78.2  -2.201  0.2492
 Commercial - Park                      -0.1801 0.0819 82.1  -2.200  0.2490
 Commercial - (Public Right-Of-Way)     -0.2421 0.0733 79.3  -3.304  0.0173
 Commercial - Residential               -0.1260 0.0733 79.3  -1.720  0.5229
 Commercial - Vacant Lot                 0.0289 0.0897 80.6   0.322  0.9995
 Institutional - Park                   -0.0121 0.0788 83.6  -0.153  1.0000
 Institutional - (Public Right-Of-Way)  -0.0741 0.0692 77.2  -1.071  0.8914
 Institutional - Residential             0.0420 0.0692 77.2   0.606  0.9903
 Institutional - Vacant Lot              0.1969 0.0871 82.3   2.261  0.2219
 Park - (Public Right-Of-Way)           -0.0620 0.0749 79.7  -0.828  0.9615
 Park - Residential                      0.0541 0.0749 79.7   0.721  0.9788
 Park - Vacant Lot                       0.2090 0.0921 85.5   2.269  0.2184
 (Public Right-Of-Way) - Residential     0.1161 0.0653 74.6   1.778  0.4860
 (Public Right-Of-Way) - Vacant Lot      0.2710 0.0843 82.8   3.216  0.0221
 Residential - Vacant Lot                0.1549 0.0843 82.8   1.838  0.4475

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 6 estimates 

##OPTION 2: glmmTMBm with logit link

#Try beta regression, mixed effects
# Small constant to shrink 0s and 1s inward, because beta regression only handles values between 0 and 1. 
eps <- 0.001
subsite_table$Cov_pie_adj <- subsite_table$Cov_pie
# Shrink 0s up to eps
subsite_table$Cov_pie_adj[subsite_table$Cov_pie_adj == 0] <- eps

modelgst_eve_beta <- glmmTMB(Cov_pie_adj ~ GreenSpace + (1 | Plot), family = beta_family(link = "logit"), data = subsite_table)
summary(modelgst_eve_beta) #Random effects not zero

 Family: beta  ( logit )
Formula:          Cov_pie_adj ~ GreenSpace + (1 | Plot)
Data: subsite_table

      AIC       BIC    logLik -2*log(L)  df.resid 
    -41.4     -20.5      28.7     -57.4        93 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Plot   (Intercept) 0.08983  0.2997  
Number of obs: 101, groups:  Plot, 22

Dispersion parameter for beta family (): 2.24 

Conditional model:
                              Estimate Std. Error z value Pr(>|z|)   
(Intercept)                     0.1924     0.2848   0.675  0.49941   
GreenSpaceInstitutional         0.7297     0.3794   1.923  0.05448 . 
GreenSpacePark                  0.8070     0.4074   1.981  0.04761 * 
GreenSpacePublic Right-Of-Way   0.9835     0.3665   2.684  0.00729 **
GreenSpaceResidential           0.4765     0.3588   1.328  0.18413   
GreenSpaceVacant Lot           -0.3792     0.4339  -0.874  0.38220   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Check model fit
# Simulate residuals from the beta GLMM
sim_res1 <- simulateResiduals(fittedModel = modelgst_eve_beta)
# Plot residual diagnostics
plot(sim_res1)

# Perform formal tests for uniformity, dispersion, zero-inflation etc.
testUniformity(sim_res1)

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.25719, p-value = 3.149e-06
alternative hypothesis: two-sided

testDispersion(sim_res1)

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 0.72351, p-value = 0.008
alternative hypothesis: two.sided

#NOT A GOOD FIT?

Anova(modelgst_eve_beta)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: Cov_pie_adj
            Chisq Df Pr(>Chisq)   
GreenSpace 15.745  5    0.00761 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise.evebeta <- emmeans(modelgst_eve_beta, pairwise ~ GreenSpace, type = "response")
summary(Pairwise.evebeta) 

$emmeans
 GreenSpace          response     SE  df asymp.LCL asymp.UCL
 Commercial             0.548 0.0705 Inf     0.410     0.679
 Institutional          0.715 0.0552 Inf     0.596     0.811
 Park                   0.731 0.0612 Inf     0.596     0.833
 Public Right-Of-Way    0.764 0.0460 Inf     0.663     0.842
 Residential            0.661 0.0537 Inf     0.550     0.757
 Vacant Lot             0.453 0.0862 Inf     0.296     0.621

Confidence level used: 0.95 
Intervals are back-transformed from the logit scale 

$contrasts
 contrast                              odds.ratio    SE  df null z.ratio
 Commercial / Institutional                 0.482 0.183 Inf    1  -1.923
 Commercial / Park                          0.446 0.182 Inf    1  -1.981
 Commercial / (Public Right-Of-Way)         0.374 0.137 Inf    1  -2.683
 Commercial / Residential                   0.621 0.223 Inf    1  -1.328
 Commercial / Vacant Lot                    1.461 0.634 Inf    1   0.874
 Institutional / Park                       0.926 0.366 Inf    1  -0.196
 Institutional / (Public Right-Of-Way)      0.776 0.271 Inf    1  -0.726
 Institutional / Residential                1.288 0.444 Inf    1   0.734
 Institutional / Vacant Lot                 3.031 1.300 Inf    1   2.585
 Park / (Public Right-Of-Way)               0.838 0.314 Inf    1  -0.471
 Park / Residential                         1.392 0.520 Inf    1   0.885
 Park / Vacant Lot                          3.274 1.460 Inf    1   2.652
 (Public Right-Of-Way) / Residential        1.660 0.545 Inf    1   1.546
 (Public Right-Of-Way) / Vacant Lot         3.907 1.610 Inf    1   3.309
 Residential / Vacant Lot                   2.353 0.961 Inf    1   2.095
 p.value
  0.3879
  0.3533
  0.0785
  0.7694
  0.9528
  1.0000
  0.9788
  0.9777
  0.1008
  0.9971
  0.9502
  0.0852
  0.6343
  0.0120
  0.2895

P value adjustment: tukey method for comparing a family of 6 estimates 
Tests are performed on the log odds ratio scale 

##OPTION 3: zero inflated beta regression, but this is complicated and I am not sure if it is the right fit

modelgst_eve_beta_01 <- glmmTMB(Cov_pie ~ GreenSpace + (1 | Plot), 
                                family = beta_family(link = "logit"),
                                ziformula = ~1,  
                                data = subsite_table
)

summary(modelgst_eve_beta_01)

 Family: beta  ( logit )
Formula:          Cov_pie ~ GreenSpace + (1 | Plot)
Zero inflation:           ~1
Data: subsite_table

      AIC       BIC    logLik -2*log(L)  df.resid 
    -75.7     -52.2      46.9     -93.7        92 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Plot   (Intercept) 0.1244   0.3527  
Number of obs: 101, groups:  Plot, 22

Dispersion parameter for beta family (): 8.74 

Conditional model:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     0.6746     0.1958   3.445 0.000571 ***
GreenSpaceInstitutional         0.9978     0.2670   3.737 0.000186 ***
GreenSpacePark                  0.7120     0.2728   2.610 0.009051 ** 
GreenSpacePublic Right-Of-Way   1.0286     0.2517   4.087 4.37e-05 ***
GreenSpaceResidential           0.6084     0.2454   2.479 0.013168 *  
GreenSpaceVacant Lot            0.4114     0.3155   1.304 0.192233    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Zero-inflation model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -2.9549     0.4587  -6.442 1.18e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Diagnostics with DHARMa
sim_res <- simulateResiduals(modelgst_eve_beta_01)
plot(sim_res)

testUniformity(sim_res)

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.13406, p-value = 0.05301
alternative hypothesis: two-sided

testDispersion(sim_res)

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 1.0927, p-value = 0.672
alternative hypothesis: two.sided

Anova(modelgst_eve_beta_01)

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: Cov_pie
            Chisq Df Pr(>Chisq)    
GreenSpace 21.595  5  0.0006251 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise.evezero <- emmeans(modelgst_eve_beta_01, pairwise ~ GreenSpace, type = "response")
summary(Pairwise.evezero) 

$emmeans
 GreenSpace          response     SE  df asymp.LCL asymp.UCL
 Commercial             0.663 0.0438 Inf     0.572     0.742
 Institutional          0.842 0.0285 Inf     0.778     0.890
 Park                   0.800 0.0352 Inf     0.722     0.860
 Public Right-Of-Way    0.846 0.0253 Inf     0.790     0.889
 Residential            0.783 0.0312 Inf     0.716     0.838
 Vacant Lot             0.748 0.0517 Inf     0.634     0.835

Confidence level used: 0.95 
Intervals are back-transformed from the logit scale 

$contrasts
 contrast                              odds.ratio     SE  df null z.ratio
 Commercial / Institutional                 0.369 0.0984 Inf    1  -3.737
 Commercial / Park                          0.491 0.1340 Inf    1  -2.610
 Commercial / (Public Right-Of-Way)         0.358 0.0900 Inf    1  -4.087
 Commercial / Residential                   0.544 0.1340 Inf    1  -2.479
 Commercial / Vacant Lot                    0.663 0.2090 Inf    1  -1.304
 Institutional / Park                       1.331 0.3810 Inf    1   0.998
 Institutional / (Public Right-Of-Way)      0.970 0.2540 Inf    1  -0.118
 Institutional / Residential                1.476 0.3820 Inf    1   1.505
 Institutional / Vacant Lot                 1.798 0.5930 Inf    1   1.778
 Park / (Public Right-Of-Way)               0.729 0.1950 Inf    1  -1.186
 Park / Residential                         1.109 0.2880 Inf    1   0.398
 Park / Vacant Lot                          1.351 0.4490 Inf    1   0.903
 (Public Right-Of-Way) / Residential        1.522 0.3610 Inf    1   1.771
 (Public Right-Of-Way) / Vacant Lot         1.854 0.5730 Inf    1   1.996
 Residential / Vacant Lot                   1.218 0.3720 Inf    1   0.645
 p.value
  0.0026
  0.0947
  0.0006
  0.1302
  0.7830
  0.9189
  1.0000
  0.6611
  0.4801
  0.8439
  0.9987
  0.9458
  0.4846
  0.3444
  0.9876

P value adjustment: tukey method for comparing a family of 6 estimates 
Tests are performed on the log odds ratio scale