Presentación

En este documento encontraremos una serie de procedimientos analíticos realizados para entender los conceptos teóricos vistos en clase acerca de los análisis de datos desde la perspectiva Bayesiana.

Caso de estudio

Tenemos un set de datos que consta de temperaturas corporales, ambientales y de sustrato tomadas para un estudio de la ecología térmica de una especie de serpiente (Atractus marthae) semifosorial altoandina (datos propios).

Preguntas de estudio

1- ¿La serpiente es termorreguladora o termoconformista?

2- ¿La serpiente usa una estrategia heliotérmica o tigmotérmica?

Razonamiento

de los datos

En este caso de estudio tenemos un único tipo de datos que son temperaturas, las cuales, están dadas en grados Celsius para las 3 variables tomadas en campo (temperatura corporal, temperatura ambiental y temperatura del sustrato). Todos los valores de las temperaturas están dados en número reales positivos.

En la siguiente tabla tenemos los datos tomados en campo en donde adicionalmente se observan otros datos asociados que para efectos de este ejercicio vamos a ignorar.

atractus=read.csv("temp_atractus.csv", header=T, sep = ";")
atractus

del modelo

  • Queremos entonces evaluar una relación matemática entre una variable aleatoria de respuesta (temperatura corporal) y dos variables explicativas numéricas (temperaturas ambientales y de sustrato).

  • En este caso la variable de respuesta posee valores continuos por lo tanto entendemos que tiene asociada una función de densidad de probabilidades, ejemplo: Normal, que está definida por dos parámetros (media = μ, varianza = σ2).

  • Entendemos entonces, que el modelo estadístico que aplica en este caso es el lineal general (GLM), el cual, postula que las variables explicativas determinan o afectan a estos parámetros de la función de densidad de probabilidades asociada a la variable de repuesta.

aspectos clave del GLM

  • La variable de respuesta tiene distribución Normal (debemos evaluarlo).

  • Hay una relación lineal entre los parámetros y las variables explicativas.

  • Las variables explicativas pueden ser numéricas y/o categóricas.

Implementación

Entendiendo ya previamente los tipos de variables que vamos a manejar, las preguntas que queremos resolver y el modelo, en este caso lineal, que aplicaría para responder a estas preguntas, definimos la forma de proceder a continuación:

1- Definir la distribución de probabilidad de la variable respuesta

2- Examinar la relación que existe entre la variable respuesta y las variables explicativas

3- Formular el modelo, lo cual incluye definir las distribuciones previas

4- Ajustar el modelo

5- Evaluar la calidad del ajuste

1-) Distribución de la variable respuesta (t_cor):

Con lo cual se define la Verosimilitud del modelo, de este modo, para evaluarlo utilicé la librería (fitdistrplus)

library(fitdistrplus)

Vamos a considerar 3 distribuciones candidatas, para las cuales, ajustaremos cada distribución estimando sus dos parámetros, y con esto, calculamos la distribución acumulada teórica esperada si los datos tuvieran la distribución lognormal, normal o gamma.

lognor.t_cor=fitdist(atractus$t_cor,"lnorm")
normal.t_cor=fitdist(atractus$t_cor,"norm")
gamma.t_cor=fitdist(atractus$t_cor,"gamma")

luego de guardar estos valores estimados los graficamos para interpretarlos de mejor manera y elegir la distribución que esté más cerca de la distribución de los datos de nuestra variable respuesta (t_cor). Para esto usamos dos gráficos complementarios usando las funciones cdfcomp y qqcomp de la librería (fitdistrplus) de la siguiente manera:

cdf.cor2=cdfcomp(list(lognor.t_cor,normal.t_cor,gamma.t_cor),xlogscale = T,plotstyle = "ggplot") 
qq.cor2=qqcomp(list(lognor.t_cor,normal.t_cor,gamma.t_cor),
               plotstyle = "ggplot")

Realizamos el plot de estos dos gráficos usando la función grid.arrange de la librería (gridExtra) así:

library(gridExtra)
grid.arrange(cdf.cor2,qq.cor2,ncol=2)

Estos dos gráficos nos muestran las 3 distribuciones candidatas ajustadas a la distribución empírica de los datos.

El gráfico de la izquierda nos muestra las distribuciones acumuladas, donde, los puntos negros corresponden a la distribución empírica de los datos y las líneas y puntos de colores corresponden a cada una de las distribuciones candidatas como se indíca en la figura.

