Modelo de regresion lineal multiple

Modelo de regresion lineal multiple

Este es el uso de la regresión lineal con múltiples variables, y la ecuación es:

\(Y = b_0 + b_1X_1 + b_2X_2 + b_3X_3 + … + b_nX_n + e\)

\(Y\) y \(b_0\) son los mismos que en el modelo de regresión lineal simple.

\(b_1X_1\) representa el coeficiente de regresión ( b1) sobre la primera variable independiente ( X1).El mismo análisis se aplica a todos los coeficientes y variables de regresión restantes.

\(e\) es el error del modelo (residuales), que define cuánta variación se introduce en el modelo al estimar Y.

Es posible que no siempre obtengamos una línea recta para un caso de regresión múltiple. Sin embargo, podemos controlar la forma de la línea ajustando un modelo más apropiado.

Los supuestos de regresión lineal múltiple

1. Los valores residuales se distribuyen normalmente . Esto puede comprobarse utilizando un gráfico de probabilidad normal o un histograma .

2. Debe haber una relación lineal entre las variables dependientes e independientes . Esto se puede ilustrar mediante diagramas de dispersión que muestran una relación lineal o curvilínea.

3. la multicolinealidad es otro supuesto, lo que significa que las variables independientes no están altamente correlacionadas entre sí. La multicolinealidad dificulta identificar qué variables explican mejor la variable dependiente. Esta suposición se verifica calculando una matriz de correlaciones bivariadas de Pearson entre todas las variables independientes. Si no hay colinealidad en los datos, todos los valores deben ser inferiores a 0,8.

4. La homocedasticidad asume que la varianza de los errores residuales es similar a través del valor de cada variable independiente. Una forma de comprobarlo es a través de un gráfico de los valores pronosticados frente a los valores residuales estandarizados para ver si los puntos se distribuyen por igual entre todos los valores de las variables independientes.

Una guía paso a paso para la regresión lineal múltiple en R

1. Visualizando los datos
2. Desarrollo del modelo
3. Identificar variables estadisticamente significativa
4. Predicción

Modelo 1

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(knitr)
library(readr)
library(ggplot2)
library(patchwork) # PAra varias gráficos en el mismo renglón

## Warning: package 'patchwork' was built under R version 4.2.3

Se llevó a cabo un conjunto de ensayos experimentales con un horno para determinar una forma de predecir el tiempo de cocción, y, a diferentes niveles de ancho del horno, \(x_1\), y a diferentes temperaturas, \(x_2\). Se registraron los siguientes datos:

yp <-c(6.40, 15.05, 18.75, 30.25, 44.85, 48.85, 51.55, 61.50, 100.44, 111.42)
x1 <-c(1.32, 2.69, 3.56, 4.41, 5.35, 6.20, 7.12, 8.87, 9.80, 10.65)
x2 <-c(1.15, 3.40, 4.10, 8.75, 14.82, 15.15, 15.32, 18.18, 35.19, 40.40)
datos<-data.frame(yp, x1, x2)
kable(datos, caption = "Factores que influyen en el tiempo de coccion segun diferentes niveles de ancho del horno y diferentes temperaturas")

Factores que influyen en el tiempo de coccion segun diferentes niveles de ancho del horno y diferentes temperaturas

yp	x1	x2
6.40	1.32	1.15
15.05	2.69	3.40
18.75	3.56	4.10
30.25	4.41	8.75
44.85	5.35	14.82
48.85	6.20	15.15
51.55	7.12	15.32
61.50	8.87	18.18
100.44	9.80	35.19
111.42	10.65	40.40

Variable dependiente \(y\) = tiempo de cocción Variable independiente \(x_1\) = ancho del horno Variable independiente \(x_2\) = diferentes temperaturas

Visualizando los datos

g1 <- ggplot(data = datos, mapping = aes(x = x1, y = yp)) +
  geom_point(color = "forestgreen", size = 2) +
  labs(title  =  'yp ~ x1', x  =  'x1') +
  geom_smooth(method = "lm", se = FALSE, color = "black") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) 

g2 <- ggplot(data = datos, mapping = aes(x = x2, y = yp)) +
  geom_point(color = "orange", size = 2) +
  labs(title  =  'yp ~ x2', x  =  'x2') +
  geom_smooth(method = "lm", se = FALSE, color = "black") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

