Ejemplo 01: Treinta aleaciones del tipo 90/10 Cu-Ni, cada una con un contenido específico de hierro son estudiadas bajo un proceso de corrosión.Tras un período de 60 días se obtiene la pérdida de peso (en miligramos al cuadrado por decímetro y día) de cada una de las aleaciones debido al proceso de corrosión.El objetivo es estudiar el nivel de corrosión en función del contenido de hierro.A continuación se presenta el banco de datos y se realiza la primera inspección gráfica.

hierro <- c(0.01, 0.48, 0.71, 0.95, 1.19, 0.01, 0.48, 1.44, 0.71, 
            1.96, 0.01, 1.44, 1.96)
peso <- c(127.6, 124, 110.8, 103.9, 101.5, 130.1, 122, 92.3, 113.1, 
          83.7, 128, 91.4, 86.2)
corrosion <- data.frame(hierro,peso)
library(ggplot2)
ggplot(corrosion, aes(x = hierro, y = peso)) +
  geom_point() +
  labs(x = "Contenido en hierro", y = "Pérdida de peso")

Conclusión: Se observa cómo al ir aumentando el contenido en hierro de la aleación disminuye linealmente la pérdida de peso.El modelo estadístico que propongamos deberá ser capaz de explicar dicho comportamiento.

Ejemplo 2: Queremos estudiar la relación existente entre la concentración de madera contenida en la pulpa a partir de la que se elabora papel (madera), y la resistencia (tension, en términos de tensión que soporta) del papel resultante. El objetivo del análisis es describir la tendencia observada. A continuación se presenta el banco de datos y se realiza la primera inspección gráfica.

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")

Conclusión: Podemos ver cómo la resistencia del papel crece al aumentar la concentración de madera hasta llegar a valores de 9 y disminuye a partir de ese valor. En este caso la relación apreciada es de tipo parabólico (descrita por una parábola). Este hecho se debe tener en cuenta en la propuesta de un modelo preliminar.

Ejemplo 3: Se ha realizado un experimento para tratar de conocer la viscosidad de cierto compuesto en función de la cantidad de un tipo der aceite que se usa en su fabricación. Se asume una relación de tipo lineal entre la viscosidad y la cantidad de aceite utilizada.

aceite <- c(0, 12, 24, 36, 48, 60, 0, 12, 24, 36, 48, 60, 0, 12, 24,
            36, 48, 60, 12, 24, 36, 48, 60)
viscosidad <- c(26, 38, 50, 76, 108, 157, 17, 26, 37, 53, 83, 124, 13,
                20, 27, 37, 57, 87, 15, 22, 27, 41, 63)
aceites<-data.frame(aceite, viscosidad)
ggplot(aceites, aes(x = aceite, y = viscosidad)) +
  geom_point() +
  labs(x = "Cantidad de aceite", y = "Viscosidad")

Estimación R Tomaremos el ejemplo 1 como referencia,para los datos de Corrosión se propone un modelo de regresión lineal simple para estudiar la relación entre la pérdida de peso debida a la corrosión y el contenido de hierro de la forma siguiente: # peso = Bo + B1 * hierro + error

# Ajuste del modelo
fit <- lm(peso ~ hierro, data = corrosion)
# Solución con tidy
#install.packages("tidymodels")
library(tidymodels)
## ── Attaching packages ────────────────────────────────────── tidymodels 1.1.0 ──
## ✔ broom        1.0.5     ✔ rsample      1.1.1
## ✔ dials        1.2.0     ✔ tibble       3.2.1
## ✔ dplyr        1.1.2     ✔ tidyr        1.3.0
## ✔ infer        1.0.4     ✔ tune         1.1.1
## ✔ modeldata    1.1.0     ✔ workflows    1.1.3
## ✔ parsnip      1.1.0     ✔ workflowsets 1.0.1
## ✔ purrr        1.0.1     ✔ yardstick    1.2.0
## ✔ recipes      1.0.6
## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ purrr::discard() masks scales::discard()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ✖ recipes::step()  masks stats::step()
## • Learn how to get started at https://www.tidymodels.org/start/
tidy(fit)
# Solución con glm_coef
#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: ggridges
## 
## New to ggformula?  Try the tutorials: 
##  learnr::run_tutorial("introduction", package = "ggformula")
##  learnr::run_tutorial("refining", package = "ggformula")
## Loading required package: gtsummary
## 
## Attaching package: 'gtsummary'
## The following objects are masked from 'package:recipes':
## 
##     all_double, all_factor, all_integer, all_logical, all_numeric
## Loading required package: huxtable
## 
## Attaching package: 'huxtable'
## The following object is masked from 'package:gtsummary':
## 
##     as_flextable
## The following object is masked from 'package:dplyr':
## 
##     add_rownames
## The following object is masked from 'package:scales':
## 
##     number_format
## The following object is masked from 'package:ggplot2':
## 
##     theme_grey
## Loading required package: magrittr
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:tidyr':
## 
##     extract
## The following object is masked from 'package:purrr':
## 
##     set_names
glm_coef(fit)
ParameterCoefficientPr(>|t|)
(Intercept)129.79 (126.7, 132.87)< 0.001
hierro-24.02 (-26.84, -21.2)< 0.001
# Por lo tanto la forma del modelo seria la siguiente: 
# Gráfico del ajuste
sjPlot::plot_model(fit,"pred", terms = ~hierro, 
           ci.lvl = NA, 
           show.data = TRUE, 
           axis.title = c("Contenido en hierro", "Peso"),
           title = " ")

Predicción del modelo En primer lugar recuperamos los datos y modelo obtenido para los datos de corrosión.Estamos interesados en conocer la pérdida de peso medio (estimación de la media) y la pérdida de peso específica (predicción de una futura observación) para una barra en particular cuando el contenido en hierro es de 0.5, 1, y 1.5.

# cargamos datos de predicción
newpred <- data.frame(hierro = c(0.5, 1, 0.5))
# Predicción para la media de la respuesta
# Opción interval = "confidence" 
newdata <- data.frame(newpred, 
                      predict(fit, 
                              newpred, 
                              interval = "confidence")) 
round(newdata, 2)
hierrofitlwrupr
0.5118116120
1  106104108
0.5118116120
# cargamos datos de predicción
newpred <- data.frame(hierro = c(0.5, 1, 0.5))
# Predicción de un único valor
# Opción interval = "prediction" 
newdata <- data.frame(newpred, 
                      predict(fit, 
                              newpred, 
                              interval = "prediction")) 
round(newdata, 2)
hierrofitlwrupr
0.5118111  125
1  10698.8113
0.5118111  125

Como ya habíamos comentado la estimación que obtenemos es la misma pero los intervalos de confianza son más amplios.

# Gráfico del ajuste intervalo de confianza 95%
sjPlot::plot_model(fit, "pred", terms = ~hierro, 
           ci.lvl = 0.95, 
           show.data = TRUE, 
           axis.title = c("Contenido en hierro", "Peso"),
           title = " ")