###load libraries and read truncated indices data####
suppressPackageStartupMessages(library(tidyverse))
## Warning: package 'tidyverse' was built under R version 3.3.3
## Warning: package 'ggplot2' was built under R version 3.3.3
## Warning: package 'purrr' was built under R version 3.3.3
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
## 
##     date
library(broom)
library(forcats)
library(readr)
library(knitr)

setwd ("P:/MTM2 Oranga Taiao/EMP Data/EMP Data/Biotic Index calculator/Index calculation/MTM new Data")
source('theme_javier.R')
theme_set(theme_javier(8))
options(digits = 3)
all_ind<- read_csv('data_multiple_regression.csv',col_types =  cols(Date = 'D')) %>% 
  mutate(Type = factor(Type),
         logTOC = log(TOC+1))
options(digits = 3)

###data exploration####
##plot months
ggplot(all_ind, aes(x = factor(Month))) + geom_bar()

ggplot(all_ind, aes(x = Season)) + geom_bar()

##plot year
ggplot(all_ind, aes(x = year)) + geom_bar()

all_ind %>% 
  dplyr::select(AMBI:mud,Date) %>%
  gather(key,value,AMBI:TBI) %>% 
ggplot(.,aes(Date,value))+
  geom_point(alpha=.3)+
  geom_smooth()+
  facet_wrap(~key,scales = 'free')
## `geom_smooth()` using method = 'loess'

##grain size method
ggplot(all_ind, aes(x = grain.size.method)) + geom_bar()

all_ind %>% 
  dplyr::select(AMBI:TBI,-n,sqrtmud,grain.size.method) %>%
  gather(key,value,AMBI:TBI) %>% 
ggplot(.,aes(color= grain.size.method,x = sqrtmud,y = value))+
  geom_point()+
  geom_smooth(method='lm')+
  facet_wrap(~key,scales = 'free')

##core size and diameter
ggplot(all_ind, aes(x = factor(core.diameter))) + geom_bar()

all_ind %>% 
  dplyr::select(AMBI:TBI,-n,sqrtmud,core.diameter) %>%
  gather(key,value,AMBI:TBI) %>% 
  ggplot(.,aes(x = factor(core.diameter), y = value))+
  geom_boxplot()+
  facet_wrap(~key,scales = 'free')

ggplot(all_ind, aes(x = factor(core.depth))) + geom_bar()

all_ind %>% 
  dplyr::select(AMBI:TBI,-n,sqrtmud,core.depth) %>%
  gather(key,value,AMBI:TBI) %>% 
  ggplot(.,aes(x = factor(core.depth), y = value))+
  geom_boxplot()+
  facet_wrap(~key,scales = 'free')

##tides
ggplot(all_ind, aes(x = tidal.height)) + geom_bar()

all_ind %>% 
  dplyr::select(AMBI:TBI,-n,sqrtmud,tidal.height) %>%
  gather(key,value,AMBI:TBI) %>% 
  ggplot(.,aes(x = factor(tidal.height), y = value))+
  geom_boxplot()+
  facet_wrap(~key,scales = 'free')

##Typology
ggplot(all_ind, aes(x =  Type)) + geom_bar()

all_ind %>% 
  dplyr::select(AMBI:TBI,-n,sqrtmud,Type) %>%
  gather(key,value,AMBI:TBI) %>% 
  ggplot(.,aes(x = Type, y = value))+
  geom_boxplot()+
  facet_wrap(~key,scales = 'free')

all_ind %>% 
  dplyr::select(mud:TP,Type) %>%
  gather(key,value,mud:TP) %>% 
  ggplot(.,aes(x = Type, y = value))+
  geom_boxplot()+
  facet_wrap(~key,scales = 'free_y')+
  theme(axis.text.x = element_text(angle = 45,hjust = 1))
## Warning: Removed 368 rows containing non-finite values (stat_boxplot).

kable(table(all_ind$Type,all_ind$Region))
Cook Strait Eastern North Island North Eastern Southern Western North Island
7 14 48 22 41 0
8 20 0 67 22 3
##region
ggplot(all_ind, aes(x =  Region)) + geom_bar() +coord_flip()