g1+g2

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Multicolinealidad

¿Qué problemas causa la multicolinealidad? La multicolinealidad causa los siguientes dos tipos básicos de problemas:

1. Las estimaciones de los coeficientes pueden variar enormemente en función de qué otras variables independientes se encuentran en el modelo. Los coeficientes se vuelven muy sensibles a pequeños cambios en el modelo.

2. La multicolinealidad reduce la precisión de los coeficientes estimados, lo que debilita el poder estadístico de su modelo de regresión. Es posible que no pueda confiar en los valores p para identificar variables independientes que sean estadísticamente significativas.

data_independiente <- data.frame(x1,x2)
cor(data_independiente,method = "pearson")

##           x1        x2
## x1 1.0000000 0.9375592
## x2 0.9375592 1.0000000

Podemos notar fuerte correlación porque su valor es superior a 0,8.

La multicolinealidad ocurre cuando las variables independientes en un modelo de regresión están correlacionadas. Esta correlación es un problema porque las variables independientes deberían ser independientes . Si el grado de correlación entre las variables es lo suficientemente alto, puede causar problemas al ajustar el modelo e interpretar los resultados.

resaltaré los problemas que puede causar la multicolinealidad, le mostraré cómo probar su modelo y destacaré algunas formas de resolverlo. En algunos casos, la multicolinealidad no es necesariamente un problema y le mostraré cómo tomar esta determinación.

¿Por qué la multicolinealidad es un problema potencial?

Un objetivo clave del análisis de regresión es aislar la relación entre cada variable independiente y la variable dependiente. La interpretación de un coeficiente de regresión es que representa el cambio medio en la variable dependiente por cada cambio de 1 unidad en una variable independiente cuando se mantienen constantes todas las demás variables independientes . Esa última parte es crucial para nuestra discusión sobre la multicolinealidad.

La idea es que puedes cambiar el valor de una variable independiente y no de las otras. Sin embargo, cuando las variables independientes están correlacionadas, indica que los cambios en una variable están asociados con cambios en otra variable. Cuanto más fuerte es la correlación, más difícil es cambiar una variable sin cambiar otra. Se vuelve difícil para el modelo estimar la relación entre cada variable independiente y la variable dependiente de forma independiente porque las variables independientes tienden a cambiar al unísono.

Hay dos tipos básicos de multicolinealidad:

1. Multicolinealidad estructural : Este tipo ocurre cuando creamos un término modelo usando otros términos. En otras palabras, es un subproducto del modelo que especificamos en lugar de estar presente en los datos mismos. Por ejemplo, si eleva al cuadrado el término X para modelar la curvatura, claramente existe una correlación entre X y X 2 .

2. Multicolinealidad de datos : este tipo de multicolinealidad está presente en los propios datos en lugar de ser un artefacto de nuestro modelo. Es más probable que los experimentos de observación muestren este tipo de multicolinealidad.

.

¿Tengo que corregir la multicolinealidad?

La necesidad de reducir la multicolinealidad depende de su gravedad y de su objetivo principal para su modelo de regresión. Tenga en cuenta los siguientes tres puntos:

1. La gravedad de los problemas aumenta con el grado de multicolinealidad. Por lo tanto, si solo tiene una multicolinealidad moderada, es posible que no necesite resolverla.

2. La multicolinealidad afecta solo a las variables independientes específicas que están correlacionadas. Por lo tanto, si la multicolinealidad no está presente para las variables independientes que le interesan particularmente, es posible que no necesite resolverla. Suponga que su modelo contiene las variables experimentales de interés y algunas variables de control. Si existe una alta multicolinealidad para las variables de control pero no para las variables experimentales, entonces puede interpretar las variables experimentales sin problemas.

3. La multicolinealidad afecta los coeficientes y los valores p , pero no influye en las predicciones, la precisión de las predicciones y las estadísticas de bondad de ajuste . Si su objetivo principal es hacer predicciones y no necesita comprender el rol de cada variable independiente, no necesita reducir la multicolinealidad severa.

Desarrollo de dos modelos

El primer modelo, es un modelo considerando las dos variables independientes \(x_1\) y \(x_2\)

modelo <- lm(formula = yp ~ ., data = datos)

summary(modelo)

