multiple_regressions_options.R

###Multiple stressors models####
###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(broom)
library(knitr)
library(forcats)

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('all_indices_truncated.csv', col_types =  cols(Date = 'D')) %>%
  mutate(Type = factor(Type))

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


###Model 1 interaction Type and stressors####
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
  )

###Model 2 interaction stressor vs Region and Type####
EG_models2 <-
  dat_mult %>%
  gather(index, value, AMBI:MEDOCC, BQI:TBI, logN, S) %>%
  group_by(index) %>%
  do(
    fit_index2 = lm(
      value ~  (metals + sqrtmud + sqrtTP) * Type + (metals + sqrtmud + sqrtTP) * Region,
      data = .,
      weights = sqrt(n)
    ),
    direction = 'backward',
    trace = 0
  )

index_coef_all2  <-
  tidy(EG_models2, fit_index2) %>%
  data.frame(.) %>%
  filter(term != '(Intercept)') %>%
  mutate(term = factor(term))

##Model summary
index_Summ2  <-
  glance(EG_models2, fit_index2)
kable(index_Summ2)

index	r.squared	adj.r.squared	sigma	statistic	p.value	df	logLik	AIC	BIC	deviance	df.residual
AMBI	0.351	0.284	0.867	5.25	0	23	-180	408	491	160.93	214
AMBI_S	0.454	0.398	0.604	8.10	0	23	-94	236	319	77.96	214
BENTIX	0.370	0.305	1.296	5.71	0	23	-275	598	681	359.38	214
BQI	0.582	0.539	1.718	13.53	0	23	-342	732	815	631.29	214
ITI	0.403	0.342	16.285	6.58	0	23	-875	1798	1881	56754.52	214
logN	0.635	0.597	1.009	16.91	0	23	-216	480	563	217.97	214
M_AMBI	0.482	0.429	0.127	9.07	0	23	276	-504	-421	3.43	214
MEDOCC	0.339	0.271	1.132	4.98	0	23	-243	534	617	273.98	214
S	0.574	0.530	4.890	13.11	0	23	-590	1228	1311	5116.24	214
TBI	0.593	0.552	0.123	14.19	0	23	284	-519	-436	3.22	214

ggplot(index_coef_all2,
       aes(
         x = estimate,
         y = fct_rev(fct_inorder(term)),
         color = term
       )) +
  geom_point() +
  geom_errorbarh(aes(
    xmin = estimate - std.error * 1.96,
    xmax = estimate + std.error * 1.96,
    height = .2
  )) +
  facet_wrap(~ index, scales = 'free_x') +
  geom_vline(xintercept = 0,
             lty = 2,
             col = 'gray60') +
  ylab('') +
  xlab('Coefficients') +
  scale_color_discrete(guide = F) +
  theme_javier()

###residuals
index_residuals2  <-
  augment(EG_models2, fit_index2)

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

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

## Warning: Removed 30 rows containing non-finite values (stat_qq).

###Model 3 Stressors only####
EG_models3 <-
  dat_mult %>%
  gather(index, value, AMBI:MEDOCC, BQI:TBI, logN, S) %>%
  group_by(index) %>%
  do(fit_index3 = step(
    lm(
      value ~  metals + sqrtmud + sqrtTP,
      data = .,
      weights = sqrt(n)
    ),
    direction = 'backward',
    trace = 0
  ))

index_coef_all3  <-
  tidy(EG_models3, fit_index3, conf.int = T) %>%
  data.frame(.) %>%
  filter(term != '(Intercept)') %>%
  mutate(term = factor(term))

kable(index_coef_all3)

