---
# System prepended metadata

title: CCFRP Data Analysis + Statistics and R Notes

---

# CCFRP Data Analysis + Statistics and R Notes
Modeling in Ecology:
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5970551/
## Stats Notes
**ANOVA**- one-way analysis of variance used to determine whether there are any statistically significant differences between the means of three or more independent (unrelated) groups

*can tell you that two groups are statistically different from each other but not which two groups for that you need a post hoc test. 

y = B0 + B1[x] + e 
Categorical value for X 

**T-test**- compare two groups only 
`t.test(DependentVariable~category, data=dataset)`
When Category has only two options 
DependentVariable is continuous

**One-way analysis of means** (not assuming equal variances) 
`oneway.test(DependentVariable~category, data=dataset)`
When Category has more than two options 

*or* 
**lm test** (assuming equal variances)
`summary(lm(DependentVariable~category,data=dataset)`
Gives: 
```
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                     4.1737     1.0461   3.990 7.27e-05 ***
Vertical_DistWC                 4.1273     1.4007   2.947  0.00331 ** 
siteREF                        -1.3938     1.4470  -0.963  0.33575    
locSP                          -1.3506     1.3976  -0.966  0.33418    
siteREF:locSP                  -0.6799     2.0008  -0.340  0.73409    
Vertical_DistWC:siteMPA:locBH  -4.3670     2.3781  -1.836  0.06670 .  
Vertical_DistWC:siteREF:locBH  -1.3121     2.2106  -0.594  0.55298    
Vertical_DistWC:siteMPA:locSP   8.7787     1.9130   4.589 5.23e-06 ***
Vertical_DistWC:siteREF:locSP       NA         NA      NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.2 on 743 degrees of freedom
Multiple R-squared:  0.1817,	Adjusted R-squared:  0.174 
F-statistic: 23.57 on 7 and 743 DF,  p-value: < 2.2e-16
```



Analysis of variance table: displays the significance of categories
`anova(lm(DependentVariable~category,data=dataset)`

```
Analysis of Variance Table

Response: CPUE_by_Guild
                        Df Sum Sq Mean Sq F value    Pr(>F)    
Vertical_Dist            1   7363  7362.7 70.8250 < 2.2e-16 ***
site                     1   3523  3522.8 33.8872 8.687e-09 ***
loc                      1    144   144.2  1.3874 0.2392192    
site:loc                 1   1564  1564.3 15.0475 0.0001141 ***
Vertical_Dist:site:loc   3   4556  1518.7 14.6093 2.974e-09 ***
Residuals              743  77239   104.0                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```
Use a tukey test to determine which interactions are significant 
```
m1 <- lm(CPUE_by_Guild ~ Vertical_Dist + site + loc + site:loc + site:loc:Vertical_Dist, data = df.sum.across.guilds.per.drift)
m1_tukey <- TukeyHSD(aov(m1))
m1_tukey
```

