library(tidyverse)
## Warning: package 'tidyverse' was built under R version 3.3.3
## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr
## Warning: package 'ggplot2' was built under R version 3.3.3
## Warning: package 'tibble' was built under R version 3.3.3
## Warning: package 'tidyr' was built under R version 3.3.3
## Warning: package 'readr' was built under R version 3.3.3
## Warning: package 'purrr' was built under R version 3.3.3
## Warning: package 'dplyr' was built under R version 3.3.3
## Conflicts with tidy packages ----------------------------------------------
## filter(): dplyr, stats
## lag(): dplyr, stats
library(modelr)
## Warning: package 'modelr' was built under R version 3.3.3
library(viridis)
library(lme4)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following object is masked from 'package:tidyr':
##
## expand
library(forcats)
library(sjPlot)
## Warning: package 'sjPlot' was built under R version 3.3.3
## Visit http://strengejacke.de/sjPlot for package-vignettes.
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)
###read data###
dat_mult <-
read_csv('dat_mult.csv') %>%
mutate(Type = factor(Type))
## Parsed with column specification:
## cols(
## .default = col_double(),
## Region = col_character(),
## n = col_integer(),
## council = col_character(),
## estuary = col_character(),
## site = col_character(),
## cesy = col_character(),
## ces = col_character(),
## cesym = col_character(),
## year = col_integer(),
## Type = col_integer()
## )
## See spec(...) for full column specifications.
## Warning: package 'bindrcpp' was built under R version 3.3.3
#####Grup data by index#####
by_index <-
dat_mult %>%
gather(index, value, AMBI:S) %>%
group_by(index) %>%
mutate(value = scale(value)) %>%
group_by(index) %>%
nest()
###model formula####
index_model <-
function(df) {
lmer(
value ~
Mud +
Metals +
TP +
(1 | Type) +
(0 + Mud|Type) +
(0 + TP|Type)+
(0 + Metals|Type) +
(1 |year) +
(1 |Region/estuary),
data =df,
REML = FALSE,
weights = sqrt(n),
na.action = "na.fail")
}
models_by_index <-
by_index %>%
mutate(model = map(data, index_model))
summary <-
models_by_index %>%
mutate(summary = map(model, broom::glance)) %>%
unnest(summary, .drop = TRUE)
res <-
models_by_index %>%
mutate(residuals = map(model, broom::augment)) %>%
unnest(residuals, .drop = TRUE)
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
## Warning: Deprecated: please use `purrr::possibly()` instead
res %>%
gather(stressor, value1, Mud:TP) %>%
group_by(stressor) %>%
ggplot(., aes(y = value, x = value1, color = Type)) +
geom_smooth(aes(y = .fitted), method = 'lm', size = 0.8) +
geom_point(alpha = 0.3) +
theme_javier() +
facet_grid(index~fct_relevel(stressor,c('Mud')))+
labs(x = '', y = '')
res %>%
gather(stressor, value1, Mud:TP) %>%
group_by(stressor) %>%
ggplot(., aes(y = value, x = value1, color = Region)) +
geom_smooth(aes(y = .fitted), method = 'lm', size = 0.8) +
geom_point(alpha = 0.3) +
theme_javier() +
facet_grid(index~fct_relevel(stressor,c('Mud')))+
labs(x = '', y = '')+
theme(legend.position = 'top')
coef <-
models_by_index %>%
mutate(coef = map(model, ~ broom::tidy(.x, conf.int = T))) %>%
unnest(coef, .drop = TRUE) %>%
data.frame(.) %>%
mutate(
term = factor(term),
term = fct_recode (
term,
Intercept = '(Intercept)',
Type = 'sd_(Intercept).Type.3',
`Mud x Type` = 'sd_Mud.Type.2',
`Metals x Type` = 'sd_Metals.Type',
`TP x Type` = 'sd_TP.Type.1',
`Estuary (Region)` = 'sd_(Intercept).estuary.Region',
Year = 'sd_(Intercept).year',
Region = 'sd_(Intercept).Region',
Residual = 'sd_Observation.Residual'
)
) %>%
write_csv(.,'lmm_coefficients.csv')
###Coefficient plots#####
ggplot(coef,
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) +
geom_vline(xintercept = 0,
lty = 2,
col = 'gray60') +
ylab('') +
xlab('Coefficients') +
theme_javier()+
scale_color_discrete(guide = F)
## Warning: Removed 72 rows containing missing values (geom_errorbarh).
####residuals plots#####
ggplot(res) +
geom_point(aes(x = .fitted, y = .resid), alpha = .3) +
facet_wrap(~ index, scale = 'free') +
geom_hline(yintercept = 0,
lty = 2,
col = 2)
ggplot(res) +
stat_qq(aes(sample = .wtres), alpha = .3) +
facet_wrap(~ index, scale = 'free') +
geom_abline(
intercept = 0,
slope = 1,
lty = 2,
col = 2
)
###Pseudo R-squared
library(MuMIn)
## Warning: package 'MuMIn' was built under R version 3.3.3
pseudo_rsq <-
models_by_index %>%
mutate(r_sq = map(model, MuMIn::r.squaredGLMM)) %>%
unnest(r_sq, .drop = TRUE) %>%
mutate(Type = rep(c('Marginal', 'Conditional'), 9),
index = factor(index)) %>%
arrange(Type,r_sq) %>%
print()
## # A tibble: 18 x 3
## index r_sq Type
## <fctr> <dbl> <chr>
## 1 AMBI 0.21360 Conditional
## 2 ITI 0.22641 Conditional
## 3 BENTIX 0.25547 Conditional
## 4 NZ-AMBI 0.32355 Conditional
## 5 M-AMBI 0.37759 Conditional
## 6 S 0.42459 Conditional
## 7 BQI 0.45602 Conditional
## 8 TBI 0.45629 Conditional
## 9 logN 0.51810 Conditional
## 10 logN 0.00859 Marginal
## 11 ITI 0.04574 Marginal
## 12 S 0.07560 Marginal
## 13 BQI 0.08811 Marginal
## 14 BENTIX 0.09985 Marginal
## 15 TBI 0.11051 Marginal
## 16 AMBI 0.12117 Marginal
## 17 M-AMBI 0.15110 Marginal
## 18 NZ-AMBI 0.21946 Marginal
###Pseudo R-square plot######
ggplot(pseudo_rsq) +
geom_col(aes(
x = fct_reorder(index, r_sq),
y = r_sq,
fill = Type
), position = 'stack') +
labs(y = 'Pseudo R - square', x = '') +
coord_flip() +
theme(legend.position = c(.8, .1))
###Plot firtted values (redundant with pots facet_grid above)####
####Fitted values by region###
ggplot(res,aes(Mud, y = value, color=Type)) +
facet_wrap(~index) +
geom_smooth(aes(y=.fitted),method = 'lm', size=0.8) +
geom_point(alpha = 0.3) +
theme_javier()
ggplot(res,aes(Metals, y = value, color=Type)) +
facet_wrap(~index) +
geom_smooth(aes(y=.fitted),method = 'lm', size=0.8) +
geom_point(alpha = 0.3) +
theme_javier()
ggplot(res,aes(TP, y = value, color=Type)) +
facet_wrap(~index) +
geom_smooth(aes(y=.fitted),method = 'lm', size=0.8) +
geom_point(alpha = 0.3) +
theme_javier()
###Fitted values
ggplot(res,aes(Mud, y = value, color=Type)) +
facet_wrap(~index) +
geom_smooth(aes(y=.fitted),method = 'lm', size=0.8) +
geom_point(alpha = 0.3) +
theme_javier()
ggplot(res,aes(Metals, y = value, color=Type)) +
facet_wrap(~index) +
geom_smooth(aes(y=.fitted),method = 'lm', size=0.8) +
geom_point(alpha = 0.3) +
theme_javier()
ggplot(res,aes(TP, y = value, color=Type)) +
facet_wrap(~index) +
geom_smooth(aes(y=.fitted),method = 'lm', size=0.8) +
geom_point(alpha = 0.3) +
theme_javier()
# ###plot random effects###
# ran_ef <-
# models_by_index %>%
# mutate(ran_ef = map(model, ~sjp.lmer(.x, type = "re", show.ci = F, sort.est = "(Intercept)")),
# ran_ef1 = map(model,~sjp.lmer(.x, type = "ri.slope")))