all_ind %>% 
  dplyr::select(AMBI:TBI,-n,Region) %>%
  gather(key,value,AMBI:TBI) %>% 
  ggplot(.,aes(x = Region, y = value))+
  geom_boxplot()+
  facet_wrap(~key,scales = 'free_y')+
  theme(axis.text.x = element_text(angle = 45,hjust = 1))

all_ind %>% 
  dplyr::select(mud:TP,Region) %>%
  gather(key,value,mud:TP) %>% 
  ggplot(.,aes(x = Region, y = value))+
  geom_boxplot()+
  facet_wrap(~key,scales = 'free_y')+
  theme(axis.text.x = element_text(angle = 45,hjust = 1))
## Warning: Removed 368 rows containing non-finite values (stat_boxplot).

##Nitrogen
ggplot(all_ind, aes(x =  nitrogen.type)) + geom_bar()

all_ind %>% 
  dplyr::select(AMBI:TBI,-n,nitrogen.type) %>%
  gather(key,value,AMBI:TBI) %>% 
  ggplot(.,aes(x = nitrogen.type, y = value))+
  geom_boxplot()+
  facet_wrap(~key,scales = 'free_y')+
  theme(axis.text.x = element_text(angle = 45,hjust = 1))

##Taxonomy
ggplot(all_ind, aes(x =  taxonomy.by.)) + geom_bar() + coord_flip()

all_ind %>% 
  dplyr::select(AMBI:TBI,-n,taxonomy.by.) %>%
  gather(key,value,AMBI:TBI) %>% 
  ggplot(.,aes(x = taxonomy.by., y = value))+
  geom_boxplot()+
  facet_wrap(~key,scales = 'free_y')+
  theme(axis.text.x = element_text(angle = 45,hjust = 1))

###Geographic distrution of samples####
library(maps)
## 
## Attaching package: 'maps'
## The following object is masked from 'package:purrr':
## 
##     map
library(mapdata)
nz <- map_data("nzHires")

all_ind %>%
  dplyr::select(lat,lon,cesym,n,Type) %>%
  group_by(cesym) %>% 
  ggplot(aes(x = lon, y = lat)) +
  geom_polygon(
    data = nz,
    aes(x = long, y = lat, group = group),
    colour = "black",
    fill = 'grey80'
  ) +
  theme_bw() +
  ylab("Longitude") +
  xlab("Longitude") +
  coord_cartesian(
    xlim = c(165, 180),
    ylim = c(-47.5,-34),
    expand = FALSE
  ) +
  geom_point(aes(color = Type), 
             size= 2)
## Warning: Removed 98 rows containing missing values (geom_point).

###Environmental variables#####
barplot(colSums(is.na(all_ind[, 14:24])) / nrow(all_ind) * 100, 
        ylab = 'Proportion of missing values (%)',
        col = terrain.colors(10))

long_env <- 
  all_ind %>% 
  dplyr::select (logAFDW:metals) %>%
  gather(key, value) 

ggplot(long_env) +
  geom_histogram(aes(x = value), color = 1, fill = 'gray80') +
  xlab('') +
  facet_wrap( ~ key, scales = 'free')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 47 rows containing non-finite values (stat_bin).

