Marcelo Cardillo IMHICIHU-CONICET Presentado en el 3er Congreso Argentino de Estudios Líticos. Puerto Madryn.2025. Argentina

Palabras clave: Asociación, Causalidad, Diagramas Acíclicos Dirigidos, simulacion

Establecer relaciones causales (causa-efecto entre fenómenos) es a menudo el objetivo último de una investigación, aunque en estudios observacionales y naturalistas suele ser factible identificar sólo asociaciones. La experimentación es el método ideal para establecer causalidad, pero en estudios sin control directo de las variables intervinientes, las posibilidades se limitan. Aquí proponemos que los modelos causales, en particular los diagramas acíclicos dirigidos (DAGs), permiten clarificar relaciones causa-efecto al identificar y controlar factores de confusión en análisis tanto experimentales como observacionales. Asimismo, sirven de base para generar modelos generativos y explorar asi nuestros supuestos experimentales. Aqui presento el codigo empleado en esta presentacion

#Paquetes para la generacion del DAG  y la simulacion
library(dagitty)
library(ggdag)
## 
## Adjuntando el paquete: 'ggdag'
## The following object is masked from 'package:stats':
## 
##     filter
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
library(truncnorm)
## Warning: package 'truncnorm' was built under R version 4.4.3
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Ejemplo de lascas en un experimento de pisoteo

#DAG
DAG_Pisoteo <- dagify(
  Fractura~Espesor,
  Espesor~ Peso,
  Peso ~ Sustrato,
  Sustrato ~ Calzado,
  Materia_prima ~ Pisoteo,
  Calzado~Materia_prima,
  exposure = "Pisoteo",
  outcome = "Fractura",
  latent = "Peso",
  labels = c(
    "Pisoteo" = "Pisoteo",
    "Calzado" = "Calzado",
    "Sustrato" = "Sustrato",
    "Espesor" = "Espesor",
    "Materia_prima" = "Materia_prima",
    "Fractura" = "Fractura",
    "Peso" = "Peso"))

El pisoteo es la variable de exposicion, mientras que la fractura es la variable respuesta o resultado de la exposicion (outcome). Las demas variables, median entre exposicion y resultado, sugiriendo que el pisoteo no tiene un efecto directo sobre la Fractura

# Asignamos coordenadas manualmente (opcional para obtener mismo formato)
coordinates(DAG_Pisoteo) <- list(
  x = c(Pisoteo=0, Materia_prima=1, Calzado=2, Sustrato=3, Peso=4, Espesor=5, Fractura=6),
  y = c(Pisoteo=0, Materia_prima=1, Calzado=2, Sustrato=3, Peso=2, Espesor=1, Fractura=0)
)

plot(DAG_Pisoteo)

En este caso, el pisoteo no actúa como una variable sino como una constante en el experimento; Materia prima, Calzado, Sustrato, Peso y Espesor transfieren el efecto del pisoteo, por lo tanto, no estimamos su efecto (ya sea total o indirecto), pero sí podemos modelar las relaciones causales condicionales entre las demás variables y la fractura en el contexto de pisoteo usando el DAG como base para un modelo generativo

###Simulacion de variables del DAG para explorar la operatividad del modelo experimental bajo distintas condiciones

set.seed(42)#para repetibilidad

n =1000

# 1. Variables categóricas
Calzado = rbinom(n, 1, 0.5)     # 1 = duro, 0 = blando
Sustrato = rbinom(n, 1, 0.5)    # 1 = duro, 0 = blando

# 2. Espesor (normal truncada entre 1.5 y 5 mm)

Espesor =rtruncnorm(n, a = 1.5, b = 5, mean = 3, sd = 0.7)

# 3. Log-odds base (blando/blando, espesor medio): 30% → logit(0.3)
logit_base =qlogis(0.3)

# 4. Efectos:
# Calzado duro → duplica probabilidad: logit(0.6) - logit(0.3)
beta_calzado= qlogis(0.6) - qlogis(0.3)

# Sustrato duro → +30% (logit(0.3*1.3) - logit(0.3))
beta_sustrato= qlogis(0.3 * 1.3) - qlogis(0.3)

# Espesor: cada mm reduce 5% → -log(0.95) ≈ -0.051
beta_espesor = log(0.7)

# 5. Línea predictiva y probabilidad
linpred =logit_base + beta_calzado * Calzado + beta_sustrato * Sustrato + beta_espesor * Espesor
p_fractura= plogis(linpred)

# 6. Resultado (binario)
Fractura=rbinom(n, 1, p_fractura)

# 7. Dataset
Modelo_Pisoteo= data.frame(Calzado, Sustrato, Espesor, Fractura)

