PEP725 explorations
Beni Stocker
12/2/2020
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.4 ✓ dplyr 1.0.2
## ✓ tidyr 1.1.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(broom)
library(rbeni)##
## Attaching package: 'rbeni'## The following object is masked _by_ '.GlobalEnv':
##
## add_alphalibrary(lme4) #lmer## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(MuMIn) #r.squaredGLMM
library(lmerTest) #p-values##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(effects) #plot effects## Loading required package: carData## lattice theme set by effectsTheme()
## See ?effectsTheme for details.library(sjPlot) #plot effects## #refugeeswelcomedo_eval <- TRUEAnalysis
Read the data complemented with drivers, obtained from Constantin Zohner (7.12.2020).
df <- data.table::fread("~/data/pep/processed/DataMeta_3_Drivers_20_11_10.csv") %>%
as_tibble() %>%
rename(lon = LON, lat = LAT, year = YEAR, off = DoY_off, on = DoY_out,
anom_off_zani = autumn_anomaly, anom_on_zani = spring_anomaly,
species = Species, id_site = PEP_ID, id_species_site = timeseries) %>%
## use the on-water-stressed version of A
mutate(cA_tot = `cA_tot-w`)Select which one to use for further analysis.
df_anl <- dfInterannual
Get correlation between on and off for withing species-sites => interannual temporal variation.
get_coef_cA_tot <- function(df){df %>% filter(term == "cA_tot") %>% pull(estimate)}
get_p_cA_tot <- function(df){df %>% filter(term == "cA_tot") %>% pull(p.value)}
get_coef_year <- function(df){df %>% filter(term == "year") %>% pull(estimate)}
get_p_year <- function(df){df %>% filter(term == "year") %>% pull(p.value)}
df_temporal <- df_anl %>%
## for each species and site
ungroup() %>%
group_by(id_species_site) %>%
nest() %>%
## EOS ~ Assim
mutate(linmod = purrr::map(data, ~lm(off ~ cA_tot, data = .))) %>%
mutate(summ = purrr::map(linmod, ~summary(.))) %>%
mutate(df_coef = purrr::map(linmod, ~tidy(.))) %>%
mutate(coef_cA_tot = purrr::map_dbl(df_coef, ~get_coef_cA_tot(.))) %>%
mutate(p_value_cA_tot = purrr::map_dbl(df_coef, ~get_p_cA_tot(.))) %>%
select(-linmod) %>%
mutate(rsq = purrr::map_dbl(summ, "r.squared"),
adj_rsq = purrr::map_dbl(summ, "adj.r.squared")) %>%
select(-summ) %>%
## SOS temporal trend
mutate(linmod = purrr::map(data, ~lm(on ~ year, data = .))) %>%
mutate(df_coef = purrr::map(linmod, ~tidy(.))) %>%
mutate(coef_year_on = purrr::map_dbl(df_coef, ~get_coef_year(.))) %>%
mutate(p_value_year_on = purrr::map_dbl(df_coef, ~get_p_year(.))) %>%
select(-linmod, -df_coef) %>%
## EOS temporal trend
mutate(linmod = purrr::map(data, ~lm(off ~ year, data = .))) %>%
mutate(df_coef = purrr::map(linmod, ~tidy(.))) %>%
mutate(coef_year_off = purrr::map_dbl(df_coef, ~get_coef_year(.))) %>%
mutate(p_value_year_off = purrr::map_dbl(df_coef, ~get_p_year(.))) %>%
select(-linmod, -df_coef)
save(df_temporal, file = "~/data/df_interannual_assim.RData")Plot distribution of slopes
load("~/data/df_interannual_assim.RData")
df_temporal %>%
ggplot(aes(coef_year_on, ..density..)) +
geom_histogram() +
xlim(-2.5, 2.5) +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_vline(xintercept = median(df_temporal$coef_year_on, na.rm = TRUE), color = "red", linetype = "dashed") +
labs(title = "Temporal trend of SOS", subtitle = paste(nrow(df_temporal), "sites and species"))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
df_temporal %>%
ggplot(aes(coef_year_off, ..density..)) +
geom_histogram() +
xlim(-2.5, 2.5) +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_vline(xintercept = median(df_temporal$coef_year_off, na.rm = TRUE), color = "red", linetype = "dashed") +
labs(title = "Temporal trend of EOS", subtitle = paste(nrow(df_temporal), "sites and species"))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 11 rows containing non-finite values (stat_bin).
df_temporal %>%
ggplot(aes(coef_cA_tot, ..density..)) +
geom_histogram() +
xlim(-0.25, 0.25) +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_vline(xintercept = median(df_temporal$coef_cA_tot, na.rm = TRUE), color = "red", linetype = "dashed") +
labs(title = "Slope of EOS ~ cAtot", subtitle = paste(nrow(df_temporal), "sites and species"))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 1 rows containing non-finite values (stat_bin).## Warning: Removed 2 rows containing missing values (geom_bar).
# out <- df_temporal %>%
# analyse_modobs2("coef_year_on", "coef_year_off")
#
# out$gg +
# labs(title = "Trend in EOS vs. trend in SOS")Interannual temporal variation with LMMs
# Both models are equivalent
Fit_IAV_cA_tot = lmer(off ~ cA_tot + (1|id_site) + (1|species) , data = df_anl, na.action = "na.exclude")## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00312966 (tol = 0.002, component 1)summary(Fit_IAV_cA_tot)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: off ~ cA_tot + (1 | id_site) + (1 | species)
## Data: df_anl
##
## REML criterion at convergence: 3185265
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -10.1942 -0.6397 0.0026 0.6544 4.7960
##
## Random effects:
## Groups Name Variance Std.Dev.
## id_site (Intercept) 69.03 8.309
## species (Intercept) 37.55 6.128
## Residual 86.43 9.297
## Number of obs: 434226, groups: id_site, 3855; species, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.661e+02 2.513e+00 5.096e+00 145.7 2.01e-10 ***
## cA_tot -5.735e-02 1.387e-04 4.342e+05 -413.6 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cA_tot -0.079
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00312966 (tol = 0.002, component 1)r.squaredGLMM(Fit_IAV_cA_tot)## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.## R2m R2c
## [1,] 0.2392869 0.6593589plot(allEffects(Fit_IAV_cA_tot))
Fit_IAV_cGSI = lmer(off ~ cGSI + (1|id_site) + (1|species) , data = df_anl, na.action = "na.exclude")
summary(Fit_IAV_cGSI)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: off ~ cGSI + (1 | id_site) + (1 | species)
## Data: df_anl
##
## REML criterion at convergence: 3231234
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -11.5983 -0.6249 0.0287 0.6661 4.9373
##
## Random effects:
## Groups Name Variance Std.Dev.
## id_site (Intercept) 120.14 10.961
## species (Intercept) 26.31 5.129
## Residual 95.70 9.783
## Number of obs: 434226, groups: id_site, 3855; species, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.820e+02 2.122e+00 5.276e+00 180.0 3.36e-11 ***
## cGSI -6.822e-01 2.042e-03 4.320e+05 -334.1 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cGSI -0.138r.squaredGLMM(Fit_IAV_cGSI)## R2m R2c
## [1,] 0.2094307 0.6875532plot(allEffects(Fit_IAV_cGSI))
# Analysis of residuals
ResG <- residuals(Fit_IAV_cA_tot)
FittedG <- fitted(Fit_IAV_cA_tot)
par(mfrow=c(2,2))
plot(ResG ~ FittedG, xlab="Fitted values", ylab="Residuals", main="Residuals vs. fitted")
abline(h=0) ## Homocedasticity
plot(df_anl$off ~ FittedG, xlab="Fitted", ylab="Observed", main = "xyplot")
abline(0, 1)
hist(ResG, main="Histogram of residuals", xlab="Residuals") ## Normality
boxplot(ResG,ylab = "Residuals")
Weird pattern. Why is this?
df_anl <- df_anl %>%
ungroup() %>%
mutate(res = residuals(Fit_IAV_cA_tot), fit = fitted(Fit_IAV_cA_tot))
## location of weird points
plot_map_simpl() +
geom_point(
data = df_anl %>%
filter(fit > 310 & res < -40),
aes(x = lon, y = lat),
color = "red") +
xlim(-10,30) + ylim(30, 60)## Loading required package: sp## ### Welcome to rworldmap ##### For a short introduction type : vignette('rworldmap')## Regions defined for each Polygons## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.## Warning: Removed 4624 row(s) containing missing values (geom_path).
## distribution of A
ggplot() +
geom_histogram(
data = df_anl %>%
filter(fit > 310 & res < -40),
aes(x = cA_tot, y = ..density..),
fill = "red",
alpha = 0.5) +
geom_histogram(
data = df_anl %>%
filter(!(fit > 310 & res < -40)),
aes(x = cA_tot, y = ..density..),
alpha = 0.5) +
labs(title = "Assimilation")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## distribution of EOS
ggplot() +
geom_histogram(
data = df_anl %>%
filter(fit > 310 & res < -40),
aes(x = off, y = ..density..),
fill = "red",
alpha = 0.5) +
geom_histogram(
data = df_anl %>%
filter(!(fit > 310 & res < -40)),
aes(x = off, y = ..density..),
alpha = 0.5) +
labs(title = "EOS")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## distribution of years
ggplot() +
geom_histogram(
data = df_anl %>%
filter(fit > 310 & res < -40),
aes(x = year, y = ..density..),
fill = "red",
alpha = 0.5) +
geom_histogram(
data = df_anl %>%
filter(!(fit > 310 & res < -40)),
aes(x = year, y = ..density..),
alpha = 0.5) +
labs(title = "Year of observation")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
These are only very few observations with weirdly low A.
They are clearly off the linear regression between A and EOS.
df_anl <- df_anl %>%
mutate(off_fit = FittedG)
df_anl %>%
analyse_modobs2("off_fit", "off")## Loading required package: LSD##
## Attaching package: 'LSD'## The following objects are masked from 'package:rbeni':
##
## colorpalette, complementarycolor, convertcolor, convertgrey,
## daltonize, disco, distinctcolors, heatpairs, heatscatter,
## heatscatterpoints## Loading required package: ggthemes## Loading required package: RColorBrewer## $df_metrics
## # A tibble: 10 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 rmse standard 9.26e+ 0
## 2 rsq standard 5.35e- 1
## 3 mae standard 7.25e+ 0
## 4 n standard 4.34e+ 5
## 5 slope standard 1.01e+ 0
## 6 mean_obs standard 2.82e+ 2
## 7 prmse standard 3.28e- 2
## 8 pmae standard 2.57e- 2
## 9 bias standard -4.91e-12
## 10 pbias standard 1.10e- 3
##
## $gg## `geom_smooth()` using formula 'y ~ x'
##
## $linmod
##
## Call:
## lm(formula = obs ~ mod, data = df)
##
## Coefficients:
## (Intercept) mod
## -2.192 1.008
##
##
## $results
## # A tibble: 1 x 6
## rsq rmse mae bias slope n
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.535 9.26 7.25 -4.91e-12 1.01 434226Identify them using k-means clustering.
df_anl <- df_anl %>%
# mutate(cluster = kmeans(select(., off, off_fit), centers = 2)$cluster) %>%
rowwise() %>%
mutate(cluster = ifelse( (abs(off - off_fit) > 43) & off_fit > 300, "bad", "good")) %>%
ungroup()
df_anl %>%
ggplot(aes(x = off_fit, y = off, color = cluster)) +
geom_point()
Remove these sites.
print(nrow(df_anl))## [1] 434226df_anl <- df_anl %>%
filter(cluster == "good")
print(nrow(df_anl))## [1] 434017Run LMM for interannual temporal variation again.
# Both models are equivalent
Fit_IAV_cA_tot = lmer(off ~ cA_tot + (1|id_site) + (1|species) , data = df_anl, na.action = "na.exclude")
summary(Fit_IAV_cA_tot)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: off ~ cA_tot + (1 | id_site) + (1 | species)
## Data: df_anl
##
## REML criterion at convergence: 3173886
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.9966 -0.6499 -0.0011 0.6592 4.9611
##
## Random effects:
## Groups Name Variance Std.Dev.
## id_site (Intercept) 71.47 8.454
## species (Intercept) 38.47 6.202
## Residual 84.45 9.189
## Number of obs: 434017, groups: id_site, 3855; species, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.704e+02 2.544e+00 5.012e+00 145.6 2.77e-10 ***
## cA_tot -6.033e-02 1.403e-04 4.340e+05 -430.0 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cA_tot -0.079r.squaredGLMM(Fit_IAV_cA_tot)## R2m R2c
## [1,] 0.2527448 0.6753698plot(allEffects(Fit_IAV_cA_tot))
Fit_IAV_cGSI = lmer(off ~ cGSI + (1|id_site) + (1|species) , data = df_anl, na.action = "na.exclude")## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0276771 (tol = 0.002, component 1)summary(Fit_IAV_cGSI)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: off ~ cGSI + (1 | id_site) + (1 | species)
## Data: df_anl
##
## REML criterion at convergence: 3214513
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.9010 -0.6382 0.0253 0.6749 4.9173
##
## Random effects:
## Groups Name Variance Std.Dev.
## id_site (Intercept) 136.80 11.696
## species (Intercept) 25.67 5.067
## Residual 92.28 9.606
## Number of obs: 434017, groups: id_site, 3855; species, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.944e+02 2.099e+00 5.391e+00 187.9 1.7e-11 ***
## cGSI -7.686e-01 2.120e-03 4.322e+05 -362.5 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cGSI -0.145
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0276771 (tol = 0.002, component 1)r.squaredGLMM(Fit_IAV_cGSI)## R2m R2c
## [1,] 0.2351649 0.7229505plot(allEffects(Fit_IAV_cGSI))
# Analysis of residuals
ResG <- residuals(Fit_IAV_cA_tot)
FittedG <- fitted(Fit_IAV_cA_tot)
par(mfrow=c(2,2))
plot(ResG ~ FittedG, xlab="Fitted values", ylab="Residuals", main="Residuals vs. fitted")
abline(h=0) ## Homocedasticity
plot(df_anl$off ~ FittedG, xlab="Fitted", ylab="Observed", main = "xyplot")
abline(0, 1)
hist(ResG, main="Histogram of residuals", xlab="Residuals") ## Normality
boxplot(ResG,ylab = "Residuals")
The residuals show a systematic relationship with years. As mixed effects models show, there is a long-term trend (relationship with ‘year’ as a fixed effect, see below).
df_anl %>%
ungroup() %>%
mutate(res = residuals(Fit_IAV_cA_tot)) %>%
ggplot(aes(x = year, y = res)) +
geom_point() +
geom_smooth()## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
We can re-create Fig. 1A in Zani et al. This shows correlation of points, pooling data from all species, and after “removing” the random factor ‘site’. This essentially removes the spatial variation! This uses code adopted from ModelAnalyses_1_Correlation&PathAnalyses_Figures1.R.
remove_randomfactor <- function(df){
# fit model
f1 <- lmer(off ~ cA_tot + (1|id_site), data = df)
# Predict the data, removing the grouping value
fix.pred1 <- predict(f1, re.form = NA)
# To remove the random effects term and view the original spread of the data
# with just the random noise added, add the residual error back
# onto the predicted values
y.adj1 <- fix.pred1 + resid(f1)
# Store results
df$off_norandom <- y.adj1
return(df)
}
df_anl <- df_anl %>%
group_by(species) %>%
nest() %>%
mutate(data = map(data, ~remove_randomfactor(.))) %>%
unnest(data)
out <- df_anl %>%
analyse_modobs2("cA_tot", "off_norandom", type = "heat")
out$gg +
coord_cartesian(ylim=c(200,350))## `geom_smooth()` using formula 'y ~ x'
Ok. This shows the same. The point is that any spatial variation is treated as a random effect and effectively removed by this method!
Now, how does this relationship look like after aggregating to year bins?
remove_randomfactor_agg <- function(df, nbins){
## aggregate
df <- df %>%
mutate(yearbin = cut(year, breaks = nbins, labels = FALSE)) %>%
group_by(id_site, yearbin) %>%
summarise(cA_tot = mean(cA_tot, na.rm = TRUE), off = mean(off, na.rm = TRUE)) %>%
ungroup()
# fit model
f1 <- lmer(off ~ cA_tot + (1|id_site), data = df)
# Predict the data, removing the grouping value
fix.pred1 <- predict(f1, re.form = NA)
# To remove the random effects term and view the original spread of the data
# with just the random noise added, add the residual error back
# onto the predicted values
y.adj1 <- fix.pred1 + resid(f1)
# Store results
df$off_norandom <- y.adj1
return(df)
}
df_anl_agg <- df_anl %>%
group_by(species) %>%
nest() %>%
mutate(data = map(data, ~remove_randomfactor_agg(., nbins = 5))) %>%
unnest(data)## `summarise()` regrouping output by 'id_site' (override with `.groups` argument)
## `summarise()` regrouping output by 'id_site' (override with `.groups` argument)
## `summarise()` regrouping output by 'id_site' (override with `.groups` argument)
## `summarise()` regrouping output by 'id_site' (override with `.groups` argument)
## `summarise()` regrouping output by 'id_site' (override with `.groups` argument)
## `summarise()` regrouping output by 'id_site' (override with `.groups` argument)out <- df_anl_agg %>%
analyse_modobs2("cA_tot", "off_norandom", type = "heat")
out$gg +
coord_cartesian(ylim=c(200,350))## `geom_smooth()` using formula 'y ~ x'
Interesting. The negative relationship persists.
Long-term
Separating the temporal trends from the relationship with A using LMMs. Without interactions between A and year.
# Test without interaction
Fit_LT_on = lmer(off ~ scale(on) + scale(year) + (1|id_site) + (1|species), data = df_anl, na.action = "na.exclude")
summary(Fit_LT_on) ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: off ~ scale(on) + scale(year) + (1 | id_site) + (1 | species)
## Data: df_anl
##
## REML criterion at convergence: 3320943
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -10.0660 -0.6250 0.0232 0.6465 6.6071
##
## Random effects:
## Groups Name Variance Std.Dev.
## id_site (Intercept) 43.00 6.558
## species (Intercept) 40.12 6.334
## Residual 119.38 10.926
## Number of obs: 434017, groups: id_site, 3855; species, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.846e+02 2.588e+00 4.993e+00 109.96 1.21e-09 ***
## scale(on) 1.647e+00 2.396e-02 4.334e+05 68.72 < 2e-16 ***
## scale(year) 1.161e+00 2.230e-02 4.317e+05 52.10 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(n)
## scale(on) 0.001
## scale(year) 0.000 0.270plot(allEffects(Fit_LT_on))## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
## predictors scale(on), scale(year) are one-column matrices that were converted to
## vectors
## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
## predictors scale(on), scale(year) are one-column matrices that were converted to
## vectors
plot_model(Fit_LT_on, type = "eff", terms = c("on","year[2020,1960,1890]"),title = "",axis.title = c("on","off"),legend.title = "year") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),legend.position = "bottom") 
r.squaredGLMM(Fit_LT_on)## R2m R2c
## [1,] 0.01417371 0.4188375# Interaction scale(cA_tot):scale(year) is significant and positive.
Fit_LT_cA_tot = lmer(off ~ scale(cA_tot) + scale(year) + (1|id_site) + (1|species), data = df_anl, na.action = "na.exclude")
summary(Fit_LT_cA_tot) ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: off ~ scale(cA_tot) + scale(year) + (1 | id_site) + (1 | species)
## Data: df_anl
##
## REML criterion at convergence: 3118708
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.3689 -0.6582 0.0023 0.6642 5.1657
##
## Random effects:
## Groups Name Variance Std.Dev.
## id_site (Intercept) 80.64 8.980
## species (Intercept) 22.82 4.777
## Residual 74.20 8.614
## Number of obs: 434017, groups: id_site, 3855; species, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.838e+02 1.956e+00 5.077e+00 145.1 2.21e-10 ***
## scale(cA_tot) -9.953e+00 1.923e-02 4.340e+05 -517.5 < 2e-16 ***
## scale(year) 4.489e+00 1.849e-02 4.336e+05 242.8 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sc(A_)
## scal(cA_tt) 0.001
## scale(year) 0.000 -0.392plot(allEffects(Fit_LT_cA_tot))## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, :
## the predictors scale(cA_tot), scale(year) are one-column matrices that were
## converted to vectors## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, :
## the predictors scale(cA_tot), scale(year) are one-column matrices that were
## converted to vectors
plot_model(Fit_LT_cA_tot, type = "eff", terms = c("cA_tot","year[2020,1960,1890]"),title = "",axis.title = c("cA_tot","off"),legend.title = "year") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),legend.position = "bottom") 
r.squaredGLMM(Fit_LT_cA_tot)## R2m R2c
## [1,] 0.3307218 0.720474Fit_LT_cGSI = lmer(off ~ scale(cGSI) + scale(year) + (1|id_site) + (1|species), data = df_anl, na.action = "na.exclude")
summary(Fit_LT_cGSI) ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: off ~ scale(cGSI) + scale(year) + (1 | id_site) + (1 | species)
## Data: df_anl
##
## REML criterion at convergence: 3199502
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -10.1911 -0.6376 0.0277 0.6760 5.1099
##
## Random effects:
## Groups Name Variance Std.Dev.
## id_site (Intercept) 148.02 12.166
## species (Intercept) 17.91 4.232
## Residual 89.05 9.437
## Number of obs: 434017, groups: id_site, 3855; species, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.832e+02 1.739e+00 5.127e+00 162.9 1.01e-10 ***
## scale(cGSI) -9.533e+00 2.460e-02 4.328e+05 -387.5 < 2e-16 ***
## scale(year) 2.364e+00 1.912e-02 4.329e+05 123.7 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) s(GSI)
## scale(cGSI) 0.002
## scale(year) -0.001 -0.219plot(allEffects(Fit_LT_cGSI))## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
## predictors scale(cGSI), scale(year) are one-column matrices that were converted
## to vectors## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
## predictors scale(cGSI), scale(year) are one-column matrices that were converted
## to vectors
plot_model(Fit_LT_cGSI, type = "eff", terms = c("cGSI","year[2020,1960,1890]"),title = "",axis.title = c("cGSI","off"),legend.title = "year") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),legend.position = "bottom") 
r.squaredGLMM(Fit_LT_cGSI)## R2m R2c
## [1,] 0.2597083 0.7414438Let’s try to separate annual anomalies and site-level mean (of years 1970-2000).
separate_anom <- function(df){
df_mean <- df %>%
filter(year %in% 1970:2000) %>%
summarise(mean_on = mean(on, na.rm = TRUE),
mean_off = mean(off, na.rm = TRUE),
mean_cA_tot = mean(cA_tot, na.rm = TRUE),
mean_cGSI = mean(cGSI, na.rm = TRUE))
df %>%
mutate(anom_on = on - df_mean$mean_on,
anom_off = off - df_mean$mean_off,
mean_on = df_mean$mean_on,
mean_off = df_mean$mean_off,
mean_cA_tot = df_mean$mean_cA_tot,
anom_cA_tot = cA_tot - df_mean$mean_cA_tot,
mean_cGSI = df_mean$mean_cGSI,
anom_cGSI = cGSI - df_mean$mean_cGSI
)
}
df_anl <- df_anl %>%
group_by(id_species_site) %>%
nest() %>%
mutate(data = map(data, ~separate_anom(.))) %>%
unnest(data)
Fit_LT_anom = lmer(off ~ scale(mean_cA_tot) + scale(anom_cA_tot) + scale(year) + (1|id_site) + (1|species), data = df_anl, na.action = "na.exclude")
summary(Fit_LT_anom)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: off ~ scale(mean_cA_tot) + scale(anom_cA_tot) + scale(year) +
## (1 | id_site) + (1 | species)
## Data: df_anl
##
## REML criterion at convergence: 3087573
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.5488 -0.6577 0.0022 0.6614 5.6943
##
## Random effects:
## Groups Name Variance Std.Dev.
## id_site (Intercept) 164.85 12.840
## species (Intercept) 34.08 5.838
## Residual 71.95 8.482
## Number of obs: 431164, groups: id_site, 3789; species, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.837e+02 2.392e+00 5.090e+00 118.6 5.87e-10 ***
## scale(mean_cA_tot) -1.193e+01 4.511e-02 3.574e+05 -264.4 < 2e-16 ***
## scale(anom_cA_tot) -7.010e+00 1.439e-02 4.275e+05 -487.1 < 2e-16 ***
## scale(year) 4.381e+00 1.831e-02 4.291e+05 239.3 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(m_A_) scl(n_A_)
## scl(mn_cA_) 0.000
## scl(nm_cA_) 0.000 0.083
## scale(year) 0.000 -0.061 -0.396plot(allEffects(Fit_LT_anom))## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
## predictors scale(mean_cA_tot), scale(anom_cA_tot), scale(year) are one-column
## matrices that were converted to vectors
## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
## predictors scale(mean_cA_tot), scale(anom_cA_tot), scale(year) are one-column
## matrices that were converted to vectors
## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
## predictors scale(mean_cA_tot), scale(anom_cA_tot), scale(year) are one-column
## matrices that were converted to vectors
r.squaredGLMM(Fit_LT_anom)## R2m R2c
## [1,] 0.379972 0.835317Fit_LT_anom = lmer(off ~ scale(mean_cGSI) + scale(anom_cGSI) + scale(year) + (1|id_site) + (1|species), data = df_anl, na.action = "na.exclude")
summary(Fit_LT_anom)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: off ~ scale(mean_cGSI) + scale(anom_cGSI) + scale(year) + (1 |
## id_site) + (1 | species)
## Data: df_anl
##
## REML criterion at convergence: 3172412
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -10.2332 -0.6368 0.0273 0.6726 5.0350
##
## Random effects:
## Groups Name Variance Std.Dev.
## id_site (Intercept) 323.37 17.982
## species (Intercept) 20.90 4.572
## Residual 87.23 9.340
## Number of obs: 431164, groups: id_site, 3789; species, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.827e+02 1.889e+00 5.247e+00 149.6 9.96e-11 ***
## scale(mean_cGSI) -1.360e+01 6.753e-02 2.849e+05 -201.4 < 2e-16 ***
## scale(anom_cGSI) -5.390e+00 1.489e-02 4.273e+05 -362.1 < 2e-16 ***
## scale(year) 2.337e+00 1.899e-02 4.283e+05 123.1 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(m_GSI) scl(n_GSI)
## scl(mn_GSI) 0.004
## scl(nm_GSI) 0.001 0.076
## scale(year) -0.001 -0.043 -0.221plot(allEffects(Fit_LT_anom))## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, :
## the predictors scale(mean_cGSI), scale(anom_cGSI), scale(year) are one-column
## matrices that were converted to vectors## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, :
## the predictors scale(mean_cGSI), scale(anom_cGSI), scale(year) are one-column
## matrices that were converted to vectors
## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, :
## the predictors scale(mean_cGSI), scale(anom_cGSI), scale(year) are one-column
## matrices that were converted to vectors
r.squaredGLMM(Fit_LT_anom)## R2m R2c
## [1,] 0.3196238 0.8624611Whatever I do, there is always this weird negative relationship between A and EOS at the spatial (between-site) scale.
Spatial by species
Aggregate years, to get a “spatial” data frame, distinguishing speies.
df_spatial_sps <- df_anl %>%
group_by(id_site, species) %>%
summarise(cA_tot = mean(cA_tot, na.rm = TRUE), off = mean(off, na.rm = TRUE), cGSI = mean(cGSI, na.rm = TRUE), on = mean(on, na.rm = TRUE))## `summarise()` regrouping output by 'id_site' (override with `.groups` argument)Fit_bySps_cA_tot = lmer(off ~ cA_tot + (1|species), data = df_spatial_sps, na.action = "na.exclude")
summary(Fit_bySps_cA_tot)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: off ~ cA_tot + (1 | species)
## Data: df_spatial_sps
##
## REML criterion at convergence: 100417
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.1321 -0.6104 0.0165 0.6331 4.7047
##
## Random effects:
## Groups Name Variance Std.Dev.
## species (Intercept) 44.58 6.677
## Residual 55.94 7.479
## Number of obs: 14626, groups: species, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.132e+02 2.909e+00 6.479e+00 107.66 8.77e-12 ***
## cA_tot -1.979e-02 7.055e-04 1.462e+04 -28.05 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cA_tot -0.349plot(allEffects(Fit_bySps_cA_tot))
r.squaredGLMM(Fit_bySps_cA_tot)## R2m R2c
## [1,] 0.03559164 0.4633029Fit_bySps_cGSI = lmer(off ~ cGSI + (1|species), data = df_spatial_sps, na.action = "na.exclude")
summary(Fit_bySps_cGSI)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: off ~ cGSI + (1 | species)
## Data: df_spatial_sps
##
## REML criterion at convergence: 101074.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4586 -0.6151 0.0197 0.6363 5.1038
##
## Random effects:
## Groups Name Variance Std.Dev.
## species (Intercept) 45.01 6.709
## Residual 58.53 7.650
## Number of obs: 14626, groups: species, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.745e+02 2.916e+00 6.411e+00 94.15 2.59e-11 ***
## cGSI 7.075e-02 6.901e-03 1.462e+04 10.25 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cGSI -0.342plot(allEffects(Fit_bySps_cGSI))
r.squaredGLMM(Fit_bySps_cGSI)## R2m R2c
## [1,] 0.004510397 0.4372901Fit_bySps_on = lmer(off ~ on + (1|species), data = df_spatial_sps, na.action = "na.exclude")
summary(Fit_bySps_on)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: off ~ on + (1 | species)
## Data: df_spatial_sps
##
## REML criterion at convergence: 101177.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.7741 -0.6114 0.0212 0.6466 5.0251
##
## Random effects:
## Groups Name Variance Std.Dev.
## species (Intercept) 43.33 6.582
## Residual 58.95 7.678
## Number of obs: 14626, groups: species, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.855e+02 2.886e+00 6.639e+00 98.932 9.04e-12 ***
## on -6.776e-03 9.316e-03 1.462e+04 -0.727 0.467
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## on -0.364plot(allEffects(Fit_bySps_on))
r.squaredGLMM(Fit_bySps_on)## R2m R2c
## [1,] 4.2272e-05 0.4236526There is still a negative relationship between EOS and A!
How are these drivers distributed in space?
Get correlation between on and off for withing species-sites
getmostfrequent <- function(vec){names(sort(table(vec))[1])}
df_spatial <- df_anl %>%
group_by(id_site) %>%
summarise(cA_tot = mean(cA_tot, na.rm = TRUE),
off = mean(off, na.rm = TRUE),
cGSI = mean(cGSI, na.rm = TRUE),
on = mean(on, na.rm = TRUE),
species = getmostfrequent(species)) %>%
left_join(select(df_anl, id_site, lon, lat) %>% unique(),
by = "id_site")## `summarise()` ungrouping output (override with `.groups` argument)## Adding missing grouping variables: `id_species_site`save(df_spatial, file = "data/df_spatial_assim.RData")Photosynthesis across space.
plot_map_simpl(lonmin = 0, lonmax = 30, latmin = 40, latmax = 60) +
geom_point(data = df_spatial, aes(x = lon, y = lat, color = off)) +
scale_color_viridis_c() +
labs(title = "EOS")## Regions defined for each Polygons## Warning: Removed 4793 row(s) containing missing values (geom_path).
plot_map_simpl(lonmin = 0, lonmax = 30, latmin = 40, latmax = 60) +
geom_point(data = df_spatial, aes(x = lon, y = lat, color = cGSI)) +
scale_color_viridis_c() +
labs(title = "GSI")## Regions defined for each Polygons## Warning: Removed 4793 row(s) containing missing values (geom_path).
plot_map_simpl(lonmin = 0, lonmax = 30, latmin = 40, latmax = 60) +
geom_point(data = df_spatial, aes(x = lon, y = lat, color = cA_tot)) +
scale_color_viridis_c() +
labs(title = "Assimilation")## Regions defined for each Polygons## Warning: Removed 4793 row(s) containing missing values (geom_path).
There are a few points with very low values for A (and probably late senescence) located in Switzerland - looks like in the Canton of Graubünden… (sorry for the ugly map)
plot_map_simpl(lonmin = 5, lonmax = 12, latmin = 42, latmax = 50) +
geom_point(data = df_spatial, aes(x = lon, y = lat, color = off)) +
scale_color_viridis_c() +
labs(title = "EOS")## Regions defined for each Polygons## Warning: Removed 5112 row(s) containing missing values (geom_path).## Warning: Removed 10134 rows containing missing values (geom_point).
plot_map_simpl(lonmin = 5, lonmax = 12, latmin = 42, latmax = 50) +
geom_point(data = df_spatial, aes(x = lon, y = lat, color = cGSI)) +
scale_color_viridis_c() +
labs(title = "cGSI")## Regions defined for each Polygons## Warning: Removed 5112 row(s) containing missing values (geom_path).
## Warning: Removed 10134 rows containing missing values (geom_point).
plot_map_simpl(lonmin = 5, lonmax = 12, latmin = 42, latmax = 50) +
geom_point(data = df_spatial, aes(x = lon, y = lat, color = cA_tot)) +
scale_color_viridis_c() +
labs(title = "Assimilation")## Regions defined for each Polygons## Warning: Removed 5112 row(s) containing missing values (geom_path).
## Warning: Removed 10134 rows containing missing values (geom_point).
plot_map_simpl(lonmin = 5, lonmax = 12, latmin = 42, latmax = 50) +
geom_point(data = df_spatial, aes(x = lon, y = lat, color = species)) +
labs(title = "Species (most frequent at site)")## Regions defined for each Polygons## Warning: Removed 5112 row(s) containing missing values (geom_path).
## Warning: Removed 10134 rows containing missing values (geom_point).
Yes. They are definitely in the mountains. And mostly Larix. Where else are Larix decidua distributed?
plot_map_simpl(lonmin = 0, lonmax = 30, latmin = 40, latmax = 60) +
geom_point(data = dplyr::filter(df_spatial, species == "Larix decidua"), aes(x = lon, y = lat)) +
labs(title = "Larix decidua")## Regions defined for each Polygons## Warning: Removed 4793 row(s) containing missing values (geom_path).
Ok. A bit everywhere.
Is there an interaction effect with elevation? See next section.
What’s up with these sites that have a late season end and low A?
df_spatial_sps <- df_spatial_sps %>%
ungroup() %>%
mutate(fit = fitted(Fit_bySps_cA_tot))
df_spatial_sps %>%
ggplot() +
geom_line(aes(cA_tot, fit, color = species)) +
geom_point(aes(cA_tot, off, color = species), alpha = 0.1)
Interaction effect with elevation?
df_elv <- data.table::fread("~/data/pep/processed/DataMeta_2_PhenologyObs_PEP725_CleanData.csv") %>%
as_tibble() %>%
select(id_site = s_id, alt, alt_dem) %>%
distinct() %>%
mutate(elv = ifelse(is.na(alt), alt_dem, alt)) %>%
select(-alt, -alt_dem)
df_spatial_sps <- df_spatial_sps %>%
left_join(df_elv, by = c("id_site"))
Fit_bySps_cA_tot = lmer(off ~ cA_tot + elv + (1|species), data = df_spatial_sps, na.action = "na.exclude")
summary(Fit_bySps_cA_tot)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: off ~ cA_tot + elv + (1 | species)
## Data: df_spatial_sps
##
## REML criterion at convergence: 97866.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.4190 -0.6069 0.0243 0.6373 4.8601
##
## Random effects:
## Groups Name Variance Std.Dev.
## species (Intercept) 46.51 6.820
## Residual 53.54 7.317
## Number of obs: 14344, groups: species, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.192e+02 2.993e+00 6.665e+00 106.66 5.03e-12 ***
## cA_tot -2.325e-02 7.580e-04 1.434e+04 -30.67 < 2e-16 ***
## elv -3.920e-03 2.894e-04 1.434e+04 -13.55 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cA_tot
## cA_tot -0.365
## elv -0.035 0.025plot(allEffects(Fit_bySps_cA_tot))
r.squaredGLMM(Fit_bySps_cA_tot)## R2m R2c
## [1,] 0.04922484 0.4912277Fit_bySps_cGSI = lmer(off ~ cGSI + elv + (1|species), data = df_spatial_sps, na.action = "na.exclude")
summary(Fit_bySps_cGSI)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: off ~ cGSI + elv + (1 | species)
## Data: df_spatial_sps
##
## REML criterion at convergence: 98692.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.5369 -0.6144 0.0283 0.6461 5.1453
##
## Random effects:
## Groups Name Variance Std.Dev.
## species (Intercept) 46.23 6.799
## Residual 56.73 7.532
## Number of obs: 14344, groups: species, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.733e+02 3.099e+00 7.753e+00 88.210 6.40e-13 ***
## cGSI 8.220e-02 9.133e-03 1.434e+04 9.001 < 2e-16 ***
## elv -2.022e-03 3.515e-04 1.434e+04 -5.752 8.98e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cGSI
## cGSI -0.443
## elv -0.258 0.531plot(allEffects(Fit_bySps_cGSI))
r.squaredGLMM(Fit_bySps_cGSI)## R2m R2c
## [1,] 0.009618693 0.4543126Fit_bySps_on = lmer(off ~ on + elv + (1|species), data = df_spatial_sps, na.action = "na.exclude")
summary(Fit_bySps_on)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: off ~ on + elv + (1 | species)
## Data: df_spatial_sps
##
## REML criterion at convergence: 98751.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.0458 -0.6068 0.0311 0.6471 5.1232
##
## Random effects:
## Groups Name Variance Std.Dev.
## species (Intercept) 46.45 6.815
## Residual 56.97 7.548
## Number of obs: 14344, groups: species, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.806e+02 2.995e+00 6.704e+00 93.674 1.06e-11 ***
## on 4.668e-02 1.009e-02 1.433e+04 4.626 3.76e-06 ***
## elv -4.253e-03 3.213e-04 1.434e+04 -13.237 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) on
## on -0.369
## elv 0.111 -0.370plot(allEffects(Fit_bySps_on))
r.squaredGLMM(Fit_bySps_on)## R2m R2c
## [1,] 0.007906892 0.4535033