## 
## Call:
## lm(formula = yp ~ ., data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8475 -0.3438  0.0043  0.2554  1.1578 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.57723    0.59865   0.964    0.367    
## x1           2.70957    0.19935  13.592 2.75e-06 ***
## x2           2.05033    0.04743  43.227 9.26e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6481 on 7 degrees of freedom
## Multiple R-squared:  0.9997, Adjusted R-squared:  0.9997 
## F-statistic: 1.304e+04 on 2 and 7 DF,  p-value: 3.166e-13

El segundo modelo, es un modelo descartando una de las variables independientes, realizamos el modelo con la variable \(x_2\) la cual es estadisticamente más significativa en el primer modelo

modelo2 <- lm(formula = yp ~ x2, data = datos)

summary(modelo2)

## 
## Call:
## lm(formula = yp ~ x2, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0226 -1.7338 -0.3497  1.0695  5.8668 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.36967    1.61355   4.567  0.00183 ** 
## x2           2.65476    0.08077  32.869 8.01e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.173 on 8 degrees of freedom
## Multiple R-squared:  0.9926, Adjusted R-squared:  0.9917 
## F-statistic:  1080 on 1 and 8 DF,  p-value: 8.005e-10

¿Cómo saber cuál de los dos modelos es mejor?

Una forma de responder a esta pregunta es realizar una prueba de análisis de varianza (ANOVA) de los dos modelos. Contrasta la hipótesis nula ( H0 ), donde las variables que eliminamos anteriormente no tienen significación, contra la hipótesis alternativa ( H1 ) de que esas variables son significativas.

Si el nuevo modelo es una mejora del modelo original, entonces no podemos rechazar H0. Si ese no es el caso, significa que esas variables fueron significativas; por lo tanto rechazamos H0.

# Anova test
anova(modelo, modelo2)

## Analysis of Variance Table
## 
## Model 1: yp ~ x1 + x2
## Model 2: yp ~ x2
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1      7  2.940                                  
## 2      8 80.532 -1   -77.592 184.74 2.745e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Del resultado de ANOVA, observamos que el valor p ( 2.745e-06 ) es muy pequeño (menos de 0.05), por lo que rechazamos la hipótesis nula, lo que significa que el segundo modelo no es una mejora del primero.

Desarrollo del modelo

EL resumen del modelo con la función summary(modelo) identifica que las variables ancho del horno y diferentes temperaturas son estadísticamente significativas dado que presentan valores por debajo de 0.05 en Pr(>|t|)

\(\beta_0 = 0.57723\) \(\beta_1 = 2.70957\) \(\beta_2 = 2.05033\)

Por lo tanto el modelo nos quedaría de la siguiente forma

\(yˆ=\beta_0+\beta_1x_1+\beta_2x_2\)

Descartamos \(b_0\) por no ser estadisticmente significativo

modelo_nuevo <- lm(formula = yp ~ x1 + x2 -1, data = datos)

summary(modelo_nuevo)

## 
## Call:
## lm(formula = yp ~ x1 + x2 - 1, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8103 -0.3698  0.1963  0.3955  1.1807 
## 
## Coefficients:
##    Estimate Std. Error t value Pr(>|t|)    
## x1  2.87003    0.10927   26.27 4.74e-09 ***
## x2  2.02140    0.03657   55.28 1.27e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6452 on 8 degrees of freedom
## Multiple R-squared:  0.9999, Adjusted R-squared:  0.9999 
## F-statistic: 4.188e+04 on 2 and 8 DF,  p-value: < 2.2e-16

Valores residuales se distribuyen normalmente

En un modelo lineal de regresión múltiple, los residuos son la diferencia entre los valores observados y los valores predichos por el modelo para cada uno de los puntos de datos. Matemáticamente, los residuos se calculan como la diferencia entre el valor observado (y) y el valor predicho (ŷ) para cada punto de datos, es decir:

e = y - ŷ

Los residuos son una medida de qué tan bien se ajusta el modelo a los datos. Idealmente, los residuos deberían ser aleatorios y tener una distribución normal con media cero y una varianza constante (homocedasticidad), lo que indica que el modelo captura correctamente la relación entre las variables independientes y la variable dependiente. Sin embargo, si los residuos no se distribuyen normalmente o no son homocedásticos, puede indicar que el modelo no es adecuado para los datos o que hay variables importantes que no se han incluido en el modelo.