```
$`Vertical_Dist:site:loc`
                            diff         lwr         upr
WC:MPA:BH-BT:MPA:BH  -0.23966429  -6.0810098   5.6016813
BT:REF:BH-BT:MPA:BH  -1.39380690  -5.7921173   3.0045035
WC:REF:BH-BT:MPA:BH   1.42141188  -3.8602703   6.7030941
BT:MPA:SP-BT:MPA:BH  -1.35062139  -5.5988727   2.8976299
WC:MPA:SP-BT:MPA:BH  11.55540186   7.3298132  15.7809906
BT:REF:SP-BT:MPA:BH  -3.42431515  -7.8753324   1.0267021
WC:REF:SP-BT:MPA:BH   0.70302295  -3.6022107   5.0082566
BT:REF:BH-WC:MPA:BH  -1.15414261  -6.9201128   4.6118276
WC:REF:BH-WC:MPA:BH   1.66107617  -4.8040181   8.1261705
BT:MPA:SP-WC:MPA:BH  -1.11095710  -6.7632944   4.5413802
WC:MPA:SP-WC:MPA:BH  11.79506615   6.1597421  17.4303902
BT:REF:SP-WC:MPA:BH  -3.18465086  -8.9909261   2.6216244
WC:REF:SP-WC:MPA:BH   0.94268724  -4.7526016   6.6379760
WC:REF:BH-BT:REF:BH   2.81521878  -2.3829791   8.0134166
BT:MPA:SP-BT:REF:BH   0.04318551  -4.1008142   4.1871852
WC:MPA:SP-BT:REF:BH  12.94920876   8.8284449  17.0699726
BT:REF:SP-BT:REF:BH  -2.03050825  -6.3821343   2.3211178
WC:REF:SP-BT:REF:BH   2.09682985  -2.1055660   6.2992257
BT:MPA:SP-WC:REF:BH  -2.77203327  -7.8438934   2.2998269
WC:MPA:SP-WC:REF:BH  10.13398998   5.0810971  15.1868828
BT:REF:SP-WC:REF:BH  -4.84572703 -10.0885966   0.3971425
WC:REF:SP-WC:REF:BH  -0.71838894  -5.8380728   4.4012950
WC:MPA:SP-BT:MPA:SP  12.90602325   8.9458214  16.8662251
BT:REF:SP-BT:MPA:SP  -2.07369376  -6.2735930   2.1262055
WC:REF:SP-BT:MPA:SP   2.05364434  -1.9914312   6.0987198
BT:REF:SP-WC:MPA:SP -14.97971701 -19.1566913 -10.8027427
WC:REF:SP-WC:MPA:SP -10.85237892 -14.8736470  -6.8311109
WC:REF:SP-BT:REF:SP   4.12733809  -0.1301907   8.3848669
                        p adj
WC:MPA:BH-BT:MPA:BH 1.0000000
BT:REF:BH-BT:MPA:BH 0.9793221
WC:REF:BH-BT:MPA:BH 0.9921127
BT:MPA:SP-BT:MPA:BH 0.9789357
WC:MPA:SP-BT:MPA:BH 0.0000000
BT:REF:SP-BT:MPA:BH 0.2741201
WC:REF:SP-BT:MPA:BH 0.9996771
BT:REF:BH-WC:MPA:BH 0.9987720
WC:REF:BH-WC:MPA:BH 0.9940555
BT:MPA:SP-WC:MPA:BH 0.9989088
WC:MPA:SP-WC:MPA:BH 0.0000000
BT:REF:SP-WC:MPA:BH 0.7085774
WC:REF:SP-WC:MPA:BH 0.9996468
WC:REF:BH-BT:REF:BH 0.7218974
BT:MPA:SP-BT:REF:BH 1.0000000
WC:MPA:SP-BT:REF:BH 0.0000000
BT:REF:SP-BT:REF:BH 0.8488074
WC:REF:SP-BT:REF:BH 0.7984875
BT:MPA:SP-WC:REF:BH 0.7123254
WC:MPA:SP-WC:REF:BH 0.0000000
BT:REF:SP-WC:REF:BH 0.0941911
WC:REF:SP-WC:REF:BH 0.9998830
WC:MPA:SP-BT:MPA:SP 0.0000000
BT:REF:SP-BT:MPA:SP 0.8070926
WC:REF:SP-BT:MPA:SP 0.7836747
BT:REF:SP-WC:MPA:SP 0.0000000
WC:REF:SP-WC:MPA:SP 0.0000000
WC:REF:SP-BT:REF:SP 0.0651076
```

Only want 
WC:REF:SP-WC:MPA:SP 0.0000000
BT:REF:SP-BT:MPA:SP 0.8070926
WC:REF:BH-WC:MPA:BH 0.9940555
BT:REF:BH-BT:MPA:BH 0.9793221

Does not meet regression assumptions and is over penalising 