index	term	estimate	std.error	statistic	p.value	conf.low	conf.high
AMBI	sqrtmud	0.140	0.018	7.79	0.000	0.105	0.176
AMBI_S	metals	0.025	0.011	2.20	0.028	0.003	0.047
AMBI_S	sqrtmud	0.118	0.014	8.30	0.000	0.090	0.146
BENTIX	sqrtmud	-0.156	0.029	-5.39	0.000	-0.213	-0.099
BQI	metals	-0.193	0.050	-3.84	0.000	-0.291	-0.094
BQI	sqrtmud	-0.291	0.054	-5.43	0.000	-0.397	-0.186
BQI	sqrtTP	0.069	0.031	2.24	0.026	0.008	0.130
ITI	metals	1.626	0.340	4.78	0.000	0.956	2.297
ITI	sqrtmud	-1.099	0.429	-2.56	0.011	-1.944	-0.254
logN	metals	-0.122	0.034	-3.60	0.000	-0.189	-0.055
logN	sqrtmud	-0.109	0.036	-2.99	0.003	-0.180	-0.037
logN	sqrtTP	0.056	0.021	2.69	0.008	0.015	0.097
M_AMBI	metals	-0.007	0.003	-2.75	0.006	-0.012	-0.002
M_AMBI	sqrtmud	-0.021	0.003	-6.52	0.000	-0.027	-0.015
MEDOCC	sqrtmud	0.172	0.024	7.33	0.000	0.126	0.219
S	metals	-0.524	0.141	-3.71	0.000	-0.802	-0.246
S	sqrtmud	-0.872	0.151	-5.77	0.000	-1.170	-0.574
S	sqrtTP	0.212	0.087	2.44	0.015	0.041	0.384
TBI	metals	-0.017	0.004	-4.74	0.000	-0.024	-0.010
TBI	sqrtmud	-0.019	0.004	-4.75	0.000	-0.026	-0.011
TBI	sqrtTP	0.006	0.002	2.79	0.006	0.002	0.011

##Model summary
index_Summ3  <-
  glance(EG_models3, fit_index3)
kable(index_Summ3)

index	r.squared	adj.r.squared	sigma	statistic	p.value	df	logLik	AIC	BIC	deviance	df.residual
AMBI	0.205	0.202	0.916	60.70	0	2	-204	414	424	196.98	235
AMBI_S	0.334	0.328	0.638	58.60	0	3	-118	243	257	95.19	234
BENTIX	0.110	0.106	1.470	29.01	0	2	-316	638	648	507.81	235
BQI	0.226	0.216	2.240	22.62	0	4	-415	840	857	1168.82	233
ITI	0.089	0.082	19.242	11.48	0	3	-925	1858	1872	86636.14	234
logN	0.106	0.094	1.514	9.17	0	4	-322	654	671	533.77	233
M_AMBI	0.269	0.263	0.144	43.13	0	3	235	-463	-449	4.84	234
MEDOCC	0.186	0.182	1.198	53.68	0	2	-268	541	552	337.28	235
S	0.229	0.219	6.304	23.09	0	4	-660	1330	1348	9259.05	233
TBI	0.221	0.211	0.163	21.99	0	4	206	-403	-386	6.17	233

