Mostrar los casos habidos y mostrar la tendencia lineal.

df %>% 
  ggplot(aes(x=Fecha, y=Casos)) +
  geom_point(alpha=0.5) +
  geom_smooth(method="glm", method.args=list(family="poisson"), se=TRUE, color="blue") +
  labs(title="Modelo de Poisson ajustado a los casos Viogen",
       x="Fecha",
       y="Casos") +
  theme_minimal() +
  scale_x_date(date_labels="%Y-%m", date_breaks="6 month") +
  theme(axis.text.x = element_text(angle=45, hjust=1))
## `geom_smooth()` using formula = 'y ~ x'

La línea horizontal no entra en la zona sombreada. La tendencia descendente es clara. Eso sí, España ya encabeza las estadísticas de menor violencia contra las mujeres. SI la linea fuese descendiente siempre sería difícil de creer.

Tal vez el descenso no sea constante y vaya por épocas asociado a cualquier cosa: campañas de sensibilización, cambios en la legislación, casos particulares notables,…

# Vamos a hacer una estimación loess para ver si encontramos un patrón interesante...
df %>% 
  ggplot(aes(x=Fecha, y=Casos)) +

  geom_point(alpha=0.5) +
  geom_smooth(method="loess", se=TRUE, color="red") +
  labs(title="Estimación Loess de los casos Viogen",
       x="Fecha",
       y="Casos") +
  theme_minimal() +
  scale_x_date(date_labels="%Y-%m", date_breaks="6 month") +
  theme(axis.text.x = element_text(angle=45, hjust=1))
## `geom_smooth()` using formula = 'y ~ x'

Casi parece haber tres épocas… horizontal, descenso y horizontal de nuevo. Realmente hubo una época de descenso hasta la situación actual.

Permitamos algo de overfitting por si surge algo interesante.

desde=ymd("2009-11-01")
hasta=ymd("2013-05-01")

df %>% 
  ggplot(aes(x=Fecha, y=Casos)) +
  geom_point(alpha=0.5) +
  geom_smooth(method="loess", se=TRUE, color="red", span = 0.3) +
  labs(title="Estimación Loess de los casos Viogen",
       x="Fecha",
       y="Casos") +
  theme_minimal() +
  geom_vline(xintercept=desde, linetype="dashed", color="blue") +
  geom_vline(xintercept=hasta, linetype="dashed", color="blue") +
  scale_x_date(date_labels="%Y-%m", date_breaks="6 month") +
  theme(axis.text.x = element_text(angle=45, hjust=1))
## `geom_smooth()` using formula = 'y ~ x'

Más o menos lo mismo, pero más ruidoso.

Un poquitín menos de overfitting ahora:

df %>% 
  ggplot(aes(x=Fecha, y=Casos)) +
  geom_point(alpha=0.5) +
  geom_smooth(method="loess", se=TRUE, color="red", span = 0.4) +
  labs(title="Estimación Loess de los casos Viogen",
       x="Fecha",
       y="Casos") +
    geom_vline(xintercept=desde, linetype="dashed", color="blue") +
  geom_vline(xintercept=hasta, linetype="dashed", color="blue") +
  theme_minimal() +
  scale_x_date(date_labels="%Y-%m", date_breaks="6 month") +
  theme(axis.text.x = element_text(angle=45, hjust=1))
## `geom_smooth()` using formula = 'y ~ x'

Vamos a testear el modelo lineal de poisson que hemos ajustado antes.

modelo <- glm(Casos ~ Fecha, data = df, family = "poisson")
summary(modelo)
## 
## Call:
## glm(formula = Casos ~ Fecha, family = "poisson", data = df)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.432e+00  1.865e-01  13.037  < 2e-16 ***
## Fecha       -5.319e-05  1.163e-05  -4.573 4.81e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 341.61  on 270  degrees of freedom
## Residual deviance: 320.59  on 269  degrees of freedom
## AIC: 1221.5
## 
## Number of Fisher Scoring iterations: 4

El efecto del tiempo es significativo en el descenso de la viogen.

Voy a crear un modelo de 3 épocas, dos horizontales y una de descenso:

#Alterar el modelo anterior de forma que haya 3 épocas

df <- df %>%
  mutate(Epoca = case_when(
    Fecha < desde ~ "Antes",
     Fecha <= hasta ~ "Entre",
    TRUE ~ "Después"),
    Tiempo=Fecha-desde,
    TiempoEntre= case_when(
    Fecha < desde ~ 0L,
    Fecha <= hasta ~ as.integer(Fecha-desde),
    TRUE ~ as.integer(hasta-desde)),
    Despues= ifelse(Epoca == "Después", 1L, 0L))


modelo_epocas <- glm(Casos ~ TiempoEntre, data = df, family = "poisson")
summary(modelo_epocas)
## 
## Call:
## glm(formula = Casos ~ TiempoEntre, family = "poisson", data = df)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.744e+00  4.284e-02  40.706  < 2e-16 ***
## TiempoEntre -2.177e-04  4.601e-05  -4.731 2.24e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 341.61  on 270  degrees of freedom
## Residual deviance: 319.51  on 269  degrees of freedom
## AIC: 1220.4
## 
## Number of Fisher Scoring iterations: 4

Aparentemente este modelo es una explicación más simple de lo que pasa. Vamos a comparar el modelo que dice que la viogen baja todo el tiempo con el que dice que bajó desde un nivel al actual durante un periodo. Comparemos los dos modelos:

# Calcular AIC
aic_values <- AIC(modelo, modelo_epocas)

# Calcular deviance residual promedio
resid_dev_modelo <- mean(residuals(modelo, type = "deviance")^2)
resid_dev_modelo_epocas <- mean(residuals(modelo_epocas, type = "deviance")^2)

# Calcular pseudo-R² de McFadden
pseudo_r2_modelo <- 1 - (modelo$deviance / modelo$null.deviance)
pseudo_r2_modelo_epocas <- 1 - (modelo_epocas$deviance / modelo_epocas$null.deviance)


# Imprimir resultados
print(aic_values)
##               df      AIC
## modelo         2 1221.485
## modelo_epocas  2 1220.398
cat("\nDeviance residual promedio:\n")
## 
## Deviance residual promedio:
cat("Modelo (Fecha):", resid_dev_modelo, "\n")
## Modelo (Fecha): 1.183002
cat("Modelo (TiempoEntre):", resid_dev_modelo_epocas, "\n")
## Modelo (TiempoEntre): 1.178991
cat("\nPseudo-R² de McFadden:\n")
## 
## Pseudo-R² de McFadden:
cat("Modelo (Fecha):", pseudo_r2_modelo, "\n")
## Modelo (Fecha): 0.06151104
cat("Modelo (TiempoEntre):", pseudo_r2_modelo_epocas, "\n")
## Modelo (TiempoEntre): 0.06469352