Regresión Múltiple y Polinómica
Marcelo Callao Pimentel
29/06/2023
EJEMPLO 1: Para estimar la producción en madera de un bosque se suele realizar un muestreo previo en el que se realizan una serie de medidas no destructivas. Disponemos de mediciones para 20 árboles, así como el volumen (VOL) de madera que producen una vez cortados. Las variablesconsideradas son: HT o altura en pies, DBH el diámetro del tronco a 4 píes de altura (en pulgadas), D16 el diámetro del tronco a 16 pies de altura (en pulgadas), y VOL el volumen de madera conseguida (en pies cúbicos). El objetivo del análisis es determinar cuál es la relación entre dichas medidas y el volumen de madera, con el fin de poder predecir este último en función de las primeras.
dbh <- c(10.2, 13.72, 15.43, 14.37, 15, 15.02, 15.12, 15.24, 15.24, 15.28, 13.78,
15.67, 15.67, 15.98, 16.5, 16.87, 17.26, 17.28, 17.87, 19.13)
d16 <- c(9.3, 12.1, 13.3, 13.4, 14.2, 12.8, 14, 13.5, 14, 13.8, 13.6, 14,
13.7, 13.9, 14.9, 14.9, 14.3, 14.3, 16.9, 17.3)
ht <- c(89, 90.07, 95.08, 98.03, 99, 91.05, 105.6, 100.8, 94, 93.09, 89, 102,
99, 89.02, 95.09, 95.02, 91.02, 98.06, 96.01, 101)
vol <- c(25.93, 45.87, 56.2, 58.6, 63.36, 46.35, 68.99, 62.91, 58.13, 59.79,
56.2, 66.16, 62.18, 57.01, 65.62, 65.03, 66.74, 73.38, 82.87, 95.71)
bosque <- data.frame(vol, dbh, d16, ht)
bosque# Gráficos parciales
#install.packages("reshape")
library(reshape)
#install.packages("melt")
library(melt)
library(ggplot2)
datacomp = melt(bosque, id.vars = 'vol')
ggplot(datacomp) +
geom_jitter(aes(value, vol, colour = variable)) +
facet_wrap(~variable, scales = "free_x") +
labs(x = "", y = "Volumen") 
Ejemplo 02 Se ha llevado a cabo un experimento para estudiar la concentración presente de un fármaco en el hígado después de sufrir un tratamiento. Se piensa que las variables que pueden influir en la concentración son el peso del cuerpo, el peso del hígado y la dosis de fármaco administrada.
p.cuerpo <- c(176, 176, 190, 176, 200, 167, 188, 195, 176, 165, 158, 148, 149, 163,
170, 186, 146, 181, 149)
p.higado <- c(6.5, 9.5, 9.0, 8.9, 7.2, 8.9, 8.0, 10.0, 8.0, 7.9, 6.9, 7.3, 5.2, 8.4,
7.2, 6.8, 7.3, 9.0, 6.4)
dosis <- c(.88, .88, 1.0, .88, 1.0, .83, .94, .98, .88, .84, .80, .74, .75, .81, .85,
.94, .73, .90, .75)
concen <- c(.42, .25, .56, .23, .23, .32, .37, .41, .33, .38, .27, .36, .21, .28, .34,
.28, .30, .37, .46)
concentracion <- data.frame(p.cuerpo, p.higado, dosis, concen)
concentracion# Gráficos parciales
datacomp = melt(concentracion, id.vars = 'concen')
ggplot(datacomp) +
geom_jitter(aes(value, concen, colour = variable)) +
facet_wrap(~variable, scales = "free_x") +
labs(x = "", y = "Concentración del fármaco")
En este caso ninguno de los gráficos parciales muestra una gran asociación entre la concentración del fármaco y cada una de las predictoras. En todos ellos se aprecia una observación un poco más alejada del resto (concentración > 0.6) que podría ser influyente en la obtención del modelo correspondiente.
Ejemplo 03 Banco de datos de Papel de la unidad anterior, donde ya pudimos ver que la tendencia observada se comportaba más como una parábola (polinomio de grado 2) que como una recta.
madera <- c(1, 1.5, 2, 3, 4, 4.5, 5, 5.5, 6, 6.5, 7, 8, 9, 10, 11, 12, 13, 14, 15)
tension <- c(6.3, 11.1, 20.0, 24, 26.1, 30, 33.8, 34, 38.1, 39.9, 42, 46.1, 53.1,
52, 52.5, 48, 42.8, 27.8, 21.9)
papel <- data.frame(madera, tension)
ggplot(papel, aes(x = madera, y = tension)) +
geom_point() +
labs(x = "Concentración de madera", y = "Resistencia del papel")
A simple vista todas las predictoras tienen un efecto positivo en el volumen de madera obtenido, lo cual es bastante obvio, ya que cuanto más grande sea el árbol se espera que su volumen sea más grande. Sin embargo, parece que el efecto de los diámetros es superior al de la altura del árbol (pendientes más pronunciadas) aunque resulta difícildistinguir que diámetro puede ser más relevante ya que ambos se comportan de forma similar. Podemos confirmar este hecho realizando un análisis de correlación para este banco de datos.
Ahora aprenderemos como ajustar un modelo de RLM para un conjunto de datos
# Ajuste del modelo
fit.bosque <- lm(vol ~ dbh + d16 + ht, data = bosque)
# Inferencia sobre los parámetros del modelo
#install.packages("pubh")
library(pubh)## Loading required package: emmeans## Loading required package: ggformula## Loading required package: ggstance##
## Attaching package: 'ggstance'## The following objects are masked from 'package:ggplot2':
##
## geom_errorbarh, GeomErrorbarh## Loading required package: scales## Loading required package: ggridges##
## New to ggformula? Try the tutorials:
## learnr::run_tutorial("introduction", package = "ggformula")
## learnr::run_tutorial("refining", package = "ggformula")## Loading required package: gtsummary## #StandWithUkraine## Loading required package: huxtable##
## Attaching package: 'huxtable'## The following object is masked from 'package:gtsummary':
##
## as_flextable## The following object is masked from 'package:scales':
##
## number_format## The following object is masked from 'package:ggplot2':
##
## theme_grey## Loading required package: magrittrglm_coef(fit.bosque)| Parameter | Coefficient | Pr(>|t|) |
|---|---|---|
| (Intercept) | -108.58 (-138.56, -78.6) | < 0.001 |
| dbh | 1.63 (-0.55, 3.8) | 0.133 |
| d16 | 5.67 (3.12, 8.22) | < 0.001 |
| ht | 0.69 (0.35, 1.04) | < 0.001 |
Obteniendo el modelo de la forma: **vol= -108.58+ 1.63dbh+ 5.67d16+ 0.69*ht**
Representamos gráficamente la estimación e intervalo de confianza de los coeficientesdel modelo para apreciar los efectos descritos:
# Gráfico del ajuste
sjPlot::plot_model(fit.bosque,
show.values = TRUE,
vline.color = "yellow")
# Inferencia sobre los parámetros del modelo
sjPlot::tab_model(fit.bosque,
show.std = TRUE,
show.r2 = FALSE)| vol | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | std. Beta | CI | standardized CI | p |
| (Intercept) | -108.58 | 0.00 | -138.56 – -78.60 | -0.10 – 0.10 | <0.001 |
| dbh | 1.63 | 0.21 | -0.55 – 3.80 | -0.07 – 0.50 | 0.133 |
| d16 | 5.67 | 0.65 | 3.12 – 8.22 | 0.36 – 0.95 | <0.001 |
| ht | 0.69 | 0.24 | 0.35 – 1.04 | 0.12 – 0.36 | 0.001 |
| Observations | 20 | ||||
# Gráfico del ajuste
sjPlot::plot_model(fit.bosque, "pred",
ci.lvl = NA,
show.data = TRUE,
title = "Modelo ajustado")## $dbh
##
## $d16
##
## $ht
# Bondad de ajuste
glance(fit.bosque)| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.959 | 0.951 | 3.1 | 125 | 2.59e-11 | 3 | -48.7 | 107 | 112 | 153 | 16 | 20 |
# Tabla ANOVA
anova(fit.bosque)| Df | Sum Sq | Mean Sq | F value | Pr(>F) |
|---|---|---|---|---|
| 1 | 3.09e+03 | 3.09e+03 | 322 | 5.05e-12 |
| 1 | 332 | 332 | 34.6 | 2.3e-05 |
| 1 | 173 | 173 | 18.1 | 0.000606 |
| 16 | 153 | 9.58 |
Tanto el R2 como el R2 ajustado muestran porcentajes del 95% indicando que el modelo ajustado tiene buena capacidad explicativa. Además, el test F de la regresión resulta significativo (p-valor < 0.05) indicando que las predictoras consideradas pueden ser utilizadas para describir el comportamiento del volumen. Los tests individuales de la tabla ANOVA.