Interaction plot: MPA/REF vs BH/SP
```
with(df.sum.across.guilds.per.drift,interaction.plot(x.factor = loc,site,response = CPUE_by_Guild))
```
![reference link](https://i.imgur.com/X8DTSA9.png)

Interaction plot 
loc,Vertical_Dist
![](https://i.imgur.com/DOizQLQ.png)
Interaction plot 
site,Vertical_Dist
![](https://i.imgur.com/frIOAhC.png)
###   Vertical Distribution Graphs 
![](https://i.imgur.com/b4gAqeA.png)
p-values calculated in the section below:
WC_MPA_BH - WC_REF_BH 0.1733 

BT_MPA_BH - BT_REF_BH 0.0304

![](https://i.imgur.com/42zscRl.png)
p-values calculated in the section below:
WC_MPA_SP - WC_REF_SP   <.0001 

BT_MPA_SP - BT_REF_SP <.0001 
## **ANOVA and TukeyHSD post-hoc with specific interactions**
```
m1 <- glm(CPUE_by_Guild ~ Vertical_Dist + site + loc + site:loc + site:loc:Vertical_Dist, data = df.sum.across.guilds.per.drift)
anova(m1, test = "F")
```
Analysis of Deviance Table
```

                       Df Deviance Resid. Df Resid. Dev       F    Pr(>F)    
NULL                                     750      94390                      
Vertical_Dist           1   7362.7       749      87027 70.8250 < 2.2e-16 ***
site                    1   3522.8       748      83504 33.8872 8.687e-09 ***
loc                     1    144.2       747      83360  1.3874 0.2392192    
site:loc                1   1564.3       746      81796 15.0475 0.0001141 ***
Vertical_Dist:site:loc  3   4556.2       743      77239 14.6093 2.974e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```
Test if the model fits: 
`hist(residuals(m1))`
![](https://i.imgur.com/LWarhpk.png)
Does not fit well (assues normal distribution): looks like a gamma distribution which can be often log into a normal distribution with 
`glm(…, family = Gamma(link = log))`

so
```
m1 <- glm((CPUE_by_Guild + 1) ~ Vertical_Dist + site + loc + site:loc + site:loc:Vertical_Dist, family = Gamma(link = log), data =
df.sum.across.guilds.per.drift)
```
(CPUE_by_Guild + 1) becuase you can't have log(0)

To fit the model better 
![](https://i.imgur.com/RDgBMea.png)


Need only specific contrasts  for the TukeyHSD post-hoc 

Tried: 
``` 
emmeans::ref_grid(m1, at = list(Vertical_Dist = 0:1, site = 0:1, loc = 0:1))
```

error: Error in ``contrasts<-`(`*tmp*`, value = contrasts.arg[[nn]]) : 
  contrasts apply only to factors``
  
  To fix error: I changed loc, site, and Vertical_Dist into factors 
 ```
df.sum.across.guilds.per.drift$loc <- factor(df.sum.across.guilds.per.drift$loc)
df.sum.across.guilds.per.drift$site <- factor(df.sum.across.guilds.per.drift$site)
df.sum.across.guilds.per.drift$Vertical_Dist <- factor(df.sum.across.guilds.per.drift$Vertical_Dist)
```
Then ran the model and created the ref_grid 
```
#create model 
m1 <- glm((CPUE_by_Guild + 1) ~ Vertical_Dist + site + loc + site:loc + site:loc:Vertical_Dist, family = Gamma(link = log), data = df.sum.across.guilds.per.drift)
#specify contrasts 
m2<- emmeans::ref_grid(m1, at = c(Vertical_Dist = 0:1, site = 0:1, loc = 0:1))
```
  
  great referance for Anova and fitting models: https://www.r-bloggers.com/linear-models-anova-glms-and-mixed-effects-models-in-r/

Now `plot(m2)` gives: 
![](https://i.imgur.com/KULO1P6.png)

Now I am unsure how to use the ref_grid to get only the contrasts I want, I can get all of them again pairwise: 

 `pairs(m2)` gives:
``` 
contrast              estimate    SE  df z.ratio p.value
 BT,MPA,BH - WC,MPA,BH   0.0474 0.193 Inf   0.246 1.0000 
 **BT,MPA,BH - BT,REF,BH**   0.3139 0.145 Inf   2.165 0.3730 
 BT,MPA,BH - WC,REF,BH  -0.2427 0.174 Inf  -1.394 0.8602 
 BT,MPA,BH - BT,MPA,SP   0.3025 0.140 Inf   2.161 0.3759 
 BT,MPA,BH - WC,MPA,SP  -1.1736 0.139 Inf  -8.426 <.0001 
 BT,MPA,BH - BT,REF,SP   1.0843 0.147 Inf   7.391 <.0001 
 BT,MPA,BH - WC,REF,SP  -0.1274 0.142 Inf  -0.898 0.9863 
 WC,MPA,BH - BT,REF,BH   0.2665 0.190 Inf   1.402 0.8567 
 **WC,MPA,BH - WC,REF,BH**  -0.2902 0.213 Inf  -1.362 0.8746 
 WC,MPA,BH - BT,MPA,SP   0.2551 0.186 Inf   1.369 0.8713 
 WC,MPA,BH - WC,MPA,SP  -1.2210 0.186 Inf  -6.574 <.0001 
 WC,MPA,BH - BT,REF,SP   1.0369 0.191 Inf   5.418 <.0001 
 WC,MPA,BH - WC,REF,SP  -0.1748 0.188 Inf  -0.931 0.9831 
 BT,REF,BH - WC,REF,BH  -0.5566 0.171 Inf  -3.249 0.0256 
 BT,REF,BH - BT,MPA,SP  -0.0114 0.137 Inf  -0.083 1.0000 
 BT,REF,BH - WC,MPA,SP  -1.4875 0.136 Inf -10.952 <.0001 
 BT,REF,BH - BT,REF,SP   0.7704 0.143 Inf   5.372 <.0001 
 BT,REF,BH - WC,REF,SP  -0.4413 0.139 Inf  -3.186 0.0312 
 WC,REF,BH - BT,MPA,SP   0.5453 0.167 Inf   3.262 0.0245 
 WC,REF,BH - WC,MPA,SP  -0.9308 0.167 Inf  -5.589 <.0001 
 WC,REF,BH - BT,REF,SP   1.3271 0.173 Inf   7.680 <.0001 
 WC,REF,BH - WC,REF,SP   0.1153 0.169 Inf   0.683 0.9974 
 BT,MPA,SP - WC,MPA,SP  -1.4761 0.131 Inf -11.309 <.0001 
 **BT,MPA,SP - BT,REF,SP**   0.7818 0.138 Inf   5.648 <.0001 
 BT,MPA,SP - WC,REF,SP  -0.4299 0.133 Inf  -3.225 0.0276 
 WC,MPA,SP - BT,REF,SP   2.2579 0.138 Inf  16.401 <.0001 
 **WC,MPA,SP - WC,REF,SP**   1.0462 0.133 Inf   7.893 <.0001 
 BT,REF,SP - WC,REF,SP  -1.2117 0.140 Inf  -8.635 <.0001 
 ```
Only want             p-value
WC,MPA,SP - WC,REF,SP  <.0001 
BT,MPA,SP - BT,REF,SP  <.0001
WC,MPA,BH - WC,REF,BH  0.8746
BT,MPA,BH - BT,REF,BH  0.3730
*These p-values seem more reasonable are they still over penalizing because we are no longer using Tukey??? 
**Values confirmed by Connor **

Pairwise P-Value Plot:
`emmeans::pwpp(m2)`
![](https://i.imgur.com/lIoJbvV.png)
```
> multcomp::cld(m2)
 Vertical_Dist site loc prediction     SE  df .group
 BT            REF  SP       0.559 0.1027 Inf  1    
 BT            REF  BH       1.330 0.1002 Inf   2   
 BT            MPA  SP       1.341 0.0929 Inf   2   
 WC            MPA  BH       1.596 0.1615 Inf   23  
 BT            MPA  BH       1.644 0.1048 Inf   23  
 WC            REF  SP       1.771 0.0957 Inf    3  
 WC            REF  BH       1.886 0.1390 Inf    3  
 WC            MPA  SP       2.817 0.0917 Inf     4 
```

Bonforoni Test 
Papers to refer to: 
Mace Morgan 2006, 2011 
invertabrate recruitment in low at the headland and higher in the lee of the headland 

Review paper 2014 Morgan 
advances in oceanography 

### Code for Vertical_Dist Direct Comparsions
```
refGrid <- emmeans::ref_grid(m1, at = c(Vertical_Dist = 0:1, site = 0:1, loc = 0:1))
refGrid.s <- emmeans::emmeans(refGrid, specs = c("Vertical_Dist", "site", "loc"))
refGrid.s <- emmeans::emmeans(m1, specs = c("Vertical_Dist", "site", "loc"))
refGrid.s

WC_MPA_SP = c(0, 0, 0, 0, 0, 1, 0, 0)
WC_REF_SP = c(0, 0, 0, 0, 0, 0, 0, 1)

BT_MPA_SP = c(0, 0, 0, 0, 1, 0, 0, 0)
BT_REF_SP = c(0, 0, 0, 0, 0, 0, 1, 0)

WC_MPA_BH = c(0, 1, 0, 0, 0, 0, 0, 0)
WC_REF_BH = c(0, 0, 0, 1, 0, 0, 0, 0)

BT_MPA_BH = c(1, 0, 0, 0, 0, 0, 0, 0)
BT_REF_BH = c(0, 0, 1, 0, 0, 0, 0, 0)

emmeans::contrast(refGrid.s, method = list(WC_MPA_SP - WC_REF_SP) )
emmeans::contrast(refGrid.s, method = list(BT_MPA_SP - BT_REF_SP) )
emmeans::contrast(refGrid.s, method = list(WC_MPA_BH - WC_REF_BH) )
emmeans::contrast(refGrid.s, method = list(BT_MPA_BH - BT_REF_BH) )
```
### Output 
```
> refGrid.s
 Vertical_Dist site loc emmean     SE  df asymp.LCL asymp.UCL
 BT            MPA  BH   1.644 0.1048 Inf     1.438      1.85
 WC            MPA  BH   1.596 0.1615 Inf     1.280      1.91
 BT            REF  BH   1.330 0.1002 Inf     1.133      1.53
 WC            REF  BH   1.886 0.1390 Inf     1.614      2.16
 BT            MPA  SP   1.341 0.0929 Inf     1.159      1.52
 WC            MPA  SP   2.817 0.0917 Inf     2.637      3.00
 BT            REF  SP   0.559 0.1027 Inf     0.358      0.76
 WC            REF  SP   1.771 0.0957 Inf     1.583      1.96

Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 
> emmeans::contrast(refGrid.s, method = list(WC_MPA_SP - WC_REF_SP) )
 contrast                   estimate    SE  df z.ratio p.value
 c(0, 0, 0, 0, 0, 1, 0, -1)     1.05 0.133 Inf 7.893   <.0001 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(BT_MPA_SP - BT_REF_SP) )
 contrast                   estimate    SE  df z.ratio p.value
 c(0, 0, 0, 0, 1, 0, -1, 0)    0.782 0.138 Inf 5.648   <.0001 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(WC_MPA_BH - WC_REF_BH) )
 contrast                   estimate    SE  df z.ratio p.value
 c(0, 1, 0, -1, 0, 0, 0, 0)    -0.29 0.213 Inf -1.362  0.1733 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(BT_MPA_BH - BT_REF_BH) )
 contrast                   estimate    SE  df z.ratio p.value
 c(1, 0, -1, 0, 0, 0, 0, 0)    0.314 0.145 Inf 2.165   0.0304 
```
### North Central MPAs
![](https://i.imgur.com/Ode53MI.png)
MPA,BH - REF,BH 0.8586
MPA,SP - REF,SP <.0001
### Model Simplified for only Site and Loc 
![](https://i.imgur.com/H405t8X.png)
Pairwise
```
> pairs(m2)
 contrast        estimate    SE  df z.ratio p.value
 MPA,BH - REF,BH   0.0296 0.166 Inf  0.178  0.9980 
 MPA,BH - MPA,SP  -1.0482 0.160 Inf -6.549  <.0001 
 MPA,BH - REF,SP   0.0527 0.163 Inf  0.324  0.9883 
 REF,BH - MPA,SP  -1.0778 0.156 Inf -6.923  <.0001 
 REF,BH - REF,SP   0.0231 0.158 Inf  0.146  0.9989 
 MPA,SP - REF,SP   1.1010 0.152 Inf  7.245  <.0001 
```
MPA,BH - REF,BH 0.9980
MPA,SP - REF,SP <.0001 

#### Direct Contrasts 
```
 site loc emmean    SE  df z.ratio p.value
 MPA  BH    1.92 0.120 Inf 15.953  <.0001 
 REF  BH    1.89 0.115 Inf 16.513  <.0001 
 MPA  SP    2.97 0.105 Inf 28.170  <.0001 
 REF  SP    1.87 0.109 Inf 17.071  <.0001 

Results are given on the log(mu + 1) (not the response) scale. 
P value adjustment: bonferroni method for 4 tests 
> #pairs(refGrid.s)
> 
> MPA_BH = c(1, 0, 0, 0)
> REF_BH = c(0, 1, 0, 0)
> 
> MPA_SP = c(0, 0, 1, 0)
> REF_SP = c(0, 0, 0, 1)
> 
> 
> emmeans::contrast(refGrid.s, method = list(MPA_BH - REF_BH) )
 contrast       estimate    SE  df z.ratio p.value
 c(1, -1, 0, 0)   0.0296 0.166 Inf 0.178   0.8586 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list( MPA_SP - REF_SP) )
 contrast       estimate    SE  df z.ratio p.value
 c(0, 0, 1, -1)      1.1 0.152 Inf 7.245   <.0001 

Results are given on the log (not the response) scale.
```

## Species MPA vs. REF 

![](https://i.imgur.com/FenYIUC.png)
![](https://i.imgur.com/FQYcvy2.png)

**Key: * = significant** 

Deacon* (SP & BH)
```
emmeans::contrast(refGrid.s, method = list(DEA_MPA_SP - DEA_REF_SP) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.767
    SE  df z.ratio p.value
 0.125 Inf 6.133   <.0001 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(DEA_MPA_BH - DEA_REF_BH) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,     -0.54
    SE  df z.ratio p.value
 0.235 Inf -2.297  0.0216 

Results are given on the log (not the response) scale. 
> 
```
Blue* (SP)
```
> emmeans::contrast(refGrid.s, method = list(BLU_MPA_SP - BLU_REF_SP) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.919
    SE  df z.ratio p.value
 0.126 Inf 7.321   <.0001 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(BLU_MPA_BH - BLU_REF_BH) )
 contrast                                                       estimate
 c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.446
    SE  df z.ratio p.value
 0.271 Inf 1.646   0.0997 
```
Vermillion
```
> emmeans::contrast(refGrid.s, method = list(VER_MPA_SP - VER_REF_SP) )
 contrast                                                       estimate    SE  df
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.262 0.257 Inf
 z.ratio p.value
 1.019   0.3082 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(VER_MPA_BH - VER_REF_BH) )
 contrast                                                       estimate    SE  df
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    -0.387 0.385 Inf
 z.ratio p.value
 -1.005  0.3149
 ```
 Lingcod 
 ```
  emmeans::contrast(refGrid.s, method = list(LING_MPA_SP - LING_REF_SP) )
 contrast                                                       estimate    SE  df z.ratio
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.185 0.186 Inf 0.995  
 p.value
 0.3197 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(LING_MPA_BH - LING_REF_BH) )
 contrast                                                       estimate    SE  df z.ratio
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,     0.207 0.163 Inf 1.273  
 p.value
 0.2031 
 ```
 Black* (SP)
 ```
 > emmeans::contrast(refGrid.s, method = list(BLA_MPA_SP - BLA_REF_SP) )
  contrast                                                       estimate    SE  df z.ratio
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.412 0.154 Inf 2.675  
 p.value
 0.0075 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(BLA_MPA_BH - BLA_REF_BH) )
 contrast                                                       estimate   SE  df z.ratio
 c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    -0.326 0.27 Inf -1.208 
 p.value
 0.2272 
 ```
 Brown 
 ```
 > emmeans::contrast(refGrid.s, method = list(BWN_MPA_SP - BWN_REF_SP) )
 contrast                                                       estimate    SE  df z.ratio
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.256 0.591 Inf 0.433  
 p.value
 0.6649 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(BWN_MPA_BH - BWN_REF_BH) )
 contrast                                                       estimate    SE  df z.ratio
 c(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.187 0.177 Inf 1.056  
 p.value
 0.2909 
 ```
 China 
 ```
 > emmeans::contrast(refGrid.s, method = list(CHN_MPA_SP - CHN_REF_SP) )
 contrast                                                       estimate    SE  df z.ratio
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.233 0.166 Inf 1.399  
 p.value
 0.1618 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(CHN_MPA_BH - CHN_REF_BH) )
 contrast                                                       estimate    SE  df z.ratio
 c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.235 0.242 Inf 0.973  
 p.value
 0.3307 
 ```
 Gopher* (SP)
 ```
 > emmeans::contrast(refGrid.s, method = list(GPR_MPA_SP - GPR_REF_SP) )
 contrast                                                       estimate   SE  df z.ratio
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.453 0.14 Inf 3.246  
 p.value
 0.0012 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(GPR_MPA_BH - GPR_REF_BH) )
 contrast                                                       estimate   SE  df z.ratio
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,    -0.108 0.17 Inf -0.637 
 p.value
 0.5241 
 ```
 Olive* (SP)
 ```
 emmeans::contrast(refGrid.s, method = list(OLV_MPA_SP - OLV_REF_SP) )
 contrast                                                       estimate    SE  df z.ratio
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.463 0.156 Inf 2.967  
 p.value
 0.0030 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(OLV_MPA_BH - OLV_REF_BH) )
 contrast                                                       estimate    SE  df z.ratio
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,    -0.155 0.889 Inf -0.174 
 p.value
 0.8619 