ggplot(index_coef_all3,
       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_residuals3  <-
  augment(EG_models3, fit_index3)

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

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

#######Model 4 stressor models no interactions ####
EG_models4 <-
  dat_mult %>%
  gather(index, value, AMBI:MEDOCC, BQI:TBI, logN, S) %>%
  group_by(index) %>%
  do(fit_index4 = step(
    lm(
      value ~  (metals + sqrtmud + sqrtTP) ^ 2,
      data = .,
      weights = sqrt(n)
    ),
    direction = 'backward',
    trace = 0
  ))

index_coef_all4  <-
  tidy(EG_models4, fit_index4, conf.int = T) %>%
  data.frame(.) %>%
  filter(term != '(Intercept)') %>%
  mutate(term = factor(term))

kable(index_coef_all4)

index	term	estimate	std.error	statistic	p.value	conf.low	conf.high
AMBI	metals	0.008	0.051	0.155	0.877	-0.092	0.108
AMBI	sqrtmud	0.123	0.022	5.678	0.000	0.081	0.166
AMBI	sqrtTP	0.001	0.013	0.055	0.956	-0.024	0.026
AMBI	metals:sqrtmud	0.035	0.011	3.054	0.003	0.012	0.058
AMBI	metals:sqrtTP	-0.006	0.004	-1.624	0.106	-0.013	0.001
AMBI_S	metals	0.023	0.035	0.659	0.511	-0.046	0.092
AMBI_S	sqrtmud	0.111	0.015	7.409	0.000	0.081	0.140
AMBI_S	sqrtTP	0.002	0.009	0.181	0.857	-0.016	0.019
AMBI_S	metals:sqrtmud	0.029	0.008	3.721	0.000	0.014	0.045
AMBI_S	metals:sqrtTP	-0.005	0.002	-2.100	0.037	-0.010	0.000
BENTIX	metals	0.099	0.048	2.049	0.042	0.004	0.194
BENTIX	sqrtmud	-0.158	0.033	-4.840	0.000	-0.222	-0.094
BENTIX	metals:sqrtmud	-0.026	0.012	-2.219	0.027	-0.049	-0.003
BQI	metals	-0.193	0.050	-3.842	0.000	-0.291	-0.094
BQI	sqrtmud	-0.291	0.054	-5.427	0.000	-0.397	-0.186
BQI	sqrtTP	0.069	0.031	2.242	0.026	0.008	0.130
ITI	metals	0.719	0.632	1.138	0.256	-0.526	1.964
ITI	sqrtmud	-1.129	0.428	-2.640	0.009	-1.971	-0.286
ITI	metals:sqrtmud	0.264	0.155	1.702	0.090	-0.042	0.569
logN	metals	0.097	0.067	1.442	0.151	-0.036	0.230
logN	sqrtmud	-0.752	0.167	-4.493	0.000	-1.082	-0.422
logN	sqrtTP	-0.052	0.035	-1.478	0.141	-0.122	0.017
logN	metals:sqrtmud	-0.064	0.018	-3.640	0.000	-0.099	-0.030
logN	sqrtmud:sqrtTP	0.032	0.008	3.925	0.000	0.016	0.048
M_AMBI	metals	-0.007	0.003	-2.755	0.006	-0.012	-0.002
M_AMBI	sqrtmud	-0.021	0.003	-6.520	0.000	-0.027	-0.015
MEDOCC	metals	0.006	0.067	0.089	0.930	-0.126	0.138
MEDOCC	sqrtmud	0.159	0.029	5.539	0.000	0.102	0.215
MEDOCC	sqrtTP	0.004	0.017	0.251	0.802	-0.029	0.037
MEDOCC	metals:sqrtmud	0.041	0.015	2.751	0.006	0.012	0.071
MEDOCC	metals:sqrtTP	-0.008	0.005	-1.618	0.107	-0.017	0.002
S	metals	-0.524	0.141	-3.715	0.000	-0.802	-0.246
S	sqrtmud	-0.872	0.151	-5.769	0.000	-1.170	-0.574
S	sqrtTP	0.212	0.087	2.440	0.015	0.041	0.384
TBI	metals	-0.028	0.009	-3.091	0.002	-0.046	-0.010
TBI	sqrtmud	-0.018	0.004	-4.559	0.000	-0.025	-0.010
TBI	sqrtTP	0.007	0.002	3.072	0.002	0.003	0.011
TBI	metals:sqrtmud	-0.004	0.002	-2.178	0.030	-0.009	0.000
TBI	metals:sqrtTP	0.001	0.001	2.049	0.042	0.000	0.003

##Model summary multiple stressor interactions####
index_Summ4  <-
  glance(EG_models4, fit_index4)
kable(index_Summ4)

index	r.squared	adj.r.squared	sigma	statistic	p.value	df	logLik	AIC	BIC	deviance	df.residual
AMBI	0.244	0.228	0.900	14.94	0	6	-198	410	434	187.28	231
AMBI_S	0.377	0.363	0.621	27.95	0	6	-110	233	258	89.02	231
BENTIX	0.129	0.117	1.461	11.47	0	4	-314	637	654	497.07	233
BQI	0.226	0.216	2.240	22.62	0	4	-415	840	857	1168.82	233
ITI	0.101	0.089	19.164	8.68	0	4	-924	1857	1875	85572.35	233
logN	0.166	0.148	1.468	9.18	0	6	-314	641	666	497.84	231
M_AMBI	0.269	0.263	0.144	43.13	0	3	235	-463	-449	4.84	234
MEDOCC	0.216	0.199	1.186	12.71	0	6	-263	540	565	324.91	231
S	0.229	0.219	6.304	23.09	0	4	-660	1330	1348	9259.05	233
TBI	0.238	0.221	0.162	14.39	0	6	209	-404	-380	6.04	231

ggplot(index_coef_all4,
       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_residuals4  <-
  augment(EG_models4, fit_index4)

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

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

##Model summary multiple stressor interactions####
index_Summ4  <-
  glance(EG_models4, fit_index4)
kable(index_Summ4)

index	r.squared	adj.r.squared	sigma	statistic	p.value	df	logLik	AIC	BIC	deviance	df.residual
AMBI	0.244	0.228	0.900	14.94	0	6	-198	410	434	187.28	231
AMBI_S	0.377	0.363	0.621	27.95	0	6	-110	233	258	89.02	231
BENTIX	0.129	0.117	1.461	11.47	0	4	-314	637	654	497.07	233
BQI	0.226	0.216	2.240	22.62	0	4	-415	840	857	1168.82	233
ITI	0.101	0.089	19.164	8.68	0	4	-924	1857	1875	85572.35	233
logN	0.166	0.148	1.468	9.18	0	6	-314	641	666	497.84	231
M_AMBI	0.269	0.263	0.144	43.13	0	3	235	-463	-449	4.84	234
MEDOCC	0.216	0.199	1.186	12.71	0	6	-263	540	565	324.91	231
S	0.229	0.219	6.304	23.09	0	4	-660	1330	1348	9259.05	233
TBI	0.238	0.221	0.162	14.39	0	6	209	-404	-380	6.04	231

ggplot(index_coef_all4,
       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_residuals4  <-
  augment(EG_models4, fit_index4)

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

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

#######Model 5 stressor models no interactions, type interaction and  region####
EG_models5 <-
  dat_mult %>%
  gather(index, value, AMBI:MEDOCC, BQI:TBI, logN, S) %>%
  group_by(index) %>%
  do(fit_index5 = step(
    lm(
      value ~  (metals + sqrtmud + sqrtTP) ^ 2 +  (metals + sqrtmud + sqrtTP) * Type + Region ,
      data = .,
      weights = sqrt(n)
    ),
    direction = 'backward',
    trace = 0
  ))

index_coef_all5  <-
  tidy(EG_models5, fit_index5, conf.int = T) %>%
  data.frame(.) %>%
  filter(term != '(Intercept)') %>%
  mutate(term = factor(term))

kable(index_coef_all5)

index	term	estimate	std.error	statistic	p.value	conf.low	conf.high
AMBI	metals	0.026	0.062	0.409	0.683	-0.097	0.148
AMBI	sqrtmud	0.441	0.120	3.660	0.000	0.203	0.678
AMBI	sqrtTP	0.012	0.022	0.542	0.588	-0.032	0.056
AMBI	Type8	0.109	0.335	0.327	0.744	-0.550	0.769
AMBI	RegionEastern North Island	-0.473	0.148	-3.202	0.002	-0.764	-0.182
AMBI	RegionNorth Eastern	-0.027	0.116	-0.231	0.817	-0.255	0.202
AMBI	RegionSouthern	-0.095	0.116	-0.826	0.409	-0.323	0.132
AMBI	RegionWestern North Island	-0.140	0.298	-0.469	0.639	-0.728	0.448
AMBI	metals:sqrtmud	0.065	0.014	4.662	0.000	0.038	0.093
AMBI	metals:sqrtTP	-0.009	0.004	-2.475	0.014	-0.017	-0.002
AMBI	sqrtmud:sqrtTP	-0.012	0.005	-2.246	0.026	-0.023	-0.001
AMBI	sqrtmud:Type8	-0.135	0.048	-2.834	0.005	-0.229	-0.041
AMBI	sqrtTP:Type8	0.031	0.021	1.477	0.141	-0.010	0.072
AMBI_S	metals	0.049	0.043	1.138	0.256	-0.036	0.133
AMBI_S	sqrtmud	0.337	0.082	4.098	0.000	0.175	0.500
AMBI_S	sqrtTP	0.002	0.015	0.130	0.897	-0.028	0.032
AMBI_S	Type8	0.065	0.229	0.284	0.777	-0.386	0.516
AMBI_S	RegionEastern North Island	-0.345	0.101	-3.424	0.001	-0.544	-0.147
AMBI_S	RegionNorth Eastern	0.046	0.079	0.578	0.564	-0.110	0.202
AMBI_S	RegionSouthern	-0.009	0.079	-0.113	0.910	-0.165	0.147
AMBI_S	RegionWestern North Island	-0.032	0.204	-0.158	0.874	-0.434	0.370
AMBI_S	metals:sqrtmud	0.050	0.010	5.179	0.000	0.031	0.068
AMBI_S	metals:sqrtTP	-0.008	0.003	-3.100	0.002	-0.013	-0.003
AMBI_S	sqrtmud:sqrtTP	-0.007	0.004	-2.006	0.046	-0.015	0.000
AMBI_S	sqrtmud:Type8	-0.133	0.033	-4.093	0.000	-0.198	-0.069
AMBI_S	sqrtTP:Type8	0.027	0.014	1.917	0.057	-0.001	0.056
BENTIX	metals	-0.012	0.098	-0.120	0.905	-0.205	0.181
BENTIX	sqrtmud	-0.735	0.189	-3.890	0.000	-1.107	-0.362
BENTIX	sqrtTP	-0.004	0.035	-0.108	0.914	-0.072	0.065
BENTIX	Type8	0.408	0.525	0.777	0.438	-0.626	1.442
BENTIX	RegionEastern North Island	0.945	0.231	4.084	0.000	0.489	1.402
BENTIX	RegionNorth Eastern	-0.134	0.182	-0.738	0.461	-0.493	0.224
BENTIX	RegionSouthern	-0.120	0.181	-0.660	0.510	-0.477	0.237
BENTIX	RegionWestern North Island	-0.162	0.468	-0.345	0.730	-1.083	0.760
BENTIX	metals:sqrtmud	-0.094	0.022	-4.272	0.000	-0.137	-0.051
BENTIX	metals:sqrtTP	0.011	0.006	1.902	0.058	0.000	0.023
BENTIX	sqrtmud:sqrtTP	0.023	0.008	2.705	0.007	0.006	0.040
BENTIX	sqrtmud:Type8	0.165	0.075	2.211	0.028	0.018	0.313
BENTIX	sqrtTP:Type8	-0.068	0.033	-2.064	0.040	-0.133	-0.003
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	-4.977	1.450	-3.433	0.001	-7.833	-2.121
ITI	sqrtTP	-0.135	0.323	-0.416	0.678	-0.771	0.502
ITI	Type8	-3.422	1.728	-1.980	0.049	-6.826	-0.017
ITI	RegionEastern North Island	11.118	2.463	4.513	0.000	6.264	15.972
ITI	RegionNorth Eastern	-2.377	2.242	-1.060	0.290	-6.795	2.042
ITI	RegionSouthern	-4.665	2.249	-2.074	0.039	-9.096	-0.233
ITI	RegionWestern North Island	-3.375	5.763	-0.586	0.559	-14.729	7.980
ITI	sqrtmud:sqrtTP	0.161	0.068	2.371	0.019	0.027	0.294
logN	metals	0.111	0.062	1.774	0.077	-0.012	0.234
logN	sqrtmud	0.103	0.042	2.479	0.014	0.021	0.185
logN	sqrtTP	-0.019	0.022	-0.862	0.390	-0.062	0.024
logN	Type8	-1.798	0.397	-4.527	0.000	-2.581	-1.015
logN	RegionEastern North Island	-1.035	0.178	-5.830	0.000	-1.385	-0.685
logN	RegionNorth Eastern	1.466	0.143	10.260	0.000	1.185	1.748
logN	RegionSouthern	0.651	0.142	4.579	0.000	0.371	0.932
logN	RegionWestern North Island	0.768	0.364	2.112	0.036	0.052	1.485
logN	metals:sqrtTP	-0.006	0.003	-2.003	0.046	-0.012	0.000
logN	sqrtmud:Type8	-0.197	0.058	-3.396	0.001	-0.312	-0.083
logN	sqrtTP:Type8	0.096	0.025	3.885	0.000	0.047	0.145
M_AMBI	metals	-0.012	0.004	-3.344	0.001	-0.019	-0.005
M_AMBI	sqrtmud	-0.025	0.003	-7.283	0.000	-0.032	-0.019
M_AMBI	sqrtTP	0.003	0.002	1.460	0.146	-0.001	0.007
M_AMBI	Type8	-0.039	0.014	-2.801	0.006	-0.067	-0.012
M_AMBI	RegionEastern North Island	-0.072	0.020	-3.530	0.001	-0.112	-0.032
M_AMBI	RegionNorth Eastern	0.019	0.017	1.163	0.246	-0.014	0.052
M_AMBI	RegionSouthern	-0.075	0.017	-4.466	0.000	-0.108	-0.042
M_AMBI	RegionWestern North Island	0.039	0.043	0.909	0.365	-0.046	0.125
M_AMBI	metals:Type8	0.010	0.005	2.026	0.044	0.000	0.019
MEDOCC	metals	0.047	0.081	0.580	0.562	-0.113	0.207
MEDOCC	sqrtmud	0.539	0.157	3.434	0.001	0.230	0.849
MEDOCC	sqrtTP	0.006	0.029	0.211	0.833	-0.051	0.063
MEDOCC	Type8	0.006	0.436	0.013	0.989	-0.854	0.866
MEDOCC	RegionEastern North Island	-0.673	0.192	-3.495	0.001	-1.052	-0.293
MEDOCC	RegionNorth Eastern	-0.025	0.151	-0.163	0.870	-0.323	0.273
MEDOCC	RegionSouthern	-0.133	0.151	-0.885	0.377	-0.430	0.163
MEDOCC	RegionWestern North Island	-0.203	0.389	-0.521	0.603	-0.969	0.564
MEDOCC	metals:sqrtmud	0.080	0.018	4.359	0.000	0.044	0.115
MEDOCC	metals:sqrtTP	-0.012	0.005	-2.547	0.012	-0.022	-0.003
MEDOCC	sqrtmud:sqrtTP	-0.014	0.007	-2.004	0.046	-0.028	0.000
MEDOCC	sqrtmud:Type8	-0.174	0.062	-2.808	0.005	-0.297	-0.052
MEDOCC	sqrtTP:Type8	0.048	0.027	1.768	0.078	-0.006	0.102
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.541	0.000	-0.022	-0.009
TBI	sqrtmud	-0.007	0.010	-0.662	0.509	-0.027	0.013
TBI	sqrtTP	0.010	0.003	3.761	0.000	0.005	0.015
TBI	Type8	-0.025	0.013	-1.915	0.057	-0.051	0.001
TBI	RegionEastern North Island	-0.109	0.019	-5.669	0.000	-0.146	-0.071
TBI	RegionNorth Eastern	0.076	0.016	4.820	0.000	0.045	0.107
TBI	RegionSouthern	-0.078	0.016	-4.860	0.000	-0.110	-0.046
TBI	RegionWestern North Island	0.004	0.041	0.108	0.914	-0.076	0.085
TBI	sqrtmud:sqrtTP	-0.001	0.000	-1.454	0.147	-0.002	0.000
TBI	metals:Type8	0.010	0.005	2.292	0.023	0.001	0.019

##Model summary multiple stressor interactions####
index_Summ5  <-
  glance(EG_models5, fit_index5)
kable(index_Summ5)

index	r.squared	adj.r.squared	sigma	statistic	p.value	df	logLik	AIC	BIC	deviance	df.residual
AMBI	0.346	0.308	0.853	9.07	0	14	-180.8	392	444	162.14	223
AMBI_S	0.470	0.439	0.583	15.20	0	14	-90.6	211	263	75.75	223
BENTIX	0.301	0.261	1.337	7.40	0	14	-287.4	605	657	398.61	223
BQI	0.557	0.540	1.716	31.73	0	10	-348.6	719	757	668.37	227
ITI	0.295	0.271	17.146	11.95	0	9	-894.7	1809	1844	67026.39	228
logN	0.578	0.558	1.058	28.04	0	12	-232.9	492	537	251.74	225
M_AMBI	0.431	0.408	0.129	19.09	0	10	264.9	-508	-470	3.77	227
MEDOCC	0.335	0.296	1.112	8.64	0	14	-243.6	517	569	275.52	223
S	0.553	0.533	4.875	27.94	0	11	-595.6	1215	1257	5371.88	226
TBI	0.582	0.564	0.121	31.52	0	11	280.4	-537	-495	3.31	226

ggplot(index_coef_all5,
       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_residuals5  <-
  augment(EG_models5, fit_index5)

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

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