# Get the model residuals
model_residuals = modelo_nuevo$residuals

# Plot the residuals
qqnorm(model_residuals)
# Plot the Q-Q line
qqline(model_residuals)

El test de Shapiro-Wilks plantea la hipótesis nula que una muestra proviene de una distribución normal. Eligimos un nivel de significanza, por ejemplo 0,05, y tenemos una hipótesis alternativa que sostiene que la distribución no es normal.

Tenemos:

H0: La distribución es normal

H1: La distribución no es normal

shapiro.test(model_residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  model_residuals
## W = 0.95058, p-value = 0.6754

Vemos que el valor de probabilidad (p) es muy superios a nuestro nivel elegido (0,05), por lo que no rechazamos la hipótesis nula.

Homocedasticidad

La homocedasticidad, también conocida como homogeneidad de varianzas, es un supuesto en los modelos de regresión lineal múltiple que establece que la varianza de los errores o residuos debe ser constante para todos los valores de las variables predictoras o independientes. Esto significa que la dispersión de los residuos alrededor de la línea de regresión debe ser similar en toda la gama de valores de las variables predictoras.

Si no se cumple este supuesto, es decir, si hay una variación no constante en los residuos en diferentes valores de las variables predictoras, entonces se dice que hay heterocedasticidad en el modelo. La heterocedasticidad puede tener efectos negativos en la inferencia estadística y en la precisión de las predicciones del modelo. Por lo tanto, es importante verificar la homocedasticidad en un modelo de regresión lineal múltiple y, en caso de encontrar heterocedasticidad, tomar medidas para remediarla.

par(mfrow = c(2, 2))
plot(modelo_nuevo)

Test de homocedasticidad: La función bptest() en R es un test de Breusch-Pagan para la heterocedasticidad en modelos de regresión. Esta función toma como entrada un modelo de regresión y devuelve el resultado de la prueba de hipótesis para la homocedasticidad de los residuos.

HO: los residuos tienen varianza constante (homocedasticidad)

HA: hay heterocedasticidad en los residuos

El resultado de la prueba incluye el valor del estadístico de prueba (el valor de la prueba de Breusch-Pagan), el valor-p y el número de grados de libertad. Si el valor-p es menor que el nivel de significancia elegido, se rechaza la hipótesis nula de homocedasticidad y se concluye que hay heterocedasticidad en los residuos.

library(lmtest)

## Warning: package 'lmtest' was built under R version 4.2.3

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.2.3

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

bptest(modelo_nuevo)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_nuevo
## BP = 5.7517, df = 1, p-value = 0.01647

Como podemos ver tenemos un p-value = 0.01647, notando que este valor es mayor que 0.05, entonces rechazamos la hipótesis nula y podemos asumir que existe heterocedasticidad

Heterocedasticidad en análisis de regresión

Heterocedasticidad significa dispersión desigual. En el análisis de regresión , hablamos de heteroscedasticidad en el contexto de los residuos o término de error. Específicamente, la heteroscedasticidad es un cambio sistemático en la dispersión de los residuos sobre el rango de valores medidos. La heterocedasticidad es un problema porque la regresión de mínimos cuadrados ordinarios (OLS) asume que todos los residuos se extraen de una población que tiene una varianza constante (homocedasticidad).

Para satisfacer los supuestos de regresión y poder confiar en los resultados, los residuos deben tener una varianza constante.

¿Qué problemas causa la heterocedasticidad?

Como mencioné anteriormente, la regresión lineal asume que la dispersión de los residuos es constante a lo largo de la gráfica. Cada vez que viola una suposición, existe la posibilidad de que no pueda confiar en los resultados estadísticos.

¿Por qué solucionar este problema? Hay dos grandes razones por las que desea la homocedasticidad:

1. Si bien la heteroscedasticidad no genera sesgo en las estimaciones de los coeficientes , sí las hace menos precisas. Una menor precisión aumenta la probabilidad de que las estimaciones del coeficiente estén más alejadas del valor correcto de la población.

2. La heterocedasticidad tiende a producir valores de p que son más pequeños de lo que deberían ser. Este efecto ocurre porque la heteroscedasticidad aumenta la varianza de los coeficientes estimados pero el procedimiento OLS no detecta este aumento.

Corregir la heterocedasticidad

Si puede averiguar el motivo de la heterocedasticidad, es posible que pueda corregirlo y mejorar su modelo. Le mostraré tres enfoques comunes para convertir la heterocedasticidad en homocedasticidad.

1. Redefiniendo las variable
2. regresión ponderada
3. Transformar la variable predictora

Trasformacion de la varaible predictora

La elección de la transformación adecuada para solucionar el supuesto de homocedasticidad en un modelo de regresión lineal múltiple puede ser un proceso iterativo que implica probar diferentes transformaciones y evaluar los resultados. A continuación, te presento algunas técnicas que pueden ser útiles para seleccionar la transformación adecuada:

1. Gráficos de residuos: Los gráficos de residuos pueden ayudarte a identificar patrones en los residuos que sugieren heterocedasticidad. Si observas un patrón en forma de abanico o de embudo en los residuos, puede ser indicativo de heterocedasticidad. En este caso, puedes probar con transformaciones logarítmicas, cuadráticas o raíces cuadradas para solucionar el problema.

2. Pruebas estadísticas: Puedes utilizar pruebas estadísticas como la prueba de Breusch-Pagan o la prueba de White para evaluar la homocedasticidad en los residuos del modelo. Si la prueba indica la presencia de heterocedasticidad, puedes probar diferentes transformaciones hasta que la prueba no indique la presencia de heterocedasticidad.

3. Conocimiento del dominio: La elección de la transformación adecuada también puede basarse en el conocimiento del dominio del problema. Por ejemplo, si estás modelando el precio de una propiedad en función de su tamaño y edad, es posible que desees aplicar una transformación logarítmica a la variable de tamaño para reflejar la relación no lineal entre el tamaño y el precio.

modelo_nuevo_2 <- lm(formula = yp ~ log(x1) + x2 -1, data = datos)

summary(modelo_nuevo_2)

## 
## Call:
## lm(formula = yp ~ log(x1) + x2 - 1, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4546 -0.7961 -0.3458  0.9207  2.5270 
## 
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)    
## log(x1)  7.52724    0.72213   10.42 6.22e-06 ***
## x2       2.34013    0.06306   37.11 3.05e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.578 on 8 degrees of freedom
## Multiple R-squared:  0.9994, Adjusted R-squared:  0.9993 
## F-statistic:  6997 on 2 and 8 DF,  p-value: 1.065e-13