El gráfico de la derecha es un QQplot (gráfico de cuantiles), entre los cuantiles teóricos (eje horizontal) si la variable tuviera distribución normal, lognormal o gamma, contra, los cuantiles empíricos (eje vertical). Con lo cual, debemos elegir una de las distribuciones que esté más cerca de la línea recta que indica el ajuste perfecto.

Elegimos la distribución

luego de analizar el ajuste de las 3 distribuciones candidatas, para este ejercicio escogemos la distribución lognormal.

2-) Relación entre las variables (t_cor, t_am, t_sus):

Para esto, realicé un gráfico de correlación entre las variables con el fin de ver que tan correlacionadas están las variables explicativas. En este caso usamos la función ggpairs de la librería (GGally) de la siguiente manera:

library("GGally")
temp_1=ggpairs(atractus, mapping = aes (color = temp)) + theme_bw() + theme(legend.text = element_text(size = 14),legend.title = element_text(size = 16, face = "bold"))
temp_1

En este gráfico tenemos para efectos de simple visualización, a todas las variables iniciales tomadas durante el trabajo. Sin embargo, solo vamos a analizar las variables que hemos discutido anteriormente. t_cor, t_sus, t_am.

En la diagonal superior vemos los valores de correlación entre las variables explicativas, con los cuales debemos revisar y decidir, en este caso, usando un criterio heurístico definido por Dorman et al. 2013, el cual indica que podemos descartar a la variables que tiene una correlación en valor absoluto mayor a 0,7 por estar fuertemente correlacionadas y serían redundantes a la hora de explicar la variación de la variable de respuesta.

Por lo tanto, vemos que nuestras variables explicativas están fuertemente correlacionadas y son redundantes con respecto a la variación de la variable de respuesta, con lo cual, podríamos descartar el uso de por lo menos la variable menos correlacionada para efectos de este ejercicio.

3-) Modelo estadístico y distribución previa:

La variable de respuesta tiene distribución lognormal como definimos previamente, entonces, transformamos esta variable con un logaritmo y obtenemos que un logaritmo de una variable lognormal se convierte en Normal, de la siguiente manera:

t_cor ~ Lognormal = logt_cor ~ Normal( μ, σ)

atractus$logt_cor=log(atractus$t_cor)

con esto entonces ya podemos modelar la media de la variable respuesta t_cor(μ) como función lineal de las variables explicativas t_am y t_sus.

Esto tiene entonces, 2 pendientes parciales para la t_cor, más el intercepto, más la variabilidad al rededor del plano ajustado denotado por la desviación estándar (σ). Con lo cual, debemos especificar 4 distribuciones previas, una para cada parámetro.

t_cor(μ) = β0 + β1X1 + β2X2 : 3 parámetros + t_cor (σ)

Estandarizar las variables explicativas

Es útil para explicar los resultados y para ajustar el modelo. Acá entonces, estandarizamos a media = 0, sd = 1, y, creamos un nuevo dataframe junto con nuestra variable t_cor transformada anteriormente.

De la siguiente manera:

atra.s=scale(atractus[,c("lrc","t_sus","t_am")],center = T,scale = T) 
atra.s=as.data.frame(cbind(logt_cor=atractus$logt_cor, atra.s))
summary (atra.s)
##     logt_cor          lrc              t_sus              t_am        
##  Min.   :2.351   Min.   :-2.9211   Min.   :-1.7607   Min.   :-1.3698  
##  1st Qu.:2.614   1st Qu.:-0.1827   1st Qu.:-0.8384   1st Qu.:-0.9202  
##  Median :2.888   Median : 0.1300   Median :-0.0559   Median :-0.4707  
##  Mean   :2.807   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.:2.986   3rd Qu.: 0.7932   3rd Qu.: 0.6707   3rd Qu.: 0.9374  
##  Max.   :3.235   Max.   : 1.7029   Max.   : 3.0184   Max.   : 1.9621

Esta estandarización permite comparar los efectos relativos de cada variable explicativa y hace al intercepto completamente interpretable.

Definimos las distribuciones previas de los parámetros

Por medio de la librería (brms) y su función get_prior, podemos definir cuales son las distribución previas que podemos usar para cada uno de los parámetros, de la siguiente manera:

library(brms)

atra.prior=get_prior(formula=logt_cor~lrc+t_sus+t_am,data=atra.s,family=gaussian)
atra.prior

Como resultado, esta función nos indica los nombres de los coeficientes que corresponden a las pendientes parciales de cada variable explicativa (class: b), el intercepto del modelo y la desviación estándar del modelo estadístico (sigma).

