Regresion Lineal Multiple
Primeros pasos
Para comenzar a analizar nuestros datos debemos cargar la base de datos que vamos a utilizar, asi como las librerias que vamos a utilizar.
library(readr)
library(DataExplorer)
library(corrplot)
library(dplyr)
library(lmtest)
library(zoo)
library(het.test)
library(nortest)
datos <- read_delim("datos.csv", ";", escape_double = FALSE, trim_ws = TRUE)
datos <- data.frame(datos)Demos un pequeno vistazo a los datos.
head(datos)Ahora vamos a ver que tanta informacion tenemos en nuestra base de datos. Para ello nos apoyaremos de una funcion en R.
introduce(datos)Viendo nuestra base de datos notamos que las columnas “like”, “comment” y “share” en nuestro modelo nos representarian un problema, pues estas 3 variables son simplemente la suma de nuestra variable a estudiar, es decir, “Total interactions”. Lo que ocasionaria que nuestro modelo unicamente se fije en esas tres variables para hacer predicciones, por ende terminaria ignorando otras variables.
Por lo que procedemos a quitarlas de nuesto Data Frame.
df <- datos[,c(-18,-17,-16)]Lo siguiente seria revisar si nuestros datos tienen datos faltantes, para ello nos apoyamos de una sencilla grafica.
df_1 <- df
names(df_1) <- seq(1,16)
plot_missing(df_1)
Notamos que los datos estan relativamente completos, unicamente tenemos 1 observaciones que tienen datos incompletos, lo que representa el 0.2% de nuestra base. Por lo que procedemos a eliminarlo.
Adicional a ello vamos a modificar una variable categorica, para poder utilizarla en el modelo.
df <- df %>% na.omit()
df$Type[df$Type == "Link"] <- 1
df$Type[df$Type == "Photo"] <- 2
df$Type[df$Type == "Status"] <- 3
df$Type[df$Type == "Video"] <- 4
df$Type <- as.numeric(df$Type)
df_1 <- df
#Estamos haciendo dos DataFrame para poder visualizar uno sin los nombres de las columnas, para una mejor
#VisualizacionUna vez hecho esto procedemos a hacer un analisis de las correlaciones de las variables que posiblemente esten en nuestro modelo.
names(df_1) <- seq(1,16)
corrplot(cor(df_1[,-16]))
En este caso vemos que las primeras 7 variables no tienen correlaciones significativas, sin embargo para las ultimas 8 existen correlaciones importantes entre todas, por lo que habra que hacerse un analisis de estas para descartar algunas.
Primer modelo
En primera instancia vamos a analizar un modelo que contemple todas las variables para ver como se comporta, asi como, verificar que los supuestos del modelo de regresion se cumplen, de lo contrario no tendriamos un modelo sustentado en la teoria.
modelo1 <- lm(Total.Interactions ~ ., data = df)
summary(modelo1)##
## Call:
## lm(formula = Total.Interactions ~ ., data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -374.96 -20.32 -2.26 19.88 508.78
##
## Coefficients:
## Estimate
## (Intercept) -2.896e+01
## Page.total.likes 7.836e-04
## Type 1.283e+01
## Category -3.207e+00
## Post.Month -9.033e+00
## Post.Weekday -2.597e-01
## Post.Hour -1.103e+00
## Paid -5.170e+00
## Lifetime.Post.Total.Reach -1.428e-03
## Lifetime.Post.Total.Impressions 5.362e-04
## Lifetime.Engaged.Users 1.468e+00
## Lifetime.Post.Consumers -1.452e+00
## Lifetime.Post.Consumptions 4.798e-03
## Lifetime.Post.Impressions.by.people.who.have.liked.your.Page -3.563e-04
## Lifetime.Post.reach.by.people.who.like.your.Page -2.318e-03
## Lifetime.People.who.have.liked.your.Page.and.engaged.with.your.post 4.241e-02
## Std. Error
## (Intercept) 5.518e+01
## Page.total.likes 5.828e-04
## Type 8.842e+00
## Category 4.038e+00
## Post.Month 2.864e+00
## Post.Weekday 1.531e+00
## Post.Hour 7.320e-01
## Paid 6.989e+00
## Lifetime.Post.Total.Reach 5.993e-04
## Lifetime.Post.Total.Impressions 2.135e-04
## Lifetime.Engaged.Users 1.791e-02
## Lifetime.Post.Consumers 1.813e-02
## Lifetime.Post.Consumptions 2.230e-03
## Lifetime.Post.Impressions.by.people.who.have.liked.your.Page 2.365e-04
## Lifetime.Post.reach.by.people.who.like.your.Page 1.263e-03
## Lifetime.People.who.have.liked.your.Page.and.engaged.with.your.post 1.325e-02
## t value
## (Intercept) -0.525
## Page.total.likes 1.344
## Type 1.451
## Category -0.794
## Post.Month -3.154
## Post.Weekday -0.170
## Post.Hour -1.507
## Paid -0.740
## Lifetime.Post.Total.Reach -2.384
## Lifetime.Post.Total.Impressions 2.512
## Lifetime.Engaged.Users 81.935
## Lifetime.Post.Consumers -80.091
## Lifetime.Post.Consumptions 2.151
## Lifetime.Post.Impressions.by.people.who.have.liked.your.Page -1.507
## Lifetime.Post.reach.by.people.who.like.your.Page -1.835
## Lifetime.People.who.have.liked.your.Page.and.engaged.with.your.post 3.200
## Pr(>|t|)
## (Intercept) 0.59999
## Page.total.likes 0.17942
## Type 0.14735
## Category 0.42753
## Post.Month 0.00171
## Post.Weekday 0.86535
## Post.Hour 0.13243
## Paid 0.45985
## Lifetime.Post.Total.Reach 0.01753
## Lifetime.Post.Total.Impressions 0.01233
## Lifetime.Engaged.Users < 2e-16
## Lifetime.Post.Consumers < 2e-16
## Lifetime.Post.Consumptions 0.03194
## Lifetime.Post.Impressions.by.people.who.have.liked.your.Page 0.13258
## Lifetime.Post.reach.by.people.who.like.your.Page 0.06708
## Lifetime.People.who.have.liked.your.Page.and.engaged.with.your.post 0.00146
##
## (Intercept)
## Page.total.likes
## Type
## Category
## Post.Month **
## Post.Weekday
## Post.Hour
## Paid
## Lifetime.Post.Total.Reach *
## Lifetime.Post.Total.Impressions *
## Lifetime.Engaged.Users ***
## Lifetime.Post.Consumers ***
## Lifetime.Post.Consumptions *
## Lifetime.Post.Impressions.by.people.who.have.liked.your.Page
## Lifetime.Post.reach.by.people.who.like.your.Page .
## Lifetime.People.who.have.liked.your.Page.and.engaged.with.your.post **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 68.37 on 483 degrees of freedom
## Multiple R-squared: 0.9687, Adjusted R-squared: 0.9677
## F-statistic: 996.7 on 15 and 483 DF, p-value: < 2.2e-16Vemos que en este modelo, hay variables que no son significativas como por ejemplo la variable Post Weekday y otras como Lifetime Post Consumers, si lo son. Esto nos da una idea de que variables podriamos seleccionar para un mejor modelo. Ahora bien, la estadistica Adjusted R-squared nos da un valor muy alto, lo que nos indica que nuestro modelo cubre alrededor del 96% de la varianza. Aunque esto podria ser o no algo bueno, quizas haya variables que esten muy ligadas al total de interacciones. Ahora revisaremos los supuestos del modelo.
El primer supuesto que vamos a analizar es de la Homocedasticidad de los residuos, es decir, que los residuos tengan una varianza constante, esto entre otras cosas nos indicaria que no hay una relacion entre las diferencias y los valores ajustados.
Primero vamos a analizarlo graficamente y procederemos con su respectiva prueba.
Ajustados<-fitted(modelo1)
Res1<-stdres(modelo1)
par(mfrow=c(1,2))
plot(Res1,ylab= "Residuales")
abline(a=0,b=0)
plot(Ajustados,Res1,ylab= "Residuales",xlab="Ajustados")
bptest(modelo1)##
## studentized Breusch-Pagan test
##
## data: modelo1
## BP = 239.61, df = 15, p-value < 2.2e-16En la primer grafica, estamos graficando los residuales con su respectivo indice en el DataFrame, pareciera que los errores no siguen ningun patron… Pero realizando la prueba bptest con Hipotesis nula de Homocedasticidad tenemos informacion para rechazar la hipotesis nula, por lo que al menos una es diferente.
Vamos a analizar la normalidad de los residuales, es decir, que las diferencias entre los valores ajustados y los reales, sigan una distribucion Normal. Primero vamos a ver la distribucion de los datos en un histograma.
hist(Res1,probability = T,breaks = 70, col = "darkorange",main="Histograma de Residuales",ylim=c(0,1))
lines(density(Res1),lwd=3,col="darkblue")
Pareciera que si tenemos una forma parecida a la de una normal. Pero para estar seguros es necesario hacer mas analisis, ahora vamos a ver como es que se comporta el QQplot de estos datos.
qqnorm(Res1)
qqline(Res1)
De igual forma observamos que en el centro de los datos si se sigue una linea recta, aunque en los extremos los cuantiles se desvian mucho. Para estar seguros de que realmente distribuyen normal vamos a aplicar una prueba de bondad de ajuste, esto nos dira con mayor certeza si realmente nuestros datos se comportan de la manera que esperamos.
ks.test(Res1,"pnorm",0,sd(Res1))##
## One-sample Kolmogorov-Smirnov test
##
## data: Res1
## D = 0.18592, p-value = 2.109e-15
## alternative hypothesis: two-sidedlillie.test(Res1)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: Res1
## D = 0.18872, p-value < 2.2e-16ad.test(Res1)##
## Anderson-Darling normality test
##
## data: Res1
## A = 37.749, p-value < 2.2e-16El resultado de la prueba nos arroja un P-value muy pequeno en todas las pruebas, por lo que podemos decir que hay evidencia estadistica para decir que nuestros errores no siguen una distribucion normal. Por lo que nuestro modelo no cumple este supuesto. Puede deberse a que las desviaciones se concentran en las colas de la distribucion pero el resto de los datos se ajustan de buena manera a la distribucion normal en la grafica sin presentar asimetria para datos centrales.
Segundo Modelo
Debido a que encontramos varias variables que no eran significativas para el modelo y que no se cumplian dos de los principales supuestos para la regresion, procederemos a seleccionar diferentes varaibles para conformar un mejor modelo con variables mas significativas. Consideraremos el grafico de correlacion para apoyarnos.
corrplot(cor(df_1[,-16]))
df_2 <- df[,c(1,2,3,8,9,10,14,16)]
modelo2 <- lm(Total.Interactions ~ ., data = df_2)
summary(modelo2)##
## Call:
## lm(formula = Total.Interactions ~ ., data = df_2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1484.4 -66.8 10.8 76.9 3845.6
##
## Coefficients:
## Estimate Std. Error
## (Intercept) -3.788e+02 1.035e+02
## Page.total.likes 4.137e-03 7.599e-04
## Type -1.847e+02 3.052e+01
## Category 9.416e+01 1.396e+01
## Lifetime.Post.Total.Reach 2.913e-03 8.862e-04
## Lifetime.Post.Total.Impressions -7.205e-04 2.231e-04
## Lifetime.Engaged.Users 1.347e-01 1.676e-02
## Lifetime.Post.reach.by.people.who.like.your.Page 2.185e-02 2.550e-03
## t value Pr(>|t|)
## (Intercept) -3.661 0.000278 ***
## Page.total.likes 5.444 8.26e-08 ***
## Type -6.053 2.83e-09 ***
## Category 6.744 4.35e-11 ***
## Lifetime.Post.Total.Reach 3.288 0.001082 **
## Lifetime.Post.Total.Impressions -3.229 0.001326 **
## Lifetime.Engaged.Users 8.038 6.85e-15 ***
## Lifetime.Post.reach.by.people.who.like.your.Page 8.569 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 259.9 on 491 degrees of freedom
## Multiple R-squared: 0.5403, Adjusted R-squared: 0.5338
## F-statistic: 82.45 on 7 and 491 DF, p-value: < 2.2e-16Consideramos estas variables principalmente porque despues de realizar varias pruebas, encontramos que son las mas significativas para la variable que tratamos de explicar. Notamos que tiene sentido ya que algunas de las variables estaban muy correlacionadas y con el apoyo del nombre de la columna podemos concluir que llevaban la misma informacion. Podemos encontrar que el valor de la adjusted R-squared disminuye pero obtenemos un mayor numero de variables significativas, la diferencia en la R se debe principalmente en anadir la variable numero 11 Lifetime.Post.Consumers la cual tiene una correlacion muy alta con la variable Lifetime.Engaged.Users, para nuestro caso de estudio nos enfocaremos en utilizar Lifetime.Engaged.Users por la forma en la que calcula las interacciones para obtener su informacion como una proporcion entre la suma de los “likes”, “comentarios” y “compartir” entre el numero de seguidores. o Procedemos a revisar los supuestos del modelo, comenzamos por el supuesto de Homocedasticidad.
Ajustados2<-fitted(modelo2)
Res2<-stdres(modelo2)
par(mfrow=c(1,2))
plot(Res2,ylab= "Residuales")
abline(a=0,b=0)
plot(Ajustados2,Res2,ylab= "Residuales",xlab="Ajustados")
bptest(modelo2)##
## studentized Breusch-Pagan test
##
## data: modelo2
## BP = 138.43, df = 7, p-value < 2.2e-16La prueba nos indica que debemos rechazar Homocedasticidad, pero podemos indentificar en las graficas que existen algunos “outliers” los cuales podroan estar afectando directamente la prueba, procederemos a revisar los demas supestos para recolectar mas informacion del modelo.
Vamos a analizar la normalidad de los residuales. Primero vamos a ver la distribucion de los datos en un histograma.
hist(Res2,probability = T,breaks = 70, col = "darkorange",main="Histograma de Residuales",ylim=c(0,1))
lines(density(Res2),lwd=3,col="darkblue")
Procedemos a revisar como se ajustan los datos a una distribucion normal.
qqnorm(Res1)
qqline(Res1)
Nuevamente observamos que en el centro de los datos si se sigue una linea recta, aunque en los extremos los cuantiles se desvian mucho. Por ello podemos esperar que la prueba de normalidad nos indique que rechaza la hipotesis nula pero para este caso debemos considerar la asimetria de los datos centrales en la grafica, seguiremos recolectando informacion del modelo para tener una mejor conclusion. Por ahora podemos considerar que los “outliers” nos estan modificando el comportamiento de los datos con respecto al resto de las observaciones, al tratarse de datos de Facebook podemos considerar que se tratan de post virales los cuales tienen una interaccion fuera del comportamiento natural de la pagina y el alcance que tiene la misma.
ks.test(Res1,"pnorm",0,sd(Res1))##
## One-sample Kolmogorov-Smirnov test
##
## data: Res1
## D = 0.18592, p-value = 2.109e-15
## alternative hypothesis: two-sidedlillie.test(Res1)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: Res1
## D = 0.18872, p-value < 2.2e-16ad.test(Res1)##
## Anderson-Darling normality test
##
## data: Res1
## A = 37.749, p-value < 2.2e-16Confirmamos que las pruebas nos indican que no sigue una distribuci??n normal, pero de manera grafica podemos notar que hay que realizar un analisis mas detallado ya que existe un mejor ajuste para los datos centrales.
Ahora revisaremos la multicolinealidad del modelo para corroborar que elegimos de una mejor forma las variables a participar en el modelo.
X1=scale(df_2)
cor(X1)## Page.total.likes
## Page.total.likes 1.00000000
## Type 0.21280098
## Category -0.09102340
## Lifetime.Post.Total.Reach -0.08603765
## Lifetime.Post.Total.Impressions -0.10474427
## Lifetime.Engaged.Users -0.11456623
## Lifetime.Post.reach.by.people.who.like.your.Page -0.06390869
## Total.Interactions 0.04426378
## Type Category
## Page.total.likes 0.21280098 -0.091023405
## Type 1.00000000 0.057114576
## Category 0.05711458 1.000000000
## Lifetime.Post.Total.Reach 0.07836723 -0.141980847
## Lifetime.Post.Total.Impressions 0.04851593 -0.094295873
## Lifetime.Engaged.Users 0.35059510 0.003494811
## Lifetime.Post.reach.by.people.who.like.your.Page 0.18287762 -0.104330957
## Total.Interactions 0.05019319 0.127386695
## Lifetime.Post.Total.Reach
## Page.total.likes -0.08603765
## Type 0.07836723
## Category -0.14198085
## Lifetime.Post.Total.Reach 1.00000000
## Lifetime.Post.Total.Impressions 0.69486342
## Lifetime.Engaged.Users 0.57049773
## Lifetime.Post.reach.by.people.who.like.your.Page 0.74294243
## Total.Interactions 0.53851846
## Lifetime.Post.Total.Impressions
## Page.total.likes -0.10474427
## Type 0.04851593
## Category -0.09429587
## Lifetime.Post.Total.Reach 0.69486342
## Lifetime.Post.Total.Impressions 1.00000000
## Lifetime.Engaged.Users 0.36842206
## Lifetime.Post.reach.by.people.who.like.your.Page 0.65186891
## Total.Interactions 0.34326509
## Lifetime.Engaged.Users
## Page.total.likes -0.114566227
## Type 0.350595103
## Category 0.003494811
## Lifetime.Post.Total.Reach 0.570497730
## Lifetime.Post.Total.Impressions 0.368422064
## Lifetime.Engaged.Users 1.000000000
## Lifetime.Post.reach.by.people.who.like.your.Page 0.611945768
## Total.Interactions 0.572091190
## Lifetime.Post.reach.by.people.who.like.your.Page
## Page.total.likes -0.06390869
## Type 0.18287762
## Category -0.10433096
## Lifetime.Post.Total.Reach 0.74294243
## Lifetime.Post.Total.Impressions 0.65186891
## Lifetime.Engaged.Users 0.61194577
## Lifetime.Post.reach.by.people.who.like.your.Page 1.00000000
## Total.Interactions 0.61788736
## Total.Interactions
## Page.total.likes 0.04426378
## Type 0.05019319
## Category 0.12738670
## Lifetime.Post.Total.Reach 0.53851846
## Lifetime.Post.Total.Impressions 0.34326509
## Lifetime.Engaged.Users 0.57209119
## Lifetime.Post.reach.by.people.who.like.your.Page 0.61788736
## Total.Interactions 1.00000000A=t(X1)%*%X1
kappa=max(eigen(A)$values)/min(eigen(A)$values)
kappa## [1] 15.04929Al obtener un valor superior a 10 pero cercano al mismo, podemos decir que podriamos tener multicolinealidad, pero ya teniamos conocimiento de esta informacion al revisar el grafico de correlacion ya que existe una alta relacion entre las variables a partir de la variable 8 hasta la 19, es por ello que decidimos quitar las variables que claramente representaban Multicolinealidad para el modelo.
Ahora revisaremos los puntos de apalancamiento del modelo para tener una mayor certeza de las observaciones que afectan directamente al modelo y que podrian estar aportandonos una informacion deficiente del comportamiento de la pagina.
dfbetas(modelo2)plot(dfbetas(modelo2)[,1],type="h", ylab="Beta0")
plot(dfbetas(modelo2)[,2],type="h", ylab="Beta1")
plot(dfbetas(modelo2)[,3],type="h", ylab="Beta2")
plot(dffits(modelo2),type="h")
Encontramos que existen algunas observaciones que estan afectando la informacion que tenemos y por ello la calidad del modelo, para poder saber que accion tomar con estas observaciones seria mejor tener mayor informacion de dichas publicaciones ya que mantienen un alto numero de interacciones y elevan el comportamiento natural de usuario con la pagina de Facebook.
Despues de reunir esta informacion podemos concluir que tenemos un posible buen modelo, pero que se ve afectado principalmente por los “outlier” en las observaciones y que debemos reunir mayor informacion de dichas observaciones para obtener una mejor seleccion de las variables que debemos considerar para mejorar las predicciones de las interacciones totales.