Summary of stats issues
Leo_Ohyama
2022-12-06
library(piecewiseSEM)
library(nlme)
library(parallel)
library(tidyverse)
library(sdmTMB)
library(patchwork)
#### corHaversine - spatial correlation with haversine distance for PGLS
source("Functions/Haversine_custom.R")
#load the subsampled data
setwd("~/Google Drive/Project_Macro/")
path_df<-readRDS("Data/Cell_data/Path_Model_Cell_Data_500_DEC22.rds")
path_df<-path_df %>% drop_na(lat)Data plots of omnivory vs. average size
library(ggpointdensity)
subsampled<-ggplot(data = path_df) +
geom_pointdensity(aes(x = prop_om, y = avg_size)) +
labs(x = "Omnivory (%)", y = "Size (mm)",
title = "Plot of subsampled size vs. omnivory\nby hexbin, 127,000 points") +
scale_color_viridis_c()
averaged<-ggplot(data = analysis_df) +
geom_point(aes(x = prop_om, y = avg_size)) +
labs(x = "Omnivory (%)", y = "Size (mm)",
title = "Plot of averaged size vs.\nomnivory by hexbin")
subsampled + averaged + plot_layout()
Path model with averaged data:
m1<-nlme::gls(log(avg_size) ~ prop_om,
correlation = corHaversine(form=~lon+lat,mimic="corSpher"),
method = "ML",
data = analysis_df)
m2<-nlme::gls(prop_om ~ ATR + NPP_cv,
correlation = corHaversine(form=~lon+lat,mimic="corSpher"),
method = "ML",
data = analysis_df)
#create path model
psem1<-piecewiseSEM::psem(m1, m2)
#summary of model
summary(psem1)##
|
| | 0%
|
|=================================== | 50%
|
|======================================================================| 100%##
## Structural Equation Model of psem1
##
## Call:
## log(avg_size) ~ prop_om
## prop_om ~ ATR + NPP_cv
##
## AIC BIC
## 33.140 65.425
##
## ---
## Tests of directed separation:
##
## Independ.Claim Test.Type DF Crit.Value P.Value
## log(avg_size) ~ ATR + ... coef 267 2.2994 0.0223 *
## log(avg_size) ~ NPP_cv + ... coef 267 2.2840 0.0232 *
##
## Global goodness-of-fit:
##
## Fisher's C = 15.14 with P-value = 0.004 and on 4 degrees of freedom
##
## ---
## Coefficients:
##
## Response Predictor Estimate Std.Error DF Crit.Value P.Value Std.Estimate
## log(avg_size) prop_om 0.2366 0.0737 267 3.2102 0.0015 0.2273
## prop_om ATR 0.0054 0.0018 267 2.9851 0.0031 0.3726
## prop_om NPP_cv 0.0087 0.0175 267 0.4966 0.6199 0.0309
##
## **
## **
##
##
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05
##
## ---
## Individual R-squared:
##
## Response method R.squared
## avg_size none 0.05
## prop_om none 0.18Path model above assesses size ~ omnivory (prop_om) and omnivory ~ temperature range (ATR) + variation of NPP (NPP_cv). Low R-square for size. Model passes both tests for normality and homogeneity of variance
Spatial model with random fields:
Here we make a mesh using the sdmtmb package (https://github.com/pbs-assess/sdmTMB):
mesh<-make_mesh(analysis_df, xy_cols = c("lon", "lat"), n_knots = 75, type = "cutoff_search")## cutoff = 1.00 | knots = 377 | ↓
## cutoff = 100.00 | knots = 26 | ↑
## cutoff = 10.00 | knots = 152 | ↓
## cutoff = 31.62 | knots = 48 | ↑
## cutoff = 17.78 | knots = 77 | ↓
## cutoff = 23.71 | knots = 57 | ↑
## cutoff = 20.53 | knots = 63 | ↑
## cutoff = 19.11 | knots = 71 | ↑
## cutoff = 18.43 | knots = 78 | ↓
## cutoff = 18.77 | knots = 75 | ✔plot(mesh)
Now we can model size ~ omnivory using a gamma distribution (log link): Below we build the model where we assess size ~ omnivory, fit it, and then assess the model summary:
m1<-sdmTMB(data = analysis_df, mesh = mesh, avg_size~scale(prop_om),
family = Gamma(link = "log"), spatial = "on",
control = sdmTMBcontrol(newton_loops = 1))
summary(m1$sd_report, select = "fixed", p.value = TRUE)## Estimate Std. Error z value Pr(>|z^2|)
## b_j 0.005063344 0.03042445 0.1664235 8.678237e-01
## b_j 0.047729718 0.01239273 3.8514286 1.174308e-04
## ln_tau_O 2.330773428 0.77002560 3.0268779 2.470938e-03
## ln_kappa -2.050683342 0.49784471 -4.1191225 3.803179e-05
## ln_phi 4.113134250 0.09627197 42.7241121 0.000000e+00Above is the model output, the top two rows represent intercept and slope beta of omnivory. Size increases as omnivory increases.
Below we assess the residuals of the model and it passes:
s_gamma <- simulate(m1, nsim = 500)
pred_fixed <- m1$family$linkinv(predict(m1)$est_non_rf)
m1_resid <- DHARMa::createDHARMa(
simulatedResponse = s_gamma,
observedResponse = analysis_df$avg_size,
fittedPredictedResponse = pred_fixed
)
plot(m1_resid)
Fitted values:
plot(ggeffects::ggeffect(m1))## $prop_om
We can fit another model of omniovry ~ ATR and NPP, just like the path model example:
analysis_df2<-analysis_df %>% drop_na()
mesh2<-make_mesh(analysis_df2, xy_cols = c("lon", "lat"), n_knots = 75, type = "cutoff_search")## cutoff = 1.00 | knots = 377 | ↓
## cutoff = 100.00 | knots = 26 | ↑
## cutoff = 10.00 | knots = 152 | ↓
## cutoff = 31.62 | knots = 48 | ↑
## cutoff = 17.78 | knots = 77 | ↓
## cutoff = 23.71 | knots = 57 | ↑
## cutoff = 20.53 | knots = 63 | ↑
## cutoff = 19.11 | knots = 71 | ↑
## cutoff = 18.43 | knots = 78 | ↓
## cutoff = 18.77 | knots = 75 | ✔m2<-sdmTMB(data = analysis_df, mesh = mesh2, prop_om~scale(ATR) + scale(NPP_cv),
family = Gamma(link = "inverse"), spatial = "on")
summary(m2$sd_report, select = "fixed", p.value = TRUE)## Estimate Std. Error z value Pr(>|z^2|)
## b_j 1.38561589 0.06453069 21.472201 2.832796e-102
## b_j -0.15721846 0.03106040 -5.061700 4.155339e-07
## b_j -0.02866367 0.01912121 -1.499051 1.338604e-01
## ln_tau_O 3.49853560 0.58887951 5.941004 2.832814e-09
## ln_kappa -3.05415852 0.52736776 -5.791326 6.983295e-09
## ln_phi 3.38197329 0.09468734 35.717271 2.132517e-279Above we see significant effect from both predictors (row 2 and 3 b_j) with omnivory decreasing with increasing seasonality and more varying spatial NPP.
Below we assess the residuals of the model and we see minor deviations but nothing too bad:
s_gamma <- simulate(m2, nsim = 500)
pred_fixed <- m2$family$linkinv(predict(m2)$est_non_rf)
m2_resid <- DHARMa::createDHARMa(
simulatedResponse = s_gamma,
observedResponse = analysis_df2$prop_om,
fittedPredictedResponse = pred_fixed
)
plot(m2_resid)