Seguido a esto especifiqué las distribuciones previas para todos los parámetros. Reflajando con esto el conocimiento previo acumulado del caso de estudio reflejado en un conjunto de valores plausibles de los parámetros de un modelo estadístico a ajustar. De la siguiente manera:

definí de manera arbitraria la previa del intercepto con una distribución normal de media = 0 y varianza = 2 indicando así, que el 95% de los valores de ese intercepto de t_cor estarán dentro de -3.5;3.5 lo cual equivalen a los exponenciales [exp(-3.5) = 0.03; exp(3.5) = 33.1], definiéndose así, el intervalo amplio de los valores plausibles de la media de la variable respuesta antes de ver los datos.

Asimismo, se realizó para los valores de las pendientes y de las previas de la desviación estándar de la variable de respuesta (t_cor). Todo esto se especifica y se grafica como sigue:

library(ggplot2)
library(cowplot)
int.a2=ggplot(data.frame(x=c(-3.5,3.5)),aes(x))+
  stat_function(fun = dnorm,n=1000,args = list(mean=0,sd=2))
slope.a2=ggplot(data.frame(x=c(-2,2)),aes(x))+
  stat_function(fun = dnorm,n=1000,args = list(mean=0,sd=1))
sd.a2=ggplot(data.frame(x=c(0,5)),aes(x))+
  stat_function(fun = dcauchy,n=1000,args = list(location=0,scale=3))

plot_grid(int.a2,slope.a2,sd.a2,ncol = 3,labels = LETTERS[1:3],
          align = "hv",label_x = 0.90,label_y = 0.95)

Luego de graficar para visualizar las distribuciones previamente definidas, creamos el objeto que contiene a esas distribuciones previas con la ayuda de la función set_prior de la librería (brms). de la siguiente manera:

atra.prior.2 = c(set_prior("normal(0,2)",class = "Intercept"),
                 set_prior("normal(0,1)",class = "b"),
                 set_prior("cauchy(0,3)",class = "sigma"))
atra.prior.2

4-) Ajustamos el modelo

Habiendo definido las distribuciones previas de los parámetros en el objeto y paso anterior, utilicé la función principal de la librería (brms), la función brms, con la cual, generamos el ajuste del modelo de la siguiente manera:

atra2.brms=brm(formula = logt_cor~lrc+t_sus+t_am, data = atra.s, family = gaussian, prior = atra.prior.2, warmup = 1000,chains = 3,iter = 2000,thin = 3,future = T)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch x86_64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/x86_64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 2.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.029 seconds (Warm-up)
## Chain 1:                0.024 seconds (Sampling)
## Chain 1:                0.053 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 7e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.033 seconds (Warm-up)
## Chain 2:                0.024 seconds (Sampling)
## Chain 2:                0.057 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 8e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.026 seconds (Warm-up)
## Chain 3:                0.025 seconds (Sampling)
## Chain 3:                0.051 seconds (Total)
## Chain 3:
atra2.brms
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: logt_cor ~ lrc + t_sus + t_am 
##    Data: atra.s (Number of observations: 120) 
##   Draws: 3 chains, each with iter = 2000; warmup = 1000; thin = 3;
##          total post-warmup draws = 1002
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     2.81      0.01     2.79     2.82 1.00      918      781
## lrc          -0.02      0.01    -0.04    -0.01 1.00      959      850
## t_sus         0.17      0.01     0.15     0.20 1.00      940      949
## t_am          0.05      0.01     0.02     0.07 1.00      916      991
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.09      0.01     0.08     0.11 1.00      973      860
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

  • Primero, obtenemos el resumen del ajuste del modelo con los parámetros definidos en la formula.

  • Luego, tenemos los estadísticos descriptivos de las distribuciones posteriores marginales para cada uno de los parámetros, lo que ajusta el modelo es una distribución posterior que tiene 4 parámetros: 1 intercepto, 2 pendientes y 1 desviación estándar de t_cor.

  • Ese objeto de 4 dimensiones no es visualizable, por lo tanto, tenemos las distribuciones marginales que son las proyecciones unidimensionales de dicha distribución de 4 dimensiones.

  • Los estadísticos descriptivos son la media del intercepto (Estimate), el Est. error (sd) de esa distribución posterior del intercepto, los cuantiles de 2.5 y 97.5% que define el intervalo de credibilidad del 95%, luego tenemos el estadístico Rhat (“R_sombrero”) es un métrico de la convergencia de las 3 cadenas que utilizamos a una distribución estacionaria posterior, el estadístico Bulk_ESS (Effective Sample Size) y Tail_ESS, nos dan una idea de cuantos valores estadísticamente independientes obtenemos en cada una de las cadenas para cada uno de los parámetros.

