R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

R Markdown Ejemplo de Regresión Lineal Simple Anuncios y Carros Vendidos

Este es el archivo Markdown del ejemplo de Regresión Lineal Simple en R, de la Base de Datos Anuncios y Carros vendidos Este ejercicio es una actividad del Módulo 5 de la Concentración CD2001C Analítica para Negocios: De los datos a las decisiones del semestre Agosto - Diciembre 2022. Profesora: Dra. Ma. Margarita Orozco Gómez.

Instalación de librerías

Para iniciar la ejecución de este Markdown es necesario instalar en la consola las librerías: tidyverse, readxl, ggplot2, janitor, dplyr, dbplyr, modeest, y ggthemes.

Archivo R Script de instalación de librerías

Para ello, es recomendable ejecutar primero el archivo R Script de instalación de librerías en el que están las siguientes instrucciones: install.packages(‘tidyverse’) install.packages(‘readxl’) install.packages(‘ggplot2’) install.packages(‘janitor’) install.packages(‘dplyr’) install.packages(‘dbplyr’) install.packages(“modeest”) install.packages(“ggthemes”) install.packages(“ggthemes”) install.packages(“psych”) install.packages(“lubridate”)

Llamado de las librerías en Markdown

Después de ejecutar las librerías en el R Script o la consola, en el Markdown hay que llamarlas para que estén habilitadas las funciones que utilizaremos.

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6      ✔ purrr   0.3.5 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.4.1 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(readxl)
library(ggplot2)
library(dplyr)
library(dbplyr)
## 
## Attaching package: 'dbplyr'
## 
## The following objects are masked from 'package:dplyr':
## 
##     ident, sql
library(janitor)
## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test
library(modeest)
## Registered S3 method overwritten by 'rmutil':
##   method         from
##   print.response httr
library(ggthemes)
library(psych)
## Registered S3 method overwritten by 'psych':
##   method         from  
##   plot.residuals rmutil
## 
## Attaching package: 'psych'
## 
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
library(lubridate)
## 
## Attaching package: 'lubridate'
## 
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

Abrir el archivo de datos

library(readxl)
carros1 <- read_excel("Datos Anuncios y Carros V1.xlsx")

#View(carros1)

Cambiar nombres de variables a minúsculas y quitar acentos

Esto hará que escribir los nombres de las variables en los análisis posteriores sea más fácil, ya que no tendremos que recordar qué es mayúscula y qué no y si le habíamos o no puesto acentos. Instala el paquete janitor:

library(janitor)
carros2 <- carros1 %>% 
  janitor::clean_names()
names(carros2)
## [1] "mes"             "anuncios"        "carros_vendidos"

Escala de medición de las variables

Para analizar la escala de medición de las variables utilizaremos la función sapply(data, class)

sapply(carros2, class)
##             mes        anuncios carros_vendidos 
##       "numeric"       "numeric"       "numeric"

Otra función que nos sirve para analizar la escala de medición de las variables es str(data)

str(carros2)
## tibble [14 × 3] (S3: tbl_df/tbl/data.frame)
##  $ mes            : num [1:14] 1 2 3 4 5 6 7 8 9 10 ...
##  $ anuncios       : num [1:14] 20 36 36 50 200 195 140 126 195 158 ...
##  $ carros_vendidos: num [1:14] 7 4 8 7 18 19 9 9 13 17 ...

Recodificando variables cadena en factores

Las variables que tienen definidas categorías Ej. nominales, ordinales, o de intervalo, en R deben ser codificadas como factores. Para ello, utilizaremos la función as.data.frame en combinación con la función unclass para convertir todas las columnas que aparecen como cadena (characters/chr) en variables tipo factor.

carros3 <- as.data.frame (unclass(carros2),            
                        stringsAsFactors = TRUE)
sapply(carros3, class)
##             mes        anuncios carros_vendidos 
##       "numeric"       "numeric"       "numeric"

Quitar notación científica

Para que las escalas no aparezcan en notación científica en los gráficos

options(scipen=999)

Histograma de la variable dependiente Carros Vendidos

hist(carros3$carros_vendidos, density = 20, prob = TRUE, 
     xlab = "Carros vendidos",
     main = "Distribución de la Variable Dependiente")
curve(dnorm(x, mean = mean(carros3$carros_vendidos), sd = sd(carros3$carros_vendidos)),
      col = "red", lwd = 2, add = TRUE,
      yaxt = "n")

### Gráfico QQ de la variable dependiente Número de Carros Vendidos

# QQ Plot
qqnorm (carros3$carros_vendidos, ylab = "Carros Vendidos", col = "red")
qqline(carros3$carros_vendidos)

Gráfico de caja y bigotes de la variable dependiente Número de Carros Vendidos

boxplot(carros3$carros_vendidos, col = "green",
  xlab= " ",
    ylab= "Número de carros vendidos",
  main = "Gráfico de la Variable Dependiente")

### Histograma de la variable independiente Número de Anuncios

hist(carros3$anuncios, density = 20, prob = TRUE, 
     xlab = "Anuncios",
     main = "Histograma de la variable independiente")
curve(dnorm(x, mean = mean(carros3$anuncios), sd = sd(carros3$anuncios)),
      col = "red", lwd = 2, add = TRUE,
      yaxt = "n")

Gráfico QQ de la variable independiente Número de anuncios

# QQ Plot
qqnorm (carros3$anuncios, ylab = "Anuncios", col = "red")
qqline(carros3$anuncios)

Gráfico de caja y bigotes de la variable independiente Número de anuncios

boxplot(carros3$anuncios, col = "green",
  xlab= " ",
    ylab= "Número de anuncios",
  main = "Gráfico de la variable independiente")

### Identificación de datos atípicos con Z-score en la variable dependiente

carros3$carros_z <- (carros3$carros_vendidos - mean(carros3$carros_vendidos, na.rm=T))/sd(carros3$carros_vendidos, na.rm = TRUE)
summary(carros3$carros_z)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1.3933 -0.7773 -0.3667  0.0000  0.9166  1.6866

El resultado nos índica que no hay observaciones que definidos como OUTLIERS en la variable dependiente. Esto, porque todas las puntuaciones z están dentro de -3 a 3. Nora: tomaremos como referencia el 3 porque tenemos pocos datos.

Identificación de datos atípicos con Z-score en la variable independiente

carros3$anuncios_z <- (carros3$anuncios - mean(carros3$anuncios, na.rm=T))/sd(carros3$anuncios, na.rm = TRUE)
summary(carros3$anuncios_z)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -1.476840 -0.895038 -0.004523  0.000000  0.696015  1.372805

El resultado nos índica que no hay observaciones que definidos como OUTLIERS en la variable independiente. Esto, porque todas las puntuaciones z están dentro de -3 a 3. En caso de haber encontrado datos atípicos, se recomienda winsorizar la variable.

Gráfico de dispersión XY

ggplot(carros3, aes(anuncios, carros_vendidos)) +
    geom_point()+
  labs (x = "Número de anuncios", y = "Número de carros vendidos")+
    geom_smooth(method ="lm", color ="red", se = TRUE)
## `geom_smooth()` using formula 'y ~ x'

### Cálculo de la correlación entre X y Y

cor(carros3$anuncios, carros3$carros_vendidos)
## [1] 0.8588718

Modelo de Regresión Lineal Simple

model1<-lm(carros_vendidos~ anuncios, data = carros3) #Create a model 
summary(model1)
## 
## Call:
## lm(formula = carros_vendidos ~ anuncios, data = carros3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5548 -2.3877  0.3159  2.3770  3.2532 
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  3.28371    1.46610   2.240    0.0448 *  
## anuncios     0.06622    0.01140   5.809 0.0000836 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.596 on 12 degrees of freedom
## Multiple R-squared:  0.7377, Adjusted R-squared:  0.7158 
## F-statistic: 33.74 on 1 and 12 DF,  p-value: 0.00008358

Interpretación de los resultados del Modelo de Regresión Lineal Simple

Prueba t para el Intercepto (Beta0)

Hipótesis nula es que Beta0 = 0 “El intercepto no es significativo” Hipótesis alterna es que Beta0 =! 0 “El intercepto sí es significativo”

Con una t de 2.240 que tiene una significancia de 0.0448 menor que alfa de 0.05, se rechaza HO con una confianza de 95% (*). Es decir, el Intercepto sí es significativo en la ecuación de regresión.

Prueba t para la pendiente (Beta1)

Hipótesis nula es que Beta1 = 0 “La pendiente poblacional es igual a cero. Es decir, el número de anuncios no sirve (es independiente) para pronosticar el número de carros vendidos” Hipótesis alterna es que Beta1 =! 0 “La pendiente poblacional es diferente de cero. Es decir, el número de anuncios sí ayuda a pronosticar el número de carros vendidos”

Con una t de 5.809 que tiene una significancia de 0.0000836 menor que alfa de 0.001, se rechaza HO con una confianza de 99.9% (***). Es decir, la pendiente sí es significativa en la ecuación de regresión. Se valida estadísticamente que el coeficiente Beta (pendiente poblacional) es diferente de 0.

Prueba F para el Modelo de Regresión

Con una F de 33.74 que tiene una significancia de 0.00008358 menor que alfa de 0.001, se rechaza HO con una confianza de 99.9% (***). Es decir, la pendiente sí es significativa en la ecuación de regresión. En otras palabras, el número de carros vendidos depende del número de anuncios publicados.

Interpretación de los coeficientes de la ecuación de regresión.

Interpretación del intercepto a

Independientemente de la cantidad de anuncios que se publiquen, se espera que como mínimo se vendan 3.28 autos al mes.

Interpretación de la pendiente b

Como b = 0.06622, Por cada anuncio que se publique, se espera que las ventas mensuales se incrementen en 0.06622 autos.

Error Cuadrático Medio (MSE)

mse <- mean(residuals(model1)^2)
rmse <- sqrt(mse)
mse #Mean Square Error (ECM: Error Cuadrático Medio)
## [1] 5.778156
rmse #Root Mean Square Error
## [1] 2.40378
y.pron<-predict(model1) #Valores pronosticados
res<- resid(model1) #Residul (Y - y.pron)
stand.res<-rstandard(model1) #Residuos estandarizados 
stud.res<-rstudent(model1) #Residuos estudentizados ("Leave one out" residual)

hist(stand.res, density=20, prob=TRUE, 
     xlab="Residuos estandarizados", 
     main="Gráfico de los Residuos estandarizados")

curve(dnorm(x, mean=mean(stand.res), sd=sd(stand.res)), 
      col="darkblue", lwd=2, add=TRUE, yaxt="n")

qqnorm(stand.res, ylab = "Residuos estandarizados", col = "navyblue")
qqline(stand.res)

Prueba de Normalidad de los residuos

Recordando que, uno de los supuestos más importantes en el Análisis de Regresión, es que, los residuos (errores) tengan una distribución Normal. Es decir, que se distribuyan arriba y debajo del 0 entre -2 y 2 (los estandarizados), hacemos la prueba de Normalidad.

Hipotesis nula: Los residuos se distribuyen Normal Hipótesis alterna: Los residuos no se distribuyen Normal

Con un estadístico de Kolmogorov-Sminorv de 0.1822 y una significancia (p.value) de 0.2356 mayor que alfa de 0.10, no se rechaza H0 con una confianza de 90%. Sí se cumple el supuesto de la distribución Normal de los residuos. Cabe mencionar que esta prueba la podemos hacer tanto para los residuos como para los residuos estandarizados y obtendremos lo mismo. Hay que recordar que el estandarizar una variable no la hace Normal, solo cambia sus unidades originales a unidades de desviaciones estándar. Así que, si los residuos son Normales, los estandarizados, también lo serán. Y si no, pues, no.

library(nortest)

sw<-shapiro.test(res) #Shapiro-Wilk
ks<-lillie.test(res) #kolmogorov-Smirnov
cv<-cvm.test(res) #Cramer-von Mises
ad<-ad.test(res) #Anderson-Darling

stat<-cbind(sw$statistic,ks$statistic,cv$statistic,ad$statistic)
p.value<-cbind(sw$p.value,ks$p.value,cv$p.value,ad$p.value)

normal.table<-rbind(stat,p.value)
colnames(normal.table)<-c("Shapiro-Wilk","kolmogorov-Smirnov",
                          "Cramer-von Mises","Anderson-Darling")
rownames(normal.table)<-c("statistic", "p-value")

normal.table
##           Shapiro-Wilk kolmogorov-Smirnov Cramer-von Mises Anderson-Darling
## statistic   0.89235825          0.1822405        0.0889825        0.5773946
## p-value     0.08741804          0.2356481        0.1443133        0.1096960

Estimación puntual y de intervalo para Y, dado un valor de X

Estimación puntual

Si utilizamos la ecuación de regresión para pronosticar el número de carros que se venderían si publicáramos 50 anuncios (y.pron = a+bX = 3.2837 + 0.06622(50)), el resultado nos indica que esperaríamos vender 6.59 autos. Sin embargo, sabemos que la probabilidad de que este valor sea correcto o exacto es muy poca, tirando a 0. Por ello, es mejor calcular el intervalo de confianza para la estimación.

Estimación por intervalo para Y dado un valor de X

Si se publican 50 anuncios en un mes, se espera con un 90% de confianza que al menos se vendan 1.6 autos, y como máximo, 11.6.

nuevo <- data.frame(anuncios=50)
predict(object=model1, newdata=nuevo, interval="prediction", level=0.90)
##        fit      lwr      upr
## 1 6.594811 1.635303 11.55432

Estimación intervalo para la Media de Y

Si se publicaran 50 anuncios en diferentes meses, en promedio la venta de autos de dichos meses se espera que con un 90% de confianza se encuentre entre 4.8 y 8.4 autos.

nuevo <- data.frame(anuncios=50)
predict(object=model1, newdata=nuevo, interval="confidence", level=0.90)
##        fit      lwr      upr
## 1 6.594811 4.810709 8.378912

###Ahora vamos a obtener todos los IC para Y y los vamos a almacenar en el objeto future_y que luego luego vamos a agregar al conjunto de datos original.

future_y <- predict(object=model1, interval="prediction", level=0.95)
## Warning in predict.lm(object = model1, interval = "prediction", level = 0.95): predictions on current data refer to _future_ responses
carros4 <- cbind(carros3, future_y)

Con el código mostrado a continuación se construye el diagrama de dispersión y se agrega la línea de regresión (en azul) y los IC para las Medias de Y pronosticado (en rosado) por medio de geom_smooth. Los IC para el pronóstico de Y (en rojo) se agregan por medio de geom_line.

library(ggplot2)
ggplot(carros4, aes(x=anuncios, y=carros_vendidos))+
    geom_point() +
    geom_line(aes(y=lwr), color="red", linetype="dashed") +
    geom_line(aes(y=upr), color="red", linetype="dashed") +
    geom_smooth(method=lm, formula=y~x, se=TRUE, level=0.95, col='blue', fill='pink2') +
    theme_light()