En esta simulación se busca ejemplificar cómo el diagrama acíclico dirigido (DAG) pueden guiar el diseño y análisis de experimental. Se generaron 1.000 observaciones simuladas, considerando variables cuantitativas y categoricas. Las variables categóricas representan el tipo de calzado (duro o blando) y el sustrato (duro o blando). Luego, la variable continua Espesor, simulada como una distribución normal truncada entre 1.5 y 5 mm para reflejar rangos realistas de guijarros o artefactos. La probabilidad de fractura se calcula usando un modelo logístico, partiendo de una probabilidad base de 0.3 para el escenario de calzado blando sobre sustrato blando con espesor medio. Se incorporan efectos específicos: calzado duro incrementa la probabilidad de fractura al doble, sustrato duro aumenta la probabilidad en un 30% y cada milímetro adicional de espesor reduce la probabilidad en aproximadamente un 5%. La línea predictiva se construye sumando estos efectos en la escala del logit, y luego se transforma a probabilidades mediante la función logística. Finalmente, la variable binaria Fractura se genera a partir de estas probabilidades, produciendo un conjunto de datos que puede usarse para ilustrar cómo las relaciones causales hipotéticas entre variables físicas y de contexto determinan la ocurrencia de fracturas.

head(Modelo_Pisoteo)
##Modelo logistico

mod1=glm(Fractura ~ Espesor+Sustrato+Calzado,family = "binomial",data=Modelo_Pisoteo)

summary(mod1)
## 
## Call:
## glm(formula = Fractura ~ Espesor + Sustrato + Calzado, family = "binomial", 
##     data = Modelo_Pisoteo)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.6598     0.3658  -1.804 0.071284 .  
## Espesor      -0.3968     0.1157  -3.430 0.000604 ***
## Sustrato      0.3775     0.1495   2.526 0.011544 *  
## Calzado       1.2045     0.1526   7.891    3e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1166.5  on 999  degrees of freedom
## Residual deviance: 1082.9  on 996  degrees of freedom
## AIC: 1090.9
## 
## Number of Fisher Scoring iterations: 4

El modelo logístico describe la probabilidad de fractura en función de Espesor, Sustrato y Calzado. El intercept (-0.66) representa el logit de la probabilidad de fractura para la categoría base: calzado blando. El coeficiente de Espesor es -0.397, indicando que a medida que aumenta el espesor del material, la probabilidad de fractura disminuye significativamente; cada milímetro adicional reduce el logit de fractura en 0.39 y en escala de probabilidades, esto representa una caída aproximada del 32% en el odds ratio de fractura por cada mm extra, (exp(0.397)-1)100. El coeficiente de Sustrato es 0.378, mostrando que los sustratos duros incrementan significativamente la probabilidad de fractura comparado con sustratos blandos, aumentando el odds ratio en aproximadamente un 46% (exp(0.378)-1)100.

Finalmente, el coeficiente de Calzado es 1.205, lo que indica un efecto muy fuerte y altamente significativo: usar calzado duro multiplica el odds ratio alrededor de 3.34 veces respecto al calzado, exp(1.205).

##Grafico de la regresion logistica para espesory calzado
Esp = ggplot(Modelo_Pisoteo, aes(x = Espesor, y = Fractura, color = factor(Calzado))) +
  # Puntos observados con jitter vertical
  geom_jitter(height = 0.02, width = 0, alpha = 0.5, size = 2) +
  # Curva ajustada del modelo logístico
  geom_smooth(method = "glm", method.args = list(family = "binomial"),
              se = TRUE, formula = y ~ x, linewidth = 1.2) +
  labs(title = "Probabilidad de fractura según espesor y tipo de calzado",
       x = "Espesor (mm)",
       y = "Probabilidad de fractura",
       color = "Calzado") +
  theme_minimal(base_size = 14)

plot(Esp)

Sust = ggplot(Modelo_Pisoteo, aes(x = Espesor, y = Fractura, color = factor(Sustrato))) +
  # Puntos observados con pequeño jitter vertical
  geom_jitter(height = 0.02, width = 0, alpha = 0.5, size = 2) +
  # Curva ajustada del modelo logístico
  geom_smooth(method = "glm", method.args = list(family = "binomial"),
              se = TRUE, formula = y ~ x, linewidth = 1.2) +
  labs(title = "Probabilidad de fractura según espesor y tipo de sustrato",
       x = "Espesor (mm)",
       y = "Probabilidad de fractura",
       color = "Sustrato") +
  theme_minimal(base_size = 14)

plot(Sust)

set.seed(123)
n = 1000

# 1. Simulamos los datos
Calzado = rbinom(n, 1, 0.5)         # 1 = duro
Sustrato = rbinom(n, 1, 0.5)        # 1 = duro
Espesor = runif(n, 1, 5)            # entre 1 y 5 mm