Ahora, vamos a verificar la convergencia del algoritmo implementado mediante los siguientes gráficos:

variables(atra2.brms)
## [1] "b_Intercept" "b_lrc"       "b_t_sus"     "b_t_am"      "sigma"      
## [6] "Intercept"   "lprior"      "lp__"

visualizamos las variables del modelo con la función (variables), las seleccionamos, y generamos un gráfico en donde tenemos las distribuciones posteriores de cada parámetro mediante la función (mcmc_dens_overlay) integrada en la librería (bayesplot)

library(bayesplot)
par.names.a2=variables(atra2.brms)[1:5]

dis.names.a2=mcmc_dens_overlay(atra2.brms,pars = par.names.a2)+facet_text(on=T)+theme_bw()

dis.names.a2

Seguido a esto, podemos generar un gráfico de trazas las cuales son la secuencia de valores aceptados para cada cadena para cada uno de los parámetros, la idea es visualizar que las tres cadenas estén razonablemente bien mezcladas. Esto lo logré con la función (mcmc_trace) de la misma librería anterior.

tras.names.a2=mcmc_trace(atra2.brms,size = 0.3,
                         pars = par.names.a2)+facet_text(on=F)+theme_bw()

tras.names.a2

Luego, podemos graficar la autocorrelación que tienen los estimados de los parámetros para cada una de las cadenas, lo cual, nos debe arrojar que unas relaciones bajas indicando que los valores estimados en las cadenas son estadísticamente independientes. Lo hice usando la función (mcmc_acf) de la misma librería anterior.

afc.names.a2=mcmc_acf(atra2.brms,pars = par.names.a2)+scale_x_continuous(limits = c(2,10))+
  scale_y_continuous(limits = c(-0.2,0.2))
## 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.
afc.names.a2
## Warning: Removed 180 rows containing missing values or values outside the scale range
## (`geom_segment()`).
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

Finalmente, visualizamos las distribuciones posteriores de los parámetros usando la función (mcmc_areas) de la librería (bayesplot) después de verificar que las 3 cadenas convergieron a la misma distribución, de la siguiente manera:

library(ggbreak)
dp.names.a2=mcmc_areas(atra2.brms,pars = par.names.a2,
                       point_est = "mean",prob = 0.9)+scale_x_break(c(0.30,2.75))

dp.names.a2

Gelman et al. 2018 proponen un métrico que evalúa la calidad del ajuste con un R2 de la distribución posterior, así:

R2 = data.frame(R2=bayes_R2(atra2.brms,summary = F))
summary(R2)
##        R2        
##  Min.   :0.7858  
##  1st Qu.:0.8390  
##  Median :0.8462  
##  Mean   :0.8444  
##  3rd Qu.:0.8525  
##  Max.   :0.8656
ggplot(data = R2, aes(x=R2)) + geom_density(size=1)+labs(x=expression(paste("R"^2)), y="Dens.probabilidad")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

5-) Evaluamos la calidad del ajuste

vamos a calcular los residuales de tipo Pearson, para cada dato vamos a tener 1000 valores (distribución de residuales) para obtener una distribución de residuales para cada dato, así:

res.atractus.brms=as.data.frame(residuals(atra2.brms, type = "pearson",ndraws = 1000,summary = T))
## Warning: Type 'pearson' is deprecated and will be removed in the future.
str(res.atractus.brms)
## 'data.frame':    120 obs. of  4 variables:
##  $ Estimate : num  1.056 1.766 2.071 -0.243 0.182 ...
##  $ Est.Error: num  1.017 1.04 0.934 1.088 0.963 ...
##  $ Q2.5     : num  -0.927 -0.401 0.351 -2.414 -1.807 ...
##  $ Q97.5    : num  2.98 3.8 3.94 1.82 2.01 ...
fit.atra2.brms=as.data.frame(fitted(atra2.brms,scale = "linear",summary = T,ndraws = 1000))

fit.atra2.brms
  • Lo que obtenemos con esto entonces primero es un Estimado que es la media de los 1000 valores generados para cada dato, el error estándar para cada uno de las 120 observaciones y los intervalos de credibilidad del 95% para cada uno de los residuales.