library(lmtest)
bptest(modelo_nuevo_2)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_nuevo_2
## BP = 0.13878, df = 1, p-value = 0.7095

Notando que este valor es mayor que 0.05, entonces no rechazamos la hipótesis nula y podemos asumir que existe homocedasticidad.

Modelo final

\(\beta_0 = 0\) \(\beta_1 = 7.52724\) \(\beta_2 = 2.34013\)

Por lo tanto el modelo nos quedaría de la siguiente forma

\(y=\beta_0+\beta_1x_1+\beta_2x_2\)

\(y= 7.52724x_1+2.34013x_2\)

par(mfrow = c(2, 2))
plot(modelo_nuevo_2)

#calcular los residuos
standard_res <- rstandard (modelo_nuevo_2)

#ver los residuos estandarizados
standard_res

##          1          2          3          4          5          6          7 
##  1.0298111 -0.2384782 -0.2831889 -0.9576857 -1.6441445 -0.2294864  0.6433360 
##          8          9         10 
##  1.7823016  0.7278246 -0.8948154

#column vincula los residuos estandarizados de nuevo al marco de datos original
 final_data <- cbind (datos, standard_res)

 #ordenar residuos estandarizados descendentes
 final_data [ order (-standard_res),]

##        yp    x1    x2 standard_res
## 8   61.50  8.87 18.18    1.7823016
## 1    6.40  1.32  1.15    1.0298111
## 9  100.44  9.80 35.19    0.7278246
## 7   51.55  7.12 15.32    0.6433360
## 6   48.85  6.20 15.15   -0.2294864
## 2   15.05  2.69  3.40   -0.2384782
## 3   18.75  3.56  4.10   -0.2831889
## 10 111.42 10.65 40.40   -0.8948154
## 4   30.25  4.41  8.75   -0.9576857
## 5   44.85  5.35 14.82   -1.6441445

Prediccion

Para 2.10 de ancho del horno y una temperatura de 3.10 , ¿cuánto seria el tiempo de cocción?