cor_env <- cor(dplyr::select (all_ind,logAFDW:metals),use = 'pairwise.complete.obs')
print(cor_env,digits=3)
##         logAFDW  logPb  logCu sqrtZn sqrtmud  sqrtTP sqrtTN    logN metals
## logAFDW   1.000  0.348  0.181  0.387   0.614  0.3800  0.702  0.1620  0.378
## logPb     0.348  1.000  0.616  0.807   0.548  0.5754  0.397 -0.3604  0.832
## logCu     0.181  0.616  1.000  0.712   0.504  0.7955  0.277 -0.2601  0.766
## sqrtZn    0.387  0.807  0.712  1.000   0.427  0.7076  0.394 -0.2205  0.996
## sqrtmud   0.614  0.548  0.504  0.427   1.000  0.5156  0.629 -0.2480  0.461
## sqrtTP    0.380  0.575  0.796  0.708   0.516  1.0000  0.472 -0.0912  0.737
## sqrtTN    0.702  0.397  0.277  0.394   0.629  0.4715  1.000  0.1213  0.397
## logN      0.162 -0.360 -0.260 -0.221  -0.248 -0.0912  0.121  1.0000 -0.243
## metals    0.378  0.832  0.766  0.996   0.461  0.7373  0.397 -0.2433  1.000
corrplot::corrplot(
  cor_env,
  method = 'circle',
  type = 'lower',
  tl.cex = 1,
  tl.srt = 45,
  number.cex = .9,
  tl.col = 1,
  order = 'hclust'
)

####Indices####
##Unasigned taxa
ggplot(all_ind) +
  geom_histogram(aes(x = UA_R),
                 color = 1,
                 fill = 'gray80') +
  xlab('Unassigned taxa (%)')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

cor_indice <- cor(dplyr::select (all_ind,AMBI:TBI,logN,S,-n),use = 'pairwise.complete.obs')
print(cor_indice,digits=3)
##          AMBI AMBI_S  BENTIX  M_AMBI MEDOCC       N       S     BQI
## AMBI    1.000  0.942 -0.8344 -0.4583  0.994 -0.2019 -0.2023 -0.2483
## AMBI_S  0.942  1.000 -0.7538 -0.5260  0.924 -0.1594 -0.2654 -0.3216
## BENTIX -0.834 -0.754  1.0000  0.2747 -0.853  0.1503  0.0401  0.0868
## M_AMBI -0.458 -0.526  0.2747  1.0000 -0.427  0.0752  0.8868  0.8578
## MEDOCC  0.994  0.924 -0.8528 -0.4271  1.000 -0.1802 -0.1644 -0.2074
## N      -0.202 -0.159  0.1503  0.0752 -0.180  1.0000  0.3387  0.3159
## S      -0.202 -0.265  0.0401  0.8868 -0.164  0.3387  1.0000  0.9418
## BQI    -0.248 -0.322  0.0868  0.8578 -0.207  0.3159  0.9418  1.0000
## ITI    -0.141 -0.072  0.4921 -0.1085 -0.186 -0.0631 -0.1894 -0.2095
## TBI    -0.181 -0.253  0.0362  0.8223 -0.144  0.2746  0.9184  0.8846
## logN   -0.169 -0.149  0.0573  0.2981 -0.128  0.8527  0.5924  0.5978
##            ITI     TBI    logN
## AMBI   -0.1408 -0.1810 -0.1692
## AMBI_S -0.0720 -0.2530 -0.1489
## BENTIX  0.4921  0.0362  0.0573
## M_AMBI -0.1085  0.8223  0.2981
## MEDOCC -0.1855 -0.1440 -0.1282
## N      -0.0631  0.2746  0.8527
## S      -0.1894  0.9184  0.5924
## BQI    -0.2095  0.8846  0.5978
## ITI     1.0000 -0.2177 -0.1762
## TBI    -0.2177  1.0000  0.4871
## logN   -0.1762  0.4871  1.0000
corrplot::corrplot(
  cor_indice,
  method = 'circle',
  type = 'lower',
  tl.cex = 1,
  tl.srt = 45,
  number.cex = .9,
  tl.col = 1,
  order = 'hclust'
)

long_ind <-
  all_ind %>%
  gather(key, value, AMBI:S,BQI:TBI,logN,S,-N)

ggplot(long_ind) +
  geom_histogram(aes(x = value), color = 1, fill = 'gray80') +
  xlab('') +
  facet_wrap(~ key, scales = 'free')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

long_r_ind_env <- 
  all_ind %>% 
  gather(key, value, c(sqrtmud,sqrtTP,metals,logPb,logCu,sqrtZn,logAFDW,sqrtTN))

###AMBI
ggplot(long_r_ind_env,
       aes(x = value, y = AMBI, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('AMBI')
## Warning: Removed 47 rows containing non-finite values (stat_smooth).
## Warning: Removed 47 rows containing missing values (geom_point).

ggplot(long_r_ind_env,
       aes(x = value, y = AMBI, color = Region)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('AMBI')
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

##AMBI-S
ggplot(long_r_ind_env,
       aes(x = value, y = AMBI_S, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('AMBI-S')
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

ggplot(long_r_ind_env,
       aes(x = value, y = AMBI_S, color = Region)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('AMBI-S')
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

##M-AMBI
ggplot(long_r_ind_env,
       aes(x = value, y = M_AMBI, color =  Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('M-AMBI') 
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

ggplot(long_r_ind_env,
       aes(x = value, y = M_AMBI, color =  Region)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('M-AMBI') 
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

##BENTIX
ggplot(long_r_ind_env,
       aes(x = value, y = BENTIX, color =  Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('BENTIX') 
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

ggplot(long_r_ind_env,
       aes(x = value, y = BENTIX, color =  Region)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('BENTIX') 
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

##MEDOCC
ggplot(long_r_ind_env,
       aes(x = value, y = MEDOCC, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('MEDOCC') 
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

ggplot(long_r_ind_env,
       aes(x = value, y = MEDOCC, color = Region)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('MEDOCC') 
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

###logN
ggplot(long_r_ind_env,
       aes(x = value, y = log(N), color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('logN') 
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

ggplot(long_r_ind_env,
       aes(x = value, y = log(N), color = Region)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('logN') 
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

###richness
ggplot(long_r_ind_env,
       aes(x = value, y = S, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('Richness') 
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

ggplot(long_r_ind_env,
       aes(x = value, y = S, color = Region)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('Richness') 
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

###BQI
ggplot(long_r_ind_env,
       aes(x = value, y = BQI, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('BQI') 
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

ggplot(long_r_ind_env,
       aes(x = value, y = BQI, color = Region)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('BQI') 
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

###ITI
ggplot(long_r_ind_env,
       aes(x = value, y = ITI, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('ITI')
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

ggplot(long_r_ind_env,
       aes(x = value, y = ITI, color = Region)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('ITI')
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

###TBI
ggplot(long_r_ind_env,
       aes(x = value, y = TBI, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('TBI') 
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

ggplot(long_r_ind_env,
       aes(x = value, y = TBI, color = Region)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap( ~ key, scale = 'free') +
  xlab('') +
  ylab('TBI') 
## Warning: Removed 47 rows containing non-finite values (stat_smooth).

## Warning: Removed 47 rows containing missing values (geom_point).

###Multiple stressors models####
dat_mult <- 
  all_ind %>% 
  dplyr::select(Region,n,AMBI:MEDOCC,BQI:TBI,logN,S,Type,metals,sqrtTP,sqrtmud) %>% 
  drop_na()

###VIF####
source('corvif.R')
corvif(dplyr::select(dat_mult,metals,sqrtTP,sqrtmud,Type,Region))
## 
## 
## Variance inflation factors
## 
##         GVIF Df GVIF^(1/2Df)
## metals  2.78  1         1.67
## sqrtTP  2.67  1         1.63
## sqrtmud 1.70  1         1.30
## Type    1.58  1         1.26
## Region  2.50  4         1.12
EG_models <- 
  dat_mult %>%
  gather(index, value,AMBI:MEDOCC,BQI:TBI,logN,S) %>%
  group_by(index) %>%
  do(fit_index = step(lm(value ~  (metals + sqrtmud + sqrtTP)* Type + Region,
      data = .,
      weights = sqrt(n)),
    direction = 'backward',
    trace=0))

index_coef_all  <- 
  tidy(EG_models, fit_index,conf.int = T) %>% 
  data.frame(.) %>%
  filter(term!='(Intercept)') %>% 
  mutate(term = factor(term))

kable(index_coef_all)
index term estimate std.error statistic p.value conf.low conf.high
AMBI metals 0.033 0.016 2.045 0.042 0.001 0.066
AMBI sqrtmud 0.118 0.020 5.868 0.000 0.078 0.157
AMBI Type8 0.291 0.075 3.905 0.000 0.144 0.438
AMBI_S metals 0.046 0.016 2.907 0.004 0.015 0.076
AMBI_S sqrtmud 0.174 0.024 7.202 0.000 0.126 0.221
AMBI_S sqrtTP -0.010 0.012 -0.811 0.418 -0.035 0.014
AMBI_S Type8 0.106 0.226 0.469 0.640 -0.339 0.550
AMBI_S RegionEastern North Island -0.240 0.102 -2.350 0.020 -0.442 -0.039
AMBI_S RegionNorth Eastern 0.013 0.080 0.162 0.871 -0.145 0.171
AMBI_S RegionSouthern -0.011 0.082 -0.131 0.896 -0.171 0.150
AMBI_S RegionWestern North Island 0.166 0.211 0.787 0.432 -0.250 0.583
AMBI_S sqrtmud:Type8 -0.112 0.034 -3.302 0.001 -0.178 -0.045
AMBI_S sqrtTP:Type8 0.021 0.014 1.477 0.141 -0.007 0.049
BENTIX metals -0.054 0.027 -2.002 0.046 -0.106 -0.001
BENTIX sqrtmud -0.165 0.034 -4.841 0.000 -0.231 -0.098
BENTIX Type8 -0.285 0.140 -2.040 0.042 -0.560 -0.010
BENTIX RegionEastern North Island 0.582 0.202 2.879 0.004 0.184 0.980
BENTIX RegionNorth Eastern -0.070 0.179 -0.391 0.696 -0.422 0.282
BENTIX RegionSouthern -0.118 0.179 -0.659 0.511 -0.470 0.234
BENTIX RegionWestern North Island -0.385 0.464 -0.830 0.408 -1.298 0.529
BQI metals -0.123 0.043 -2.822 0.005 -0.208 -0.037
BQI sqrtmud -0.311 0.046 -6.732 0.000 -0.402 -0.220
BQI sqrtTP 0.025 0.033 0.756 0.450 -0.040 0.089
BQI Type8 -2.540 0.623 -4.075 0.000 -3.768 -1.312
BQI RegionEastern North Island -1.719 0.263 -6.546 0.000 -2.236 -1.201
BQI RegionNorth Eastern 0.857 0.224 3.830 0.000 0.416 1.298
BQI RegionSouthern -1.145 0.223 -5.126 0.000 -1.585 -0.705
BQI RegionWestern North Island 0.951 0.577 1.647 0.101 -0.187 2.088
BQI sqrtTP:Type8 0.098 0.034 2.936 0.004 0.032 0.164
ITI sqrtmud -1.718 0.466 -3.687 0.000 -2.637 -0.800
ITI sqrtTP 0.481 0.195 2.467 0.014 0.097 0.864
ITI Type8 -3.142 1.741 -1.805 0.072 -6.573 0.288
ITI RegionEastern North Island 11.084 2.488 4.455 0.000 6.182 15.987
ITI RegionNorth Eastern -1.691 2.246 -0.753 0.452 -6.116 2.735
ITI RegionSouthern -3.600 2.226 -1.617 0.107 -7.985 0.786
ITI RegionWestern North Island -2.303 5.802 -0.397 0.692 -13.736 9.130
logN sqrtmud 0.111 0.041 2.677 0.008 0.029 0.193
logN sqrtTP -0.029 0.019 -1.553 0.122 -0.066 0.008
logN Type8 -1.969 0.390 -5.055 0.000 -2.737 -1.202
logN RegionEastern North Island -1.089 0.170 -6.400 0.000 -1.424 -0.754
logN RegionNorth Eastern 1.393 0.139 10.044 0.000 1.119 1.666
logN RegionSouthern 0.605 0.140 4.326 0.000 0.329 0.880
logN RegionWestern North Island 0.750 0.365 2.053 0.041 0.030 1.469
logN sqrtmud:Type8 -0.202 0.058 -3.491 0.001 -0.316 -0.088
logN sqrtTP:Type8 0.104 0.024 4.275 0.000 0.056 0.153
M_AMBI metals -0.009 0.003 -2.730 0.007 -0.015 -0.002
M_AMBI sqrtmud -0.025 0.003 -7.121 0.000 -0.032 -0.018
M_AMBI sqrtTP 0.000 0.002 0.190 0.849 -0.004 0.005
M_AMBI Type8 -0.145 0.047 -3.100 0.002 -0.237 -0.053
M_AMBI RegionEastern North Island -0.076 0.020 -3.861 0.000 -0.115 -0.037
M_AMBI RegionNorth Eastern 0.023 0.017 1.363 0.174 -0.010 0.056
M_AMBI RegionSouthern -0.075 0.017 -4.460 0.000 -0.108 -0.042
M_AMBI RegionWestern North Island 0.040 0.043 0.916 0.361 -0.046 0.125
M_AMBI sqrtTP:Type8 0.005 0.003 2.105 0.036 0.000 0.010
MEDOCC sqrtmud 0.169 0.023 7.426 0.000 0.124 0.214
MEDOCC Type8 0.381 0.093 4.073 0.000 0.196 0.565
S metals -0.266 0.124 -2.141 0.033 -0.511 -0.021
S sqrtmud -0.581 0.191 -3.041 0.003 -0.958 -0.205
S sqrtTP 0.003 0.099 0.032 0.974 -0.192 0.198
S Type8 -8.189 1.788 -4.580 0.000 -11.712 -4.666
S RegionEastern North Island -5.133 0.810 -6.339 0.000 -6.729 -3.538
S RegionNorth Eastern 3.270 0.636 5.139 0.000 2.016 4.524
S RegionSouthern -2.273 0.646 -3.518 0.001 -3.546 -1.000
S RegionWestern North Island 2.322 1.676 1.386 0.167 -0.980 5.624
S sqrtmud:Type8 -0.428 0.268 -1.598 0.111 -0.955 0.100
S sqrtTP:Type8 0.385 0.112 3.431 0.001 0.164 0.607
TBI metals -0.016 0.003 -4.595 0.000 -0.022 -0.009
TBI sqrtmud -0.021 0.003 -6.363 0.000 -0.027 -0.014
TBI sqrtTP 0.007 0.002 3.852 0.000 0.003 0.011
TBI Type8 -0.026 0.013 -1.982 0.049 -0.052 0.000
TBI RegionEastern North Island -0.108 0.019 -5.618 0.000 -0.146 -0.070
TBI RegionNorth Eastern 0.073 0.016 4.660 0.000 0.042 0.104
TBI RegionSouthern -0.083 0.016 -5.242 0.000 -0.114 -0.052
TBI RegionWestern North Island 0.000 0.041 -0.009 0.992 -0.081 0.080
TBI metals:Type8 0.011 0.005 2.353 0.020 0.002 0.020
##Model summary
index_Summ  <-
  glance(EG_models, fit_index) 

kable(index_Summ)
index r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
AMBI 0.257 0.248 0.889 26.88 0 4 -196 402 419 184.14 233
AMBI_S 0.401 0.375 0.615 15.15 0 11 -105 234 276 85.54 226
BENTIX 0.231 0.208 1.384 9.84 0 8 -299 615 647 438.55 229
BQI 0.557 0.540 1.716 31.73 0 10 -349 719 757 668.37 227
ITI 0.278 0.256 17.318 12.60 0 8 -898 1813 1844 68679.25 229
logN 0.571 0.554 1.062 33.52 0 10 -235 492 530 256.23 227
M_AMBI 0.432 0.409 0.129 19.15 0 10 265 -508 -470 3.77 227
MEDOCC 0.240 0.233 1.160 36.92 0 3 -259 527 541 314.95 234
S 0.553 0.533 4.875 27.94 0 11 -596 1215 1257 5371.88 226
TBI 0.579 0.562 0.121 34.62 0 10 279 -537 -499 3.34 227
ggplot(index_coef_all,
       aes(
         x = estimate,
         y = fct_rev(fct_inorder(term)),
         color = term
       )) +
  geom_point() +
  geom_errorbarh(aes(
    xmin = conf.low,
    xmax = conf.high,
    height = .2
  )) +
  facet_wrap( ~ index,scales = 'free_x') +
  geom_vline(xintercept = 0,
             lty = 2,
             col = 'gray60') +
  ylab('') +
  xlab('Coefficients') +
  scale_color_discrete(guide = F)

###residuals
index_residuals  <- 
  augment(EG_models, fit_index)

ggplot(index_residuals)+
  geom_point(aes(x=.fitted,y=.resid),alpha = .3)+
  facet_wrap(~index,scale='free')+
  geom_hline(yintercept = 0,lty=2,col=2)

ggplot(index_residuals)+
  stat_qq(aes(sample= .std.resid),alpha = .3)+
  facet_wrap(~index,scale='free')+
  geom_abline(intercept =  0, slope = 1, lty=2,col=2)

conf.level <- 0.95
len.data <- length(dat_mult) 
conf.int <- qnorm((1 + conf.level) / 2) / sqrt(len.data)

index_residuals %>% 
  group_by(index) %>%
  do(acf1 = acf(index_residuals$.std.resid),plot = F) %>% 
  tidy(acf1) %>% 
  ggplot(., aes(x=lag, y=acf)) + 
  geom_bar(stat='identity', width=0.1) +
  geom_hline(yintercept = conf.int,
               color='blue', linetype='dashed') +
  facet_wrap(~index)+
  geom_hline(yintercept = 0,
             color=1)

####Plots of relationship with truncated data####
long_dat <- 
  dat_mult %>% 
  gather(key, value, c(sqrtmud,sqrtTP,metals))

###AMBI
ggplot(long_dat,
       aes(x = value, y = AMBI, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap(~ key, scale = 'free', ncol = 1) +
  xlab('') +
  ylab('AMBI')

###AMBI_S
ggplot(long_dat,
       aes(x = value, y = AMBI_S, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap(~ key, scale = 'free', ncol = 1) +
  xlab('') +
  ylab('AMBI-S')

###M_AMBI
ggplot(long_dat,
       aes(x = value, y = M_AMBI, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap(~ key, scale = 'free', ncol = 1) +
  xlab('') +
  ylab('M-AMBI')

##BENTIX
ggplot(long_dat,
       aes(x = value, y = BENTIX, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap(~ key, scale = 'free', ncol = 1) +
  xlab('') +
  ylab('BENTIX')

###TBI
ggplot(long_dat,
       aes(x = value, y = TBI, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap(~ key, scale = 'free', ncol = 1) +
  xlab('') +
  ylab('TBI')

###BQI
ggplot(long_dat,
       aes(x = value, y = BQI, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap(~ key, scale = 'free', ncol = 1) +
  xlab('') +
  ylab('BQI')

###MEDOCC
ggplot(long_dat,
       aes(x = value, y = MEDOCC, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap(~ key, scale = 'free', ncol = 1) +
  xlab('') +
  ylab('MEDOCC')

###TBI
ggplot(long_dat,
       aes(x = value, y = TBI, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap(~ key, scale = 'free', ncol = 1) +
  xlab('') +
  ylab('TBI')

###logN
ggplot(long_dat,
       aes(x = value, y = logN, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap(~ key, scale = 'free', ncol = 1) +
  xlab('') +
  ylab('logN')

###logN
ggplot(long_dat,
       aes(x = value, y = logN, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap(~ key, scale = 'free', ncol = 1) +
  xlab('') +
  ylab('logN')

###S
ggplot(long_dat,
       aes(x = value, y = S, color = Type)) +
  geom_point(alpha  = .3) +
  geom_smooth(method = 'lm') +
  facet_wrap(~ key, scale = 'free', ncol = 1) +
  xlab('') +
  ylab('Richness')