# Modelo base (probabilidades)
logit_base = qlogis(0.3)
beta_calzado = qlogis(0.6) - qlogis(0.3)
beta_sustrato = log(1.3)  # efecto moderado
beta_espesor = log(0.7)   # 30% reducción por mm

linpred = logit_base + beta_calzado * Calzado + 
  beta_sustrato * Sustrato + beta_espesor * Espesor

p_fractura = plogis(linpred)
Fractura = rbinom(n, 1, p_fractura)

datos = data.frame(Calzado = factor(Calzado, labels = c("Blando", "Duro")),
                    Sustrato = factor(Sustrato, labels = c("Blando", "Duro")),
                    Espesor,
                    Fractura)

# 2. Analizar el efecto del espesor en rangos crecientes
rangos_max = seq(1.5, 5, by = 0.1)

resultados = data.frame()

for (r in rangos_max) {
  sub = datos %>% filter(Espesor >= 1, Espesor <= r)
  
  if (nrow(sub) > 50) {
    modelo = glm(Fractura ~ Espesor + Calzado + Sustrato,
                  data = sub, family = binomial)
    
    beta = coef(summary(modelo))["Espesor", "Estimate"]
    pval = coef(summary(modelo))["Espesor", "Pr(>|z|)"]
    se = coef(summary(modelo))["Espesor", "Std. Error"]
    
    resultados = rbind(resultados,
                        data.frame(Rango_max = r,
                                   Beta = beta,
                                   Pvalor = pval,
                                   SE = se,
                                   N = nrow(sub)))
  }
}

La simulación muestra como el efecto del espesor del material lítico sobre la probabilidad de fractura puede variar según el rango de observación considerado (rango de variavion en el espesor). Se generan 1.000 observaciones simuladas con tres variables independientes: Calzado (blando o duro), Sustrato (blando o duro) y Espesor (uniformemente distribuido entre 1 y 5 mm). La probabilidad de fractura se calcula mediante un modelo logístico que combina un intercepto base (probabilidad inicial de 0.3 para condiciones blando/blando con espesor medio) con efectos predeterminados: calzado duro aumenta la probabilidad de fractura, sustrato duro tiene un efecto moderado positivo, y cada milímetro adicional de espesor reduce la probabilidad de fractura aproximadamente un 30%. A partir de estas probabilidades se genera la variable binaria Fractura mediante un muestreo binomial. Para estudiar la robustez del efecto del espesor, se ajusta repetidamente el modelo logístico sobre subconjuntos de datos que incluyen rangos de espesor crecientes. En cada ajuste se registran el coeficiente de espesor, su error estándar, el valor-p y el tamaño de la muestra, permitiendo analizar cómo cambian la magnitud y la significancia estadística del efecto a medida que se amplía el rango de observación.

###grafico de la simulacion
resultados = resultados %>%
  mutate(Significativo = ifelse(Pvalor < 0.05, "Sí", "No"))

# Graficamos beta y p-valor con color según significancia
ggplot(resultados, aes(x = Rango_max)) +
  geom_line(aes(y = Beta), color = "blue", linewidth = 1) +
  geom_point(aes(y = Beta, color = Significativo), size = 2) +
  scale_color_manual(values = c("Sí" = "forestgreen", "No" = "gray60")) +
  
  geom_line(aes(y = Pvalor), color = "red", linewidth = 1) +
  geom_point(aes(y = Pvalor), color = "red", size = 2) +
  geom_hline(yintercept = 0.05, linetype = "dashed", color = "darkgray") +
  
  scale_y_continuous(
    name = "Coef. beta de Espesor (azul)",
    sec.axis = sec_axis(~., name = "Valor-p (rojo)")
  ) +
  labs(
    x = "Rango máximo de espesor (desde 1 mm)",
    title = "Efecto de Espesor según el rango analizado",
    color = "Valor-p < 0.05"
  ) +
  theme_minimal(base_size = 14) +
  theme(legend.position = "top")

El grafico muestra como la estimación de un efecto del espesor junto con las demas variables puede depender fuertemente del rango de datos considerado (en este caso rangos de espesor crecientes). Conclusiones: El uso de modelos causales y DAGs nos brinda una herramienta útil para diseñar mejor nuestros experimentos, identificar fuentes de confusión y evaluar la validez de nuestras inferencias, tanto en contextos controlados como en el registro naturalista. La representación gráfica de los modelos explicita además la forma en que pensamos se generan los datos y cuáles son las variables responsables de ello. La generación de escenarios hipotéticos (contrafactuales) puede servirnos tanto para evaluar nuestro modelos bajo distintas condiciones Estas herramientas complementan a las herramientas habituales de inferencia basadas en el modelado estadístico y el contraste de hipótesis