nuevo.dato <- data.frame(x1 = 2.10, x2 = 3.10)

prediccion <- predict(modelo_nuevo_2, newdata = nuevo.dato)

paste("La cantidad estimada de tiempo de coccion es:", round(prediccion, 2))

## [1] "La cantidad estimada de tiempo de coccion es: 12.84"

Resumen del modelo 1

1. En este modelo de regresión lineal múltiple se ha realizado una transformación logarítmica a la variable \(x_1\) para cumplir con el supuesto de homocedasticidad en los residuales.

2. La salida indica que el modelo tiene un alto ajuste (Múltiple R-squared: 0.9994, Adjusted R-squared: 0.9993) y que ambas variables independientes, log(x1) y x2, son estadísticamente significativas en el modelo (p-value: < 0.001 para ambas variables).

3. Además, el valor de t para la variable log(x1) es de 10.42, lo que indica que esta variable tiene una influencia significativa en la variable dependiente yp.

4. El modelo también cumple con el supuesto de normalidad en los residuales, ya que no se informan problemas en la distribución de residuos. El modelo tiene un error estándar residual de 1.578, lo que indica una buena precisión en las predicciones.

Modelo 2

Datarium marketing es Un marco de datos que contiene el impacto de tres medios publicitarios (youtube, facebook y periódico) en las ventas. Los datos son el presupuesto publicitario en miles de dólares junto con las ventas. El experimento publicitario se ha repetido 200 veces.

Veamos una descripción general rápida del conjunto de datos para que podamos aplicar el preprocesamiento relevante antes de ajustar el modelo.

library(datarium)
head(marketing)

##   youtube facebook newspaper sales
## 1  276.12    45.36     83.04 26.52
## 2   53.40    47.16     54.12 12.48
## 3   20.64    55.08     83.16 11.16
## 4  181.80    49.56     70.20 22.20
## 5  216.96    12.96     70.08 15.48
## 6   10.44    58.68     90.00  8.64

# Check the dimension
dim(marketing)

## [1] 200   4

De los resultados anteriores, podemos observar que el conjunto de datos tiene 200 observaciones y 4 columnas.

Visualizando los datos

g1 <- ggplot(data = marketing, mapping = aes(x = youtube, y = sales)) +
  geom_point(color = "forestgreen", size = 2) +
  labs(title  =  'sales ~ youtube', x  =  'youtube') +
  geom_smooth(method = "lm", se = FALSE, color = "black") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) 

g2 <- ggplot(data = marketing, mapping = aes(x = facebook, y = sales)) +
  geom_point(color = "orange", size = 2) +
  labs(title  =  'sales ~ facebook', x  =  'facebook') +
  geom_smooth(method = "lm", se = FALSE, color = "black") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

g3 <- ggplot(data = marketing, mapping = aes(x = newspaper, y = sales)) +
  geom_point(color = "forestgreen", size = 2) +
  labs(title  =  'sales ~ newspaper', x  =  'newspaper') +
  geom_smooth(method = "lm", se = FALSE, color = "black") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) 

g1+g2+g3

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Multicolinealidad

data_independiente <- data.frame(marketing$youtube,marketing$facebook,marketing$newspaper)
cor(data_independiente,method = "pearson")

##                     marketing.youtube marketing.facebook marketing.newspaper
## marketing.youtube          1.00000000         0.05480866          0.05664787
## marketing.facebook         0.05480866         1.00000000          0.35410375
## marketing.newspaper        0.05664787         0.35410375          1.00000000

Desarrollo del modelo

modelo_2 <- lm(formula = sales ~ ., data = marketing)

summary(modelo_2)

## 
## Call:
## lm(formula = sales ~ ., data = marketing)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.5932  -1.0690   0.2902   1.4272   3.3951 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.526667   0.374290   9.422   <2e-16 ***
## youtube      0.045765   0.001395  32.809   <2e-16 ***
## facebook     0.188530   0.008611  21.893   <2e-16 ***
## newspaper   -0.001037   0.005871  -0.177     0.86    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.023 on 196 degrees of freedom
## Multiple R-squared:  0.8972, Adjusted R-squared:  0.8956 
## F-statistic: 570.3 on 3 and 196 DF,  p-value: < 2.2e-16

Vemos que la variable newspaper no es estadisticamente significativa, por lo tanto podemos evalaur el modelosin dicha variable

