Paquetes y datos

primero cargamos los paquetes

library(tidyverse)
library(piecewiseSEM)
library(janitor)
library(kableExtra)

luego leemos los datos y dejamos los nombres limpios usando clean names del paquete janitor

Data <- readxl::read_excel("HOJA PARA SEM.xlsx", 
    col_types = c("text", "text", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric")) %>% clean_names()

y eliminamos datos faltantes

## Eliminar cuando no hay datos
Data <- Data[complete.cases(Data),]

Generación de variable compuesta para Defensas fisicas

Generamos el modelo explicativo primero:

# run multiple regression
model <- lm(herbivoria ~ dureza_foliar + grosor_foliar, data = Data)

y luego obtenemos los parametros

# get loadings
beta_x1 <- summary(model)$coefficients[2, 1]

beta_x2 <-  summary(model)$coefficients[3, 1]

y con eso generamos la variable compuesta Defensas_fisicas

composite <- beta_x1 * Data$dureza_foliar + beta_x2 * Data$grosor_foliar

Data$Defensas_fisicas <- NA

Data$Defensas_fisicas <- composite

Generación modelo general

Con eso ahora podemos generar el modelo con todas las variables

modelGen <- psem(lm(herbivoria ~ n_de_herbivoros + altura + Defensas_fisicas, data = Data), 
                 lm(n_de_herbivoros ~ n_de_aranas, data = Data))

generamos el grafico

plot(modelGen)

Lo cambiamos a glm poisson por la respuesta en abundancia

modelGen <- psem(glm(herbivoria ~ n_de_herbivoros + altura + Defensas_fisicas, data = Data, "poisson"), 
                 glm(n_de_herbivoros ~ n_de_aranas, data = Data, "poisson"))

y hacemos los modelos Top-down y Bottom-up

TopDown <- psem(glm(herbivoria ~ n_de_herbivoros, data = Data, "poisson"), 
                 glm(n_de_herbivoros ~ n_de_aranas, data = Data, "poisson"))

BottomUp <- psem(glm(herbivoria ~ altura + Defensas_fisicas, data = Data, "poisson"))

y juntamos todo en una tabla para juntar por AICc

SummaryGen <- summary(modelGen)
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |======================================================================| 100%
SummaryTop <- summary(TopDown)
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%
SummaryBottom <- summary(BottomUp)


SummaryGen <- SummaryGen$IC %>% mutate(Modelo = "General")
SummaryTop <- SummaryTop$IC %>% mutate(Modelo = "Top-down")
SummaryBottom <- SummaryBottom$IC %>% mutate(Modelo = "Bottom-up") 

Summary <- list(SummaryGen, SummaryTop, SummaryBottom) %>% 
  reduce(bind_rows) %>% 
  arrange(AICc)

Vemos la tabla por AICc

kbl(Summary) %>% kable_paper()
AIC AICc BIC K n Modelo
6.000 6.126 15.804 3 194 Bottom-up
132.223 135.721 145.294 4 194 Top-down
136.357 141.461 155.964 6 194 General

En base a esto la hipótesis bottom-up tiene mas soporte estadístico

summary(BottomUp)
## 
## Structural Equation Model of BottomUp 
## 
## Call:
##   herbivoria ~ altura + Defensas_fisicas
## 
##     AIC      BIC
##  6.000   15.804
## 
## ---
## Tests of directed separation:
## 
##  No independence claims present. Tests of directed separation not possible.
## 
## Global goodness-of-fit:
## 
##   Fisher's C = 0 with P-value = 1 and on 0 degrees of freedom
## 
## ---
## Coefficients:
## 
##     Response        Predictor Estimate Std.Error  DF Crit.Value P.Value
##   herbivoria           altura  -0.1104    0.0110 191   -10.0641       0
##   herbivoria Defensas_fisicas   0.1386    0.0146 191     9.4660       0
##   Std.Estimate    
##              - ***
##              - ***
## 
##   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05
## 
## ---
## Individual R-squared:
## 
##     Response     method R.squared
##   herbivoria nagelkerke       NaN