```
Quillback 
```
> emmeans::contrast(refGrid.s, method = list(QBK_MPA_SP - QBK_REF_SP) )
 contrast                                                       estimate    SE  df z.ratio
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.419 0.965 Inf 0.434  
 p.value
 0.6640 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(QBK_MPA_BH - QBK_REF_BH) )
 contrast                                                       estimate    SE  df z.ratio
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,     0.126 0.689 Inf 0.182  
 p.value
 0.8552 
 ```
 
 ## Species Length MPA vs. REF
![](https://i.imgur.com/Tgu5S0J.png)
![](https://i.imgur.com/XLnUoQM.png)


Deacon* (SP) 
```
emmeans::contrast(refGrid.s, method = list(DEA_MPA_SP - DEA_REF_SP) )
 contrast                                                      
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
 estimate      SE  df z.ratio p.value
   0.0445 0.00846 Inf 5.260   <.0001 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(DEA_MPA_BH - DEA_REF_BH) )
 contrast                                                      
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
 estimate    SE  df z.ratio p.value
   0.0178 0.025 Inf 0.713   0.4758 
   ```
   Blue* (SP & BH)
   ```
   > emmeans::contrast(refGrid.s, method = list(BLU_MPA_SP - BLU_REF_SP) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.035
      SE  df z.ratio p.value
 0.00861 Inf 4.063   <.0001 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(BLU_MPA_BH - BLU_REF_BH) )
 contrast                                                       estimate
 c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.127
     SE  df z.ratio p.value
 0.0277 Inf 4.599   <.0001 
 ```
 Vermillion
 ```
 > emmeans::contrast(refGrid.s, method = list(VER_MPA_SP - VER_REF_SP) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    0.0328
     SE  df z.ratio p.value
 0.0534 Inf 0.615   0.5387 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(VER_MPA_BH - VER_REF_BH) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    0.0767
     SE  df z.ratio p.value
 0.0735 Inf 1.044   0.2966 
 ```
 Lingcod* (SP & BH)
 ```
 > emmeans::contrast(refGrid.s, method = list(LING_MPA_SP - LING_REF_SP) )
 contrast                                                       estimate    SE
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    0.0929 0.036
  df z.ratio p.value
 Inf 2.582   0.0098 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(LING_MPA_BH - LING_REF_BH) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,    0.0623
     SE  df z.ratio p.value
 0.0281 Inf 2.221   0.0264 
```
Black*
```
> emmeans::contrast(refGrid.s, method = list(BLA_MPA_SP - BLA_REF_SP) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    0.0436
     SE  df z.ratio p.value
 0.0161 Inf 2.701   0.0069 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(BLA_MPA_BH - BLA_REF_BH) )
 contrast                                                       estimate
 c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.102
     SE  df z.ratio p.value
 0.0361 Inf 2.819   0.0048 
 ```
 Brown
 ```
 > emmeans::contrast(refGrid.s, method = list(BWN_MPA_SP - BWN_REF_SP) )
 contrast                                                       estimate    SE
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    0.0795 0.124
  df z.ratio p.value
 Inf 0.643   0.5202 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(BWN_MPA_BH - BWN_REF_BH) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    0.0244
     SE  df z.ratio p.value
 0.0227 Inf 1.078   0.2810 
 ```
 China* (SP)
 ```
 > emmeans::contrast(refGrid.s, method = list(CHN_MPA_SP - CHN_REF_SP) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    0.0866
     SE  df z.ratio p.value
 0.0287 Inf 3.015   0.0026 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(CHN_MPA_BH - CHN_REF_BH) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    0.0399
     SE  df z.ratio p.value
 0.0379 Inf 1.053   0.2926 
 ```
 Gopher* (SP & BH)
 ```
 > emmeans::contrast(refGrid.s, method = list(GPR_MPA_SP - GPR_REF_SP) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.105
     SE  df z.ratio p.value
 0.0196 Inf 5.351   <.0001 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(GPR_MPA_BH - GPR_REF_BH) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,    0.0807
     SE  df z.ratio p.value
 0.0259 Inf 3.111   0.0019 
 ```
 Olive* (SP & BH)
 ```
 > emmeans::contrast(refGrid.s, method = list(OLV_MPA_SP - OLV_REF_SP) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    0.0935
     SE  df z.ratio p.value
 0.0228 Inf 4.105   <.0001 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(OLV_MPA_BH - OLV_REF_BH) )
 contrast                                                       estimate    SE
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,      -1.1 0.203
  df z.ratio p.value
 Inf -5.410  <.0001 
 ```
 Quillback 
 ```
 > emmeans::contrast(refGrid.s, method = list(QBK_MPA_SP - QBK_REF_SP) )
 contrast                                                       estimate    SE
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.537 0.218
  df z.ratio p.value
 Inf 2.461   0.0138 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(QBK_MPA_BH - QBK_REF_BH) )
 contrast                                                       estimate    SE
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,     0.202 0.155
  df z.ratio p.value
 Inf 1.304   0.1924 
 ```
 Canary* (SP & BH)
 ```
 > emmeans::contrast(refGrid.s, method = list(CNY_MPA_SP - CNY_REF_SP) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.165
     SE  df z.ratio p.value
 0.0279 Inf 5.926   <.0001 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(CNY_MPA_BH - CNY_REF_BH) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.238
     SE  df z.ratio p.value
 0.0529 Inf 4.501   <.0001 
 ```
 Copper* (SP & BH)
 ```
 > emmeans::contrast(refGrid.s, method = list(CPR_MPA_SP - CPR_REF_SP) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.161
     SE  df z.ratio p.value
 0.0693 Inf 2.320   0.0203 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(CPR_MPA_BH - CPR_REF_BH) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    0.0984
     SE  df z.ratio p.value
 0.0334 Inf 2.943   0.0032
 ```
 Kelp Greenling* (BH)
 ```
 > emmeans::contrast(refGrid.s, method = list(KGL_MPA_SP - KGL_REF_SP) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.048
     SE  df z.ratio p.value
 0.0573 Inf 0.838   0.4019 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(KGL_MPA_BH - KGL_REF_BH) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,     0.257
     SE  df z.ratio p.value
 0.0322 Inf 7.966   <.0001 
 ```
 Rosy 
 ```
 > emmeans::contrast(refGrid.s, method = list(RSY_MPA_SP - RSY_REF_SP) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    0.0742
     SE  df z.ratio p.value
 0.0556 Inf 1.335   0.1819 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(RSY_MPA_BH - RSY_REF_BH) )
 contrast                                                       estimate SE df
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,    nonEst NA NA
 z.ratio p.value
 NA      NA 
 ```
 Yellow Eye* (SP)
 ```
 > emmeans::contrast(refGrid.s, method = list(YLE_MPA_SP - YLE_REF_SP) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.319
     SE  df z.ratio p.value
 0.0574 Inf 5.547   <.0001 
 ```
 Yellow Tail* (SP) 
 ```
 > emmeans::contrast(refGrid.s, method = list(YTL_MPA_SP - YTL_REF_SP) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0.131
     SE  df z.ratio p.value
 0.0193 Inf 6.786   <.0001 

Results are given on the log (not the response) scale. 
> emmeans::contrast(refGrid.s, method = list(YTL_MPA_BH - YTL_REF_BH) )
 contrast                                                       estimate
 c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    0.0488
     SE  df z.ratio p.value
 0.0478 Inf 1.020   0.3077 
 ```
 