modelo_2_final <- lm(formula = sales ~ youtube + facebook, data = marketing)

summary(modelo_2_final)

## 
## Call:
## lm(formula = sales ~ youtube + facebook, data = marketing)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.5572  -1.0502   0.2906   1.4049   3.3994 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.50532    0.35339   9.919   <2e-16 ***
## youtube      0.04575    0.00139  32.909   <2e-16 ***
## facebook     0.18799    0.00804  23.382   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.018 on 197 degrees of freedom
## Multiple R-squared:  0.8972, Adjusted R-squared:  0.8962 
## F-statistic: 859.6 on 2 and 197 DF,  p-value: < 2.2e-16

Valores residuales se distribuyen normalmente

# Get the model residuals
model_residuals_2 = modelo_2_final$residuals

# Plot the residuals
qqnorm(model_residuals_2)
# Plot the Q-Q line
qqline(model_residuals_2)

Cuando observamos el gráfico QQ, podemos decir que los residuos no se distribuyen normalmente. Los residuos no se distribuyen normalmente. También vemos en el gráfico de ubicación extendida que la variabilidad (varianzas) de los puntos residuales aumenta con el valor de la variable de resultado ajustada, lo que sugiere variaciones no constantes en los errores residuales (o heteroscedasticidad). Según la prueba de normalidad de Shapiro-Wilk, el valor de p fue inferior a 0,05, por lo que se confirmó.

El test de Shapiro-Wilks plantea la hipótesis nula que una muestra proviene de una distribución normal. Eligimos un nivel de significanza, por ejemplo 0,05, y tenemos una hipótesis alternativa que sostiene que la distribución no es normal.

Tenemos:

H0: La distribución es normal

H1: La distribución no es normal

shapiro.test(model_residuals_2)

## 
##  Shapiro-Wilk normality test
## 
## data:  model_residuals_2
## W = 0.91804, p-value = 4.19e-09

Vemos que el valor de probabilidad (p) es muy inferior a nuestro nivel elegido (0,05), por lo que rechazamos la hipótesis nula.

Trasformacion en la variable

Una posible solución para reducir el problema de heteroscedasticidad es usar una transformación logarítmica o de raíz cuadrada de la variable de resultado. Optamos por la raíz cuadrada

modelo_2_completo <- lm(formula = sales ~ sqrt(youtube) + facebook, data = marketing)

summary(modelo_2_completo)

## 
## Call:
## lm(formula = sales ~ sqrt(youtube) + facebook, data = marketing)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3597 -1.0217  0.0445  1.0318  4.0953 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -1.941517   0.390781  -4.968 1.46e-06 ***
## sqrt(youtube)  1.067899   0.026249  40.683  < 2e-16 ***
## facebook       0.194496   0.006677  29.131  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.677 on 197 degrees of freedom
## Multiple R-squared:  0.929,  Adjusted R-squared:  0.9282 
## F-statistic:  1288 on 2 and 197 DF,  p-value: < 2.2e-16

Podemos ver una mejora en el F-statistic y R-squared

# Get the model residuals
model_residuals_completo = modelo_2_completo$residuals

# Plot the residuals
qqnorm(model_residuals_completo)
# Plot the Q-Q line
qqline(model_residuals_completo)

shapiro.test(model_residuals_completo)

## 
##  Shapiro-Wilk normality test
## 
## data:  model_residuals_completo
## W = 0.99076, p-value = 0.23

Vemos que el valor de probabilidad (p) es superios a nuestro nivel elegido (0,05), por lo que no rechazamos la hipótesis nula.

Homocedasticidad

par(mfrow = c(2, 2))
plot(modelo_2_completo)

Test de homocedasticidad: La función bptest() en R es un test de Breusch-Pagan para la heterocedasticidad en modelos de regresión. Esta función toma como entrada un modelo de regresión y devuelve el resultado de la prueba de hipótesis para la homocedasticidad de los residuos.

HO: los residuos tienen varianza constante (homocedasticidad)

HA: hay heterocedasticidad en los residuos

El resultado de la prueba incluye el valor del estadístico de prueba (el valor de la prueba de Breusch-Pagan), el valor-p y el número de grados de libertad. Si el valor-p es menor que el nivel de significancia elegido, se rechaza la hipótesis nula de homocedasticidad y se concluye que hay heterocedasticidad en los residuos.

