library(tidyverse) # Easily Install and Load the 'Tidyverse'
library(GGally) # Extension to 'ggplot2'
library(broom) # Convert Statistical Objects into Tidy Tibbles
library(broom.mixed) # Tidying Methods for Mixed Models
library(performance) # Assessment of Regression Models Performance
library(insight) # Easy Access to Model Information for Various Model Objects
library(MuMIn) # Multi-Model Inference
library(modelbased) # Estimation of Model-Based Predictions, Contrasts and Means
library(ggeffects)
library(betareg) # Beta Regression
library(glmmTMB) # Generalized Linear Mixed Models using Template Model Builder
library(patchwork) # The Composer of Plots
library(DHARMa) # Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models
library(kableExtra) # Construct Complex Table with 'kable' and Pipe Syntax
library(showtext) # Using Fonts More Easily in R Graphs
showtext_opts(dpi = 300)
font_add(family = "Gill Sans MT", regular = "GIL_____.TTF", italic = "GILI____.TTF")
showtext_auto()Contexto
En mamíferos carnívoros las interacciones generalmente son estudiadas por medio de patrones de co-ocurrencia y cámaras trampa.
Si la co-ocurrencia es un buen proxy de interacciones, los patrones deberían responder de igual manera a los procesos que afectan teóricamente la interacción
Paquetes
Datos
En este caso voy a usar los datos de co-ocurrencia temporal de competencia de mamíferos carnívoros. Los datos provienen de una revisión sistemática de literatura que hice en mi doctorado.
tempC_db <- read_csv("data/tempC_db.csv")[,-1] %>%
select(Ov_coeff, D_competitor, D_family, S_competitor, S_family, mass_ratio, diet_dist, Lat_abs, p_distance, Samp_dur, Locality, Label
) %>% # Selecciono variables
drop_na() # Hay unas variables sin información por lo que las elimino
dim(tempC_db)[1] 726 12Variable de respuesta
La variable de respuesta es el coeficiente de superposición temporal (Ridout y Linkie 2009). Este coeficiente toma valores de 0 cuando las especies no superponen su actividad y 1 cuando tienen superposición total.
Gráficos de las variables
Aquí algunos gráficos para poder visualizar los datos
Código
# Ov_coef vs categoric variables
cloud_Ov_plot <- function(data, var, var_string){
plot <- ggplot(data, aes(x= {{var}}, y= Ov_coeff, fill= {{var}}))+
ggdist::stat_halfeye(
adjust = .5,
width = .6,
.width = 0,
justification = -.3,
point_colour = NA)+
geom_boxplot(
width = .25,
outlier.shape = NA )+
geom_point(aes(color= {{var}}),
size = 3,
alpha = .5,
shape= 16,
fill= "black",
position = position_jitter(
seed = 1, width = .1)) +
coord_cartesian(xlim = c(1.2, NA), clip = "off")+
labs(y= "Overlap coefficient",
x= var_string) +
scale_fill_viridis_d(option = "mako", begin = 0.1,
end = 0.8)+
scale_colour_viridis_d(option = "mako", begin = 0.1,
end = 0.8)+
theme_bw()+
theme(text=element_text(family = "Gill Sans MT",size = 10),
legend.position = "none",
plot.margin= margin(t=0.5,r=0.5, b=1, l=1,
unit = "cm"))
return(plot)
}
# Dotplot function
dotplot_fun <- function(data, dimension){
d_plot <- data %>%
mutate(index = seq(n())) %>%
select_if(is.numeric) %>%
pivot_longer(cols = !index,
names_to = "Variable",
values_to = "Value") %>%
ggplot(aes(x= Value, y= index, col= Variable))+
geom_point(size= 2, alpha=0.6)+
scale_color_viridis_d(option = "mako",
begin = 0.1,
end = 0.8)+
facet_wrap(~Variable, scales = "free")+
labs(title = paste(dimension, "Dotplot" ),
y= "Order of ovservation")+
theme_bw()+
theme(legend.position = "none",
text = element_text(size=10, family = "Gill Sans MT"))
}Warning: Using the `size` aesthetic with geom_segment was deprecated in ggplot2 3.4.0.
i Please use the `linewidth` aesthetic instead.
Correlación de variables numéricas
Es importante evaluar si existe correlación fuerte entre las variables explicativas. Si hay dos variables correlacionadas fuertemente en un mismo modelo será imposible distinguir el efecto. Además, los errores pueden inflarse lo que puede generar que variables importantes pasen como no “significativas” (Dormann et al. 2013). En esta caso aplicaré una prueba de correlación de spearman.
Spat_cor <- tempC_db %>%
select_if(is.numeric) %>%
ggpairs(.,
upper = list(continuous= wrap("cor", method= "spearman",
digits=2, corSize= 80)),
lower = list( continuous= "smooth")) +
theme_bw()+
theme(text = element_text(size=8, family="Gill Sans MT"))
Spat_cor
Parece no existir una correlación preocupante
Scalar las variables
El modelo
Debido a que la variable de respuesta es una proporción que solo puede tomar valores que van de 0 a 1, usaremos un modelo lineal generalizado con familia de error Beta (Ferrari y Cribari-Neto 2004). En general no se recomienda modelar proporciones con modelo lineal general, ni realizar transformaciones de arcsin o logit (Douma y Weedon 2019).
La parametrización más común del glm beta es la de media-precisión, en donde \(\mu\) es la media esperada y \(\phi\) la precisión ( inverso de la dispersión) y se relacionan así:
\[ \frac{\mu(1- \mu)}{1+ \phi} \]
Los valores de estos dos parámetros hacen que pueda tomar muchas formas la curva de distribución, desde similar a gausiana hasta binomial. Lo que hace la beta muy flexible.
De forma algebraica la beta entenderse así:
\[ y \sim Beta (\mu, \phi) \]
\[ logit(\mu) = log (\frac{\mu}{1-\mu}) \]
Otra ventaja de la Beta es que el parámetro \(\phi\) puede ser estimado para todas las observaciones o modelado en función de variables usando un enlace log.
Estructura de efectos aleatorios
Mis datos provienen de artículos y hay localidades dentro de esto artículos. Puede que las estimaciones de co-ocurrencia sean similares entre artículos y localidades (o no) por lo que conviene checarlo. Siguiendo el protocolo de (Zuur et al. 2009), podemos comparar modelos con efectos aleatorios de otro que no.
Para ello vamos a usar el paquete glmmTMB (Brooks et al. 2017).
#Modelo sin efectos elatorios
Temp_r0 <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2,
data=tempC_data,
family=beta_family(), REML = T)
# Label aleatorio
Temp_r1 <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2+ (1|Label),
data=tempC_data,
family=beta_family(), REML = T)
# Label de cada locality
Temp_r2 <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+p_distance + diet_dist+Lat_abs+Samp_dur)^2+ (1|Label/Locality),
data=tempC_data,
family=beta_family(), REML = T)Comparando por AIC
temp_AICtab <- AICcmodavg::aictab(cand.set = list(Temp_r0, Temp_r1, Temp_r2),
sort = T,
modnames = c("no random",
"Label random",
"Label/Locality random"),
second.ord = F)
kableExtra::kable(temp_AICtab, digits = 2)| Modnames | K | AIC | Delta_AIC | ModelLik | AICWt | LL | Cum.Wt | |
|---|---|---|---|---|---|---|---|---|
| 2 | Label random | 24 | -386.96 | 0.00 | 1.00 | 0.73 | 217.48 | 0.73 |
| 3 | Label/Locality random | 25 | -384.96 | 2.00 | 0.37 | 0.27 | 217.48 | 1.00 |
| 1 | no random | 23 | -309.65 | 77.32 | 0.00 | 0.00 | 177.82 | 1.00 |
El modelo con efectos aleatorios del label es el mejor
Family: beta ( logit )
Formula:
Ov_coeff ~ (mass_ratio + I(mass_ratio^2) + p_distance + diet_dist +
Lat_abs + Samp_dur)^2 + (1 | Label)
Data: tempC_data
AIC BIC logLik deviance df.resid
-387.0 -276.9 217.5 -435.0 724
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
Label (Intercept) 0.1213 0.3484
Number of obs: 726, groups: Label, 105
Dispersion parameter for beta family (): 6.63
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.7407264 0.0967994 7.652 1.98e-14 ***
mass_ratio -0.0980909 0.0654577 -1.499 0.1340
I(mass_ratio^2) -0.0275596 0.0502955 -0.548 0.5837
p_distance -0.1290782 0.0908237 -1.421 0.1553
diet_distSame_diet -0.0338361 0.0988417 -0.342 0.7321
Lat_abs 0.0199587 0.0683350 0.292 0.7702
Samp_dur -0.0236972 0.0816522 -0.290 0.7716
mass_ratio:I(mass_ratio^2) -0.0053939 0.0179695 -0.300 0.7640
mass_ratio:p_distance -0.0192951 0.0469033 -0.411 0.6808
mass_ratio:diet_distSame_diet 0.1069080 0.0793013 1.348 0.1776
mass_ratio:Lat_abs -0.0070675 0.0389941 -0.181 0.8562
mass_ratio:Samp_dur 0.0470314 0.0379456 1.239 0.2152
I(mass_ratio^2):p_distance 0.0701428 0.0364851 1.923 0.0545 .
I(mass_ratio^2):diet_distSame_diet -0.0076185 0.0506246 -0.150 0.8804
I(mass_ratio^2):Lat_abs 0.0073499 0.0253046 0.290 0.7715
I(mass_ratio^2):Samp_dur -0.0414817 0.0271315 -1.529 0.1263
p_distance:diet_distSame_diet -0.0004413 0.0990214 -0.004 0.9964
p_distance:Lat_abs 0.0192958 0.0422704 0.456 0.6480
p_distance:Samp_dur 0.0572241 0.0477358 1.199 0.2306
diet_distSame_diet:Lat_abs 0.0739355 0.0744860 0.993 0.3209
diet_distSame_diet:Samp_dur 0.1004522 0.0991656 1.013 0.3111
Lat_abs:Samp_dur 0.0205697 0.0491088 0.419 0.6753
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Loading required namespace: carAnalysis of Deviance Table (Type II Wald chisquare tests)
Response: Ov_coeff
Chisq Df Pr(>Chisq)
mass_ratio 2.0807 1 0.14917
I(mass_ratio^2) 0.9613 1 0.32687
p_distance 4.4785 1 0.03432 *
diet_dist 0.2095 1 0.64719
Lat_abs 2.1231 1 0.14509
Samp_dur 0.0136 1 0.90727
mass_ratio:I(mass_ratio^2) 0.0901 1 0.76405
mass_ratio:p_distance 0.1692 1 0.68079
mass_ratio:diet_dist 1.8174 1 0.17762
mass_ratio:Lat_abs 0.0328 1 0.85618
mass_ratio:Samp_dur 1.5362 1 0.21518
I(mass_ratio^2):p_distance 3.6960 1 0.05454 .
I(mass_ratio^2):diet_dist 0.0226 1 0.88038
I(mass_ratio^2):Lat_abs 0.0844 1 0.77147
I(mass_ratio^2):Samp_dur 2.3376 1 0.12629
p_distance:diet_dist 0.0000 1 0.99644
p_distance:Lat_abs 0.2084 1 0.64804
p_distance:Samp_dur 1.4370 1 0.23062
diet_dist:Lat_abs 0.9853 1 0.32090
diet_dist:Samp_dur 1.0261 1 0.31107
Lat_abs:Samp_dur 0.1754 1 0.67532
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Supuestos del modelo
Este no es un modelo de conteos, pero existen métodos para saber si la variación que se obtiene es mayor o menor a la esperada por el modelo. El paquete DHARMa (Hartig 2021) ofrece una opción para modelos mixtos. DHARMa usa una aproximación con residuales simulados a partir de un proceso parecido a un bootstrap y una función de densidad empírica. Si es algo complejo pero esta mejor explicado en la página del paquete (🔗 aquí)
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 1.0825, p-value = 0.2
alternative hypothesis: two.sidedEn palabras mortales, en este test vamos a observar si la variación observada (linea roja), difiere de la que puede llegar a tener el modelo ajustado muchas veces (histograma barritas grises). Si la linea cae dentro del histograma, sugiere que el modelo no tiene sub/sobre dispersión. En nuestro caso, el modelo parece no tener problemas graves de sobre dispersión.
Aun así, vamos a provechar las ventajas de Dharma para explorar si el modelo sigue teniendo algunos problemas
Cómo se observa el modelo presenta problemas de homogeneidad de varianza. Esto puede ser causada por la mala especificación o por una o más variables. DHARMa nos puede ayudar a explorar esto
$uniformity
$uniformity$details
catPred: Diff_diet
One-sample Kolmogorov-Smirnov test
data: dd[x, ]
D = 0.080564, p-value = 0.04179
alternative hypothesis: two-sided
------------------------------------------------------------
catPred: Same_diet
One-sample Kolmogorov-Smirnov test
data: dd[x, ]
D = 0.057234, p-value = 0.1211
alternative hypothesis: two-sided
$uniformity$p.value
[1] 0.04178527 0.12110882
$uniformity$p.value.cor
[1] 0.08357054 0.12110882
$homogeneity
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 0.002 0.9645
724 
Test for location of quantiles via qgam
data: simulationOutput
p-value = 0.1511
alternative hypothesis: both
Test for location of quantiles via qgam
data: simulationOutput
p-value < 2.2e-16
alternative hypothesis: both
Test for location of quantiles via qgam
data: simulationOutput
p-value = 0.01025
alternative hypothesis: both
Test for location of quantiles via qgam
data: simulationOutput
p-value < 2.2e-16
alternative hypothesis: bothHay tres variables que pueden estar generando estos problemas. Quizás sea que la relación sea cuadrática
Recuerden que podemos modelar \(\phi\) para capturar la variación:
Parece que la formula con dispersión es a lo mejor que puedo llegar
Estructura de efectos fijos
Ahora vamos a buscar la estructura del mínimo modelo adecuado. En modelos mixtos, cuando usamos AIC para escoger la estructura de las variables fijas, necesitamos que el modelo se ajuste por MLE.
fix_temp <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+p_distance+ diet_dist+Lat_abs+Samp_dur)^2+ (1|Label),
dispformula = ~p_distance+ Samp_dur+Lat_abs,
data=tempC_data,
family=beta_family(),
REML = F,
na.action = "na.fail")
T_sel<- dredge(fix_temp,
rank = "AIC",
fixed = c("disp(p_distance)","disp(Lat_abs)","disp(Samp_dur)"),
m.lim= c(NA,6))
T_sel %>%
mutate(across(is.numeric, round, 2)) %>%
DT::datatable()Según el criterio de \(\Delta\) AIC <2, tenemos tres posibles modelos igualmente plausibles
mod1 <- glmmTMB(formula = Ov_coeff ~ Lat_abs + mass_ratio + p_distance +
(1 | Label),
dispformula = ~p_distance+ Samp_dur+Lat_abs,
data=tempC_data,
family=beta_family(),
REML = T,
na.action = "na.fail")
mod2 <- glmmTMB(formula = Ov_coeff ~ Lat_abs +mass_ratio+ I(mass_ratio^2) +
p_distance+
(1 | Label),
dispformula = ~p_distance+ Samp_dur+Lat_abs,
data=tempC_data,
family=beta_family(),
REML = T,
na.action = "na.fail")
mod3 <- glmmTMB(formula = Ov_coeff ~ mass_ratio+I(mass_ratio^2) + p_distance + (mass_ratio+I(mass_ratio^2)):p_distance+
(1 | Label),
dispformula = ~p_distance+ Samp_dur+Lat_abs,
data=tempC_data,
family=beta_family(),
REML = T,
na.action = "na.fail")
best_mods <- list(mod1, mod2, mod3) %>%
map(tidy , conf.int = T) %>%
reduce(rbind)
kbl(best_mods, caption = "Models", digits = 3) %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = F) %>%
pack_rows("Ov_coeff ~ Lat_abs +mass_ratio+ I(mass_ratio^2) +
p_distance",1,5) %>%
pack_rows("Lat_abs +mass_ratio+ I(mass_ratio^2) +
p_distance",6,11) %>%
pack_rows("Ov_coeff ~ mass_ratio+I(mass_ratio^2) + p_distance + (mass_ratio+I(mass_ratio^2)):p_distance",12,18)| effect | component | group | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|---|
| Ov_coeff ~ Lat_abs +mass_ratio+ I(mass_ratio^2) + p_distance | |||||||||
| fixed | cond | NA | (Intercept) | 0.681 | 0.051 | 13.349 | 0.000 | 0.581 | 0.780 |
| fixed | cond | NA | Lat_abs | 0.127 | 0.046 | 2.763 | 0.006 | 0.037 | 0.217 |
| fixed | cond | NA | mass_ratio | -0.067 | 0.030 | -2.238 | 0.025 | -0.126 | -0.008 |
| fixed | cond | NA | p_distance | -0.107 | 0.033 | -3.225 | 0.001 | -0.172 | -0.042 |
| ran_pars | cond | Label | sd__(Intercept) | 0.326 | NA | NA | NA | 0.249 | 0.428 |
| Lat_abs +mass_ratio+ I(mass_ratio^2) + p_distance | |||||||||
| fixed | cond | NA | (Intercept) | 0.701 | 0.056 | 12.572 | 0.000 | 0.592 | 0.810 |
| fixed | cond | NA | Lat_abs | 0.124 | 0.046 | 2.707 | 0.007 | 0.034 | 0.214 |
| fixed | cond | NA | mass_ratio | -0.051 | 0.035 | -1.442 | 0.149 | -0.120 | 0.018 |
| fixed | cond | NA | I(mass_ratio^2) | -0.020 | 0.023 | -0.897 | 0.370 | -0.065 | 0.024 |
| fixed | cond | NA | p_distance | -0.105 | 0.033 | -3.181 | 0.001 | -0.170 | -0.040 |
| ran_pars | cond | Label | sd__(Intercept) | 0.325 | NA | NA | NA | 0.247 | 0.427 |
| Ov_coeff ~ mass_ratio+I(mass_ratio^2) + p_distance + (mass_ratio+I(mass_ratio^2)):p_distance | |||||||||
| fixed | cond | NA | (Intercept) | 0.747 | 0.057 | 13.198 | 0.000 | 0.636 | 0.858 |
| fixed | cond | NA | mass_ratio | -0.048 | 0.035 | -1.376 | 0.169 | -0.118 | 0.021 |
| fixed | cond | NA | I(mass_ratio^2) | -0.043 | 0.023 | -1.832 | 0.067 | -0.089 | 0.003 |
| fixed | cond | NA | p_distance | -0.167 | 0.041 | -4.089 | 0.000 | -0.248 | -0.087 |
| fixed | cond | NA | mass_ratio:p_distance | -0.059 | 0.039 | -1.518 | 0.129 | -0.136 | 0.017 |
| fixed | cond | NA | I(mass_ratio^2):p_distance | 0.089 | 0.030 | 3.008 | 0.003 | 0.031 | 0.147 |
| ran_pars | cond | Label | sd__(Intercept) | 0.354 | NA | NA | NA | 0.274 | 0.456 |
Veamos si las estimaciones de los coeficientes son muy diferentes
CI_TC <- list(mod1, mod2, mod3) %>%
map(tidy, conf.int= T) %>%
reduce(rbind) %>%
filter(term!= "sd__(Intercept)" & term != "(phi)") %>%
mutate(model= c(rep("mod1", 4), rep("mod2", 5), rep("mod3", 6)))
ggplot(CI_TC, aes(x= estimate, y= term,
xmin= conf.low, xmax= conf.high))+
geom_pointrange( aes(col= model), position = position_dodge2(width = 0.4))+
geom_vline(xintercept = 0, linetype = "dashed")+
theme_bw()
Modo prueba de hipótesis
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: Ov_coeff
Chisq Df Pr(>Chisq)
mass_ratio 2.2707 1 0.131839
I(mass_ratio^2) 1.3229 1 0.250063
p_distance 9.6593 1 0.001884 **
diet_dist 0.0008 1 0.977330
Lat_abs 6.7078 1 0.009599 **
Samp_dur 0.0238 1 0.877446
mass_ratio:I(mass_ratio^2) 0.1527 1 0.695977
mass_ratio:p_distance 0.3021 1 0.582545
mass_ratio:diet_dist 1.6184 1 0.203317
mass_ratio:Lat_abs 0.0734 1 0.786490
mass_ratio:Samp_dur 2.4684 1 0.116155
I(mass_ratio^2):p_distance 5.0342 1 0.024852 *
I(mass_ratio^2):diet_dist 0.0367 1 0.848074
I(mass_ratio^2):Lat_abs 0.0166 1 0.897490
I(mass_ratio^2):Samp_dur 2.3275 1 0.127105
p_distance:diet_dist 0.0014 1 0.970558
p_distance:Lat_abs 2.1034 1 0.146969
p_distance:Samp_dur 0.7399 1 0.389700
diet_dist:Lat_abs 3.4185 1 0.064469 .
diet_dist:Samp_dur 0.9874 1 0.320386
Lat_abs:Samp_dur 0.1596 1 0.689519
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1p_mod <- glmmTMB(formula = Ov_coeff ~ p_distance+ Lat_abs+ I(mass_ratio^2):p_distance+
(1 | Label),
dispformula = ~p_distance+ Samp_dur+Lat_abs,
data=tempC_data,
family=beta_family(),
REML = T,
na.action = "na.fail")
p_mod %>%
tidy( conf.int = T) %>%
kbl( caption = "Models", digits = 3)| effect | component | group | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|---|
| fixed | cond | NA | (Intercept) | 0.687 | 0.052 | 13.236 | 0.000 | 0.585 | 0.789 |
| fixed | cond | NA | p_distance | -0.140 | 0.038 | -3.696 | 0.000 | -0.215 | -0.066 |
| fixed | cond | NA | Lat_abs | 0.126 | 0.047 | 2.704 | 0.007 | 0.035 | 0.218 |
| fixed | cond | NA | p_distance:I(mass_ratio^2) | 0.038 | 0.026 | 1.427 | 0.154 | -0.014 | 0.089 |
| ran_pars | cond | Label | sd__(Intercept) | 0.338 | NA | NA | NA | 0.259 | 0.441 |
Con que modelo quedarme?
Pues creo que la opción que tengo es evaluar en cual modelo los residuales se comportan mejor, así que usaremos de nuevo DHARMa
Residuales modelo 1
Residuales modelo 2
Residuales modelo 3
Residuales modelo anova
Yo me quedo con el modelo tres
Gráficos de predicción
p_distvec <- c(max(tempC_data$p_distance), mean(tempC_data$p_distance),
min(tempC_data$p_distance))
TC_pred <- ggeffect(mod3, terms =c("mass_ratio[all]", "p_distance[p_distvec]") ) %>%
mutate( mas_real= (x* attr(scale(tempC_db$mass_ratio), "scaled:scale"))+ attr(scale(tempC_db$mass_ratio), "scaled:center"),
p_disreal= (as.numeric(as.character(group))* attr(scale(tempC_db$p_distance), "scaled:scale"))+ attr(scale(tempC_db$p_distance), "scaled:center")) %>%
mutate(across(p_disreal, round, 2))
DT::datatable(TC_pred)Ahora el gráfico
(pred_plot <- ggplot(TC_pred)+
geom_point(data=tempC_db, aes(x= mass_ratio,
y= Ov_coeff),
size= 1.5, alpha= 0.4, col= "gray")+
geom_ribbon(aes(x= mas_real, y= predicted,
ymin= conf.low, ymax= conf.high,
fill= as.factor(p_disreal),
group= as.factor(p_disreal)),
alpha=0.4)+
geom_line(aes(x= mas_real, y=predicted,
group= as.factor(p_disreal)),
linewidth= 0.8 )+
labs(x= "Log (Mass ratio)",
y= "Temporal Overlap",
group= "Phylogenetic distance",
fill= "Phylogenetic distance",
tag= "A")+
ylim(c(0,1))+
theme_bw()+
theme(text = element_text(size=13, family = "Gill Sans MT"),
legend.position = c(0.5, 0.10),
legend.direction="horizontal"))
Referencias
Brooks, Mollie E., Kasper Kristensen, Koen J. van, Arni Magnusson, Casper W. Berg, Anders Nielsen, Hans J. Skaug, Martin Maechler, y Benjamin M. Bolker. 2017. «glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling» 9. https://doi.org/10.32614/RJ-2017-066.
Dormann, Carsten F., Jane Elith, Sven Bacher, Carsten Buchmann, Gudrun Carl, Gabriel Carré, Jaime R. García Marquéz, et al. 2013. «Collinearity: A Review of Methods to Deal with It and a Simulation Study Evaluating Their Performance». Ecography 36 (1): 27-46. https://doi.org/10.1111/j.1600-0587.2012.07348.x.
Douma, Jacob C., y James T. Weedon. 2019. «Analysing Continuous Proportions in Ecology and Evolution: A Practical Introduction to Beta and Dirichlet Regression». Methods in Ecology and Evolution 10 (9): 1412-30. https://doi.org/10.1111/2041-210X.13234.
Ferrari, Silvia, y Francisco Cribari-Neto. 2004. «Beta Regression for Modelling Rates and Proportions». Journal of Applied Statistics 31 (7): 799-815. https://doi.org/10.1080/0266476042000214501.
Hartig, Florian. 2021. «DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models». https://CRAN.R-project.org/package=DHARMa.
Ridout, M. S., y M. Linkie. 2009. «Estimating Overlap of Daily Activity Patterns From Camera Trap Data». Journal of Agricultural, Biological, and Environmental Statistics 14 (3): 322-37. https://doi.org/10.1198/jabes.2009.08038.
Zuur, Alain F, Elena N Ieno, Neil J Walker, Anatoly A Saveliev, y Graham M Smith. 2009. Mixed effects models and extensions in ecology with R. 1.ª ed. Statistics for Biology y Health. Springer.