Este ejemplo contiene los datos de un experimento agronómico para conocer el efecto de diferentes insecticidas (spray) sobre el número de insectos vivos (count) tras un periodo de tratamiento.

data(InsectSprays)
library(ggplot2)
insecticidas <- InsectSprays
ggplot(insecticidas,aes(x = spray, y = count)) + 
  geom_boxplot()

Se aprecia como el número de insectos vivos es superior cuando usamos los sprays A, B, y F (con medias muy similares). El resto de sprays muestran (medias) una supervivencia inferior.Se observan como dos grupos de sprays, uno con medias altas y otro con medias bajas.

Estimación y bondad de ajuste

# Ajuste del modelo
fit.insecticidas <- lm(count ~ spray, data = insecticidas)

# Inferencia sobre los parámetros del modelo
sjPlot::tab_model(fit.insecticidas, 
          show.r2 = FALSE, 
          show.p = FALSE)
  count
Predictors Estimates CI
(Intercept) 14.50 12.24 – 16.76
spray [B] 0.83 -2.36 – 4.03
spray [C] -12.42 -15.61 – -9.22
spray [D] -9.58 -12.78 – -6.39
spray [E] -11.00 -14.20 – -7.80
spray [F] 2.17 -1.03 – 5.36
Observations 72

Las ecuaciones de estimación de las medias del número de insectos vivos se pueden extraer a partir de la tabla de coeficientes estimados. Sin embargo, podemos obtener las estimaciones y representar las medias de todos los sprays haciendo uso de la función update() quitando el término de interceptación del modelo. De esta forma se prescinde de la restricción de identificabilidad y se obtienen las medias directamente.

# Modelo sin interceptación
m1 <- update(fit.insecticidas, . ~ spray - 1)
# Inferencia sin identificabilidad
sjPlot::tab_model(m1, show.r2 = FALSE, show.p = FALSE)
  count
Predictors Estimates CI
spray [A] 14.50 12.24 – 16.76
spray [B] 15.33 13.07 – 17.59
spray [C] 2.08 -0.18 – 4.34
spray [D] 4.92 2.66 – 7.18
spray [E] 3.50 1.24 – 5.76
spray [F] 16.67 14.41 – 18.93
Observations 72

En esta tabla podemos ver tanto el valor estimado con el intervalo de confianza al 95% del número de insectos vivos con cada uno de los sprays. De esta forma os sprays A, B, y Fson los menos efectivos (medias de insectos vivos más altas), mientras que los otros tres tienen medias significativamente más bajas. El spray C es el que se muestra más efectivo con una media de insectos vivos de 2 (redondeando a entero). Podemos ver gráficamente la solución obtenida mediante.

# Gráfico de medias estimadas
sjPlot::plot_model(m1, show.values = TRUE, vline.color = "yellow")

En cuanto a las medidas de bondad de ajuste utilizaremos le test F de igual de medias paraestablecer si los diferentes sprays considerados tienen un comportamiento similar o al menos hay dos de ellos que se comportan de forma distinta.

Diagnóstico

# Valores de diagnóstico
diagnostico <- fortify(fit.insecticidas)
# Gráfico de residuos vs factor
ggplot(diagnostico,aes(x = spray, y = .stdresid)) + 
  geom_boxplot() +
  geom_hline(yintercept = 0, col = "red")

En el gráfico se observa como la variabilidad dentro de cada nivel del factor parece diferente, aunque los residuos parecen tener un comportamiento lineal alrededor del cero.

# Objeto gráfico
p <- sjPlot::plot_model(fit.insecticidas, "pred", terms = "spray",
                show.stat = TRUE)
# Tabla de estimaciones
p$data

Para obtener el verdadero valor de la media debemos deshacer la transformación utilizadapara el ajuste del modelo. Representamos la media predicha así como un intervalo de predicción para ella. Ordenamos los valores obtenidos de mayor a menor predicción para una mejor visualización.

# Secuencia de valores de predicción
newdata <- data.frame(spray = unique(insecticidas$spray))
# Predicción para la media de la respuesta
newdata <- data.frame(newdata, 
                      predict(fit.insecticidas, newdata, interval = "confidence"))
# Deshacemos la transformación
#install.packages("tidyverse") 
#install.packages("modelr")
#install.packages("dplyr")    
library(tidyverse) 
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ lubridate 1.9.2     ✔ tibble    3.2.1
## ✔ purrr     1.0.1     ✔ tidyr     1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(modelr)
library(dplyr)
newdata <- newdata %>% mutate(prediccion = exp(fit) - 1,
                              lower = exp(lwr) - 1, 
                              upper = exp(upr) - 1)
# Gráfico ordenado del menor al mayor valor en función del tipo de spray
ggplot(newdata, aes(x = fct_reorder(spray, prediccion), y = prediccion)) + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = .1) +
  geom_line() +
  geom_point() +
  labs(x = "Spray", y = "Conteo") +
  coord_flip()
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

Motivación de gráficos Se ha realizado un experimento para comprobar la efectividad de diferentes antídotos (AA, AB, AC y AD) frente a diferentes venenos (VA, VB y VC). Para ello se recoge el tiempo de reacción (tiempo) que cada antídoto tarda en hacer efecto para cada veneno.

tiempo <- c(0.31, 0.45, 0.46, 0.43, 0.36, 0.29, 0.4, 0.23, 
            0.22, 0.21, 0.18, 0.23, 0.82, 1.1, 0.88, 0.72, 0.92, 0.61, 0.49, 
            1.24, 0.3, 0.37, 0.38, 0.29, 0.43, 0.45, 0.63, 0.76, 0.44, 0.35, 
            0.31, 0.4, 0.23, 0.25, 0.24, 0.22, 0.45, 0.71, 0.66, 0.62, 0.56, 
            1.02, 0.71, 0.38, 0.3, 0.35, 0.31, 0.33)
antidoto <- factor(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 
                     2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 
                     3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4), 
                   labels=c("AA", "AB", "AC", "AD"))
veneno <- factor(c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 1, 1, 1, 1, 2, 
                   2, 2, 2, 3, 3, 3, 3, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 
                   3, 3, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3), 
                 labels=c("VA", "VB", "VC"))
venenos <- data.frame(tiempo, antidoto, veneno)

Diagrama de cajas

ggplot(venenos,aes(x = antidoto, y = tiempo, color = veneno)) + 
  geom_boxplot() 

#Gráfico de interacción de medias
ggplot(venenos, 
       aes(x = antidoto, y = tiempo, group = veneno, color = veneno)) +
  stat_summary(fun = mean, geom = "point") +
  stat_summary(fun = mean, geom = "line")