library(lmtest)
bptest(modelo_2_completo)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_2_completo
## BP = 7.3121, df = 2, p-value = 0.02583

Vemos que el valor de probabilidad (p) es muy superios a nuestro nivel elegido (0,05), por lo que no rechazamos la hipótesis nula.

Modelo final

Vemos que la varianza de la homogeneidad todavía está deteriorada. Eliminemos la observación 131, que es el valor atípico.

marketing_new = marketing[-c(131),]

modelo_sales <- lm(formula = sales ~ sqrt(youtube) + facebook, data = marketing_new)

summary(modelo_sales)

## 
## Call:
## lm(formula = sales ~ sqrt(youtube) + facebook, data = marketing_new)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5602 -1.0379  0.0311  0.9549  3.9930 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -1.731194   0.380122  -4.554 9.22e-06 ***
## sqrt(youtube)  1.048788   0.025733  40.757  < 2e-16 ***
## facebook       0.196705   0.006456  30.467  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.616 on 196 degrees of freedom
## Multiple R-squared:  0.9324, Adjusted R-squared:  0.9318 
## F-statistic:  1353 on 2 and 196 DF,  p-value: < 2.2e-16

par(mfrow = c(2, 2))
plot(modelo_sales)

\(\beta_0 = -1.731194\) \(\beta_1 = 1.048788\) \(\beta_2 = 0.196705\)

Por lo tanto el modelo nos quedaría de la siguiente forma

\(y=1.731194+1.048788x_1+0.196705x_2\)

Prediccion

¿Para un presupuesto publicitario en youtube de 250 y facebook de 100 (miles de dolares) cual seria el valor de nuestras ventas?

nuevo.dato <- data.frame(youtube = 250, facebook = 100)

prediccion <- predict(modelo_sales, newdata = nuevo.dato)

paste("La cantidad estimada de tiempo de coccion es:", round(prediccion, 2))

## [1] "La cantidad estimada de tiempo de coccion es: 34.52"

Por lo tanto, dado el presupuesto publicitario en youtube de 250 y facebok de 100 ( miles de dolares) podemos esperar un valor en las ventas de 34.52 miles de dolares

Resumen del modelo 2

el modelo de regresión lineal múltiple en cuestión tiene como variable dependiente “sales” y dos variables independientes “sqrt(youtube)” y “facebook”. El modelo se ajustó a los datos de la variable dependiente utilizando la función lm() en R y se ajustó a los datos de la variable independiente utilizando la fórmula \("sales ~ sqrt(youtube) + facebook"\). El conjunto de datos utilizado para ajustar el modelo se llama “marketing_new”.

El resumen proporcionado por el modelo incluye varios resultados importantes. A continuación, se describen algunos de los más importantes:

Coeficientes: La tabla de coeficientes muestra el valor estimado de cada coeficiente para cada variable independiente. El coeficiente de “facebook” es 0.196705 y el de “sqrt(youtube)” es 1.048788. El coeficiente de intercepción es -1.731194.

Estadística t: La estadística “t” es una medida de la significancia de cada coeficiente. Un valor t grande indica que el coeficiente es significativo en la relación con la variable dependiente. El valor “t” para la intercepción es -4.554, lo que indica que el valor estimado es significativamente diferente de cero. Los valores “t” para “facebook” y “sqrt(youtube)” son significativamente diferentes de cero con valores de 30.467 y 40.757, respectivamente.

Significancia: Los valores “p” son una medida de la probabilidad de que el valor del coeficiente sea cero. Un valor de p pequeño (menos de 0.05) indica que el coeficiente es significativo. En este caso, todos los valores “p” son significativamente más pequeños que 0.05.

R-cuadrado: El R-cuadrado es una medida de la cantidad de variabilidad en la variable dependiente que se puede explicar por las variables independientes incluidas en el modelo. En este caso, el R-cuadrado es 0.9324, lo que indica que las variables independientes incluidas en el modelo explican el 93.24% de la variabilidad en la variable dependiente.

F-statistic: La estadística F se utiliza para evaluar si el modelo en su conjunto es significativo. Un valor F grande indica que el modelo es significativo. El valor F para este modelo es 1353, con un valor p muy pequeño (menos de 0.05), lo que indica que el modelo en su conjunto es significativo.