###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))
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)
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)
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')