Luego ajustamos los valores y graficamos así:

library("opdisDownsampling")
library("qqplotr")

atra.s$resid=res.atractus.brms$Estimate
atra.s$resid.Q2.5=res.atractus.brms$Q2.5
atra.s$resid.Q97.5=res.atractus.brms$Q97.5
atra.s$fit=fit.atra2.brms$Estimate
atra.s$fit.Q2.5=fit.atra2.brms$Q2.5
atra.s$fit.Q97.5=fit.atra2.brms$Q97.5

res.fit.atra.brms=ggplot(data = atra.s, aes(x=fit,y=resid))+
  geom_point(col="black",size=2)+
  geom_hline(yintercept = 0,linetype = 2, size =1.1)

res.t_sus.atra.brms=ggplot(data = atra.s, aes(x=t_sus,y=resid))+
  geom_point(col="black",size=2)+
  geom_hline(yintercept = 0,linetype=2, size=1.1)

res.t_am.atra.brms=ggplot(data = atra.s, aes(x=t_am,y=resid))+
  geom_point(col="black",size=2)+
  geom_hline(yintercept = 0,linetype=2, size=1.1)

qq.atra.brms=ggplot(data = atra.s, mapping = aes(sample=resid))+
  stat_qq_point()+stat_qq_line()+stat_qq_band(alpha=0.3)+
  labs(x="Theoretical Quantiles", y="Sample Quantiles")

plot_grid(res.fit.atra.brms,res.t_sus.atra.brms,res.t_am.atra.brms,qq.atra.brms,
          ncol = 3,labels = LETTERS[1:4],align = "hv",label_x = 0.90,label_y = 0.95)

En este paso anterior calculamos los valores ajustados generando nuevamente una distribución posterior de 1000 valores para cada uno de los valores ajustados en el objeto (fit.atra2.brms), dentro de un dataframe que generámos anteriormente (atra.s) donde tenemos a las variables transformadas y escaladas, vamos a incluir 6 columnas más de estos residuales (medias de los residuales, medias de los valores predichos y los intervalos de credibilidad) para generar las gráficas.

Estos gráficos de residuales entonces son:

1- Los residuales contra los valores predichos 2- la media delos residuales contra cada valor de la variable explicativa 3- un qqplot para verificar el ajuste de la normalidad de los residuales

lo que observamos es que hay una distribución razonablemente al azar de los residuales contra los valores ajustados con una amplia dispersión como era de esperarse, lo cual es un indicativo de un buen ajuste de los datos.

Podemos finalizar graficando las curvas condicionales así:

atra2.brms.cond.eff=conditional_effects(atra2.brms)


t_sus.cond=ggplot(data = atra2.brms.cond.eff$t_sus,aes(x=t_sus,y=estimate__))+
  geom_line(size=1)+geom_ribbon(aes(ymin = lower__,ymax=upper__),fill="grey70",alpha=0.5)

t_sus.cond

t_am.cond=ggplot(data = atra2.brms.cond.eff$t_am,aes(t_am,y=estimate__))+
  geom_line(size=1)+geom_ribbon(aes(ymin = lower__,ymax=upper__),fill="grey70",alpha=0.5)

t_am.cond

lrc.cond=ggplot(data = atra2.brms.cond.eff$lrc,aes(x=lrc,y=estimate__))+
  geom_line(size=1)+geom_ribbon(aes(ymin = lower__,ymax=upper__),fill="grey70",alpha=0.5)

lrc.cond

Estas curvas condicionales son las relaciones predichas para cada variable explicativa, con lo que vamos a obtener un gráfico de la variable de respuesta (Estimate__) contra cada variable explicativa, podemos entonces así, vizualizar las pendientes parciales desplegadas antes en el summary del modelo.

CONCLUSIÓN

Aprendí mucho con este ejercicio, el resultado estuvo de igual manera acorde a lo que habíamos obtenido previamente usando regresiones lineales simples y pruebas de anova, la serpiente en cuestión es TERMOCONFORMISTA y TIGMOTÉRMICA aunque hay que anotar que tanto la Temperatura ambiental (t_am) como la Temperatura del sustrato (t_sus) estuvieron altamente correlacionadas, esto debido claramente a que son variables anidadas. Sin embargo, desde la perspectiva de la serpiente aunque estadísticamente las dos variables sean redundates, ellas prefieren mantener su temperatura usando su microhábitat (bajo rocas) como fuente de calor.