Fecha de la ultima revisión

## [1] "2023-07-10"


Activar las librerias, averiguar que TIENE instalado estos paquetes

library(QuantPsyc)
library(car)
library(boot)
library(ggplot2)
library(tidyverse)
library(sjPlot)

Regresión Lineal

La regresión lineal es el modelo básico de evaluar si hay una relación lineal o sea recta entre dos variable. Esta relación entre las variables puedes ser positiva o negativa. Hay otro tipos de regresión, que incluye regresión no lineal tal como cuadrática \[x^2\] o cúbica \[x^3\], logarítmica \[log(x)\] entre muchas otras alternativas.

Aquí estaremos evaluando solamente la regresión lineal

Primero paso es subir los datos:

El números de barres en un bario y la cantidad de asesinados en este barrio (Datos ficticios)

library(readr)
pubs <- read_csv("Data_files_csv/pubs.csv")
pubs
pubsmortality
101e+03
202e+03
303e+03
404e+03
505e+03
606e+03
707e+03
5001e+04
# Cambiamos los nombres de las columnas

pubs=pubs%>%
  dplyr::rename(barres=pubs, mortalidad=mortality)

la función lm()

Como hacer un regresión lineal (Linear Regression) simple, usando la función lm(), para “linear model”= modelo lineal. Una regresión necesita dos variable continuas (vea el módulo de correlación). Es importante que estas variable tenga una distribución normales. La diferencia entre una correlación y una regresión es que la primera es un análisis que describe el patrón y en la segunda es que no solamente describe el patrón pero hace predicción sobre la relación entre las variables. Usando la regresión uno calcula tambien una linea que describe la relación entre las variables. Esta variable se puede describe como \(y=m_x+b\) donde la m representa la pendiente y la b representa el intercepto. También lo pueden ver en libros de la siguiente forma \(y=\alpha+\beta_x\) donde la \(\beta\) beta representa la pendiente y la \(\alpha\) el intercepto.

La formula de lm() se compone de lm(y~x, data= “df”). Nota la tilde ~. Hay dos pruebas la priemra para determinar si \(\alpha\) es distinto de cero. La hipótesis nula es

Ho: el intercepto \(\alpha\) es igual a cero Ha: el intercepto, \(\alpha\) no es igual a cero. Entonces el punto donde la linea intercepta el cero puede estar mayor de o menor de cero.

La segunda hipótesis nula es que la pendiente es diferente de cero, esto quiere decir que la pendiente no sugiere un patrón de aumentar y disminuir entre las dos variables.

Ho: la pendiente \(\beta\) es igual a cero Ha: la pendiente, \(\alpha\) no es igual a cero. Entonces la relación entre las dos variables es o positiva o negativa.

Ahora evaluamos los resultados de la regresión entre el número de barres en un vecindarios y la mortalidad en este mismo sector. Se observa que los coeficientes de la linea son \(y=3352+14.3*x\). Entonce el intercepto en cero comienza en 3352 fatalidades y por cada bar suplementario hay 14.3 más fatalidad.

Ahora para determinar si estos valores son significativo hay que evaluar el valor de p en cada linea. La hipótesis nula del intercepto tiene un valor de p =0.005, que sugiere que se debería rechazar la hipótesis nula, y por consecuencia aceptamos la hipótesis alterna. La pendiente tiene un valor de p=0.015 y también se rechaza la hipotesis nula.

pubReg <- lm(mortalidad~barres, data = pubs) 

#summary(pubReg) # Si no se acuerda de la funcíon que sigue

tab_model(
  pubReg,show.df = TRUE,
  CSS = list(
    css.depvarhead = 'color: red;',
    css.centeralign = 'text-align: left;', 
    css.firsttablecol = 'font-weight: bold;', 
    css.summary = 'color: blue;'
  )
)
  mortalidad
Predictors Estimates CI p df
(Intercept) 3351.96 1440.34 – 5263.57 0.005 6.00
barres 14.34 3.82 – 24.86 0.016 6.00
Observations 8
R2 / R2 adjusted 0.649 / 0.591

Visualización de la regresión

Se observa que hay un aumento en fatalidades con aumento en el número de barres. Pero nota el valor a la derecha que parece ser muy atípico comparado a los otros.

ggplot(pubs, aes(x=barres, y=mortalidad))+
  geom_smooth(method = lm)+ # modelo lineal
  geom_point()


El impacto de valores atípico

En cierta ocasiones valores fuera de lo normal pueden hacer grandes cambios en el resultado, en este caso la regresión. ¿Cual es el efecto del valor grande? Removemos ese valor del archivo de datos y revaluamos el modelo (la regresión lineal). Nota que ahora el modelo es sumamente diferente \(y=-163.7+103.2*x\). ¿Ahora se rechaza las dos hipótesis?

pubsnew <- pubs[ which(pubs$barres<80), ]  # remover el valor grande

pubsnew=pubsnew %>% 
  add_row(barres = 4, mortalidad = 0) # Añadiendo un par de valores


pubRegNew <- lm(mortalidad~barres, data = pubsnew)
#summary(pubRegNew)

tab_model(
  pubRegNew,show.df = TRUE,
  CSS = list(
    css.depvarhead = 'color: red;',
    css.centeralign = 'text-align: left;', 
    css.firsttablecol = 'font-weight: bold;', 
    css.summary = 'color: blue;'
  )
)
  mortalidad
Predictors Estimates CI p df
(Intercept) -163.70 -374.14 – 46.73 0.106 6.00
barres 103.20 98.18 – 108.23 <0.001 6.00
Observations 8
R2 / R2 adjusted 0.998 / 0.997
ggplot(pubsnew,aes(x=barres, y=mortalidad))+
  geom_smooth(method=lm)+
  geom_point()


La venta de discos de música

Ahora un set de datos más complejo y más realista a los que uno encontraría en un estudio de en medicina, sociología o ecológica.

Los datos representa la cantidad de dinero dedicado a la promoción de diferentes CD’s de una compañia de música y la cantidad de CD (CD/downloads) vendido. El la primera linea se observa la cantidad de libras esterlina, £ (UK) dedicado a la promoción del album de música, en la primera linea vemos que se gasto £10,256, y después la cantidad de CD vendido fue 330.

library(readr)
Album_Sales_1_new <- read_csv("Data_files_csv/Album_Sales_1_new.csv")
head(Album_Sales_1_new)
advertssales
10.3     330
986       120
1.45e+03360
1.19e+03270
575       220
569       170
length(Album_Sales_1_new$adverts)  # cuantas files de datos hay en el data frame
## [1] 200

Comenzamos con hacer un gráfico entre las dos variables. Nota que en la parte de geom_smooth(), tiene que incluir method=lm, esto significa que el método de construir la linea usara la regresión lineal. Se añade a la función lineal \(\epsilon\) que representa los errores de los valores al comparar con la linea que representa el mejor modelo.

\[Y_{ i }=\beta _{ 0 }+\beta _{ 1 }X_{ i }+\epsilon _{ i }\] Recuerda que \(\beta _{ 0 }\) es el intercepto y el \(\beta _{ 1 }X_{ i }\) es la pendiente. El área sombreada es el área de 95 de intervalo de confianza. Esto quiere decir que la mejor linea, intercepto y pendiente podría variar en este rango si repetimos el experimento. Nota aquí todas las alternativas, añadí los dos extremos, con una pendiente mayor (roja) y una menor (violeta). Cada punto representa la venta de CD con su correspondiente cantidad dedicada a la promoción. Los \(epsilon\) seria la diferencia entre la mejor linea y el valor original, esto se llama tambien los residuales.

library(ggplot2)
ggplot(Album_Sales_1_new,aes(x=adverts, y=sales))+
  geom_smooth(method=lm, se = TRUE)+
  geom_point()+
  geom_segment(aes(x=0, y=120, xend=2250, yend=380), colour="red")+
  geom_segment(aes(x=0, y=150, xend=2250, yend=320), colour="purple")

El modelo lineal con la función lm(). ¿Como interpretas los coeficientes y si estos son significativos, si se rechaza la hipótesis nula?

library(sjPlot)
head(Album_Sales_1_new)
advertssales
10.3     330
986       120
1.45e+03360
1.19e+03270
575       220
569       170
model1=lm(sales~adverts, Album_Sales_1_new)

#summary(model1)

tab_model(
  model1,show.df = TRUE,
  CSS = list(
    css.depvarhead = 'color: red;',
    css.centeralign = 'text-align: left;', 
    css.firsttablecol = 'font-weight: bold;', 
    css.summary = 'color: blue;'
  )
)
  sales
Predictors Estimates CI p df
(Intercept) 134.14 119.28 – 149.00 <0.001 198.00
adverts 0.10 0.08 – 0.12 <0.001 198.00
Observations 200
R2 / R2 adjusted 0.335 / 0.331

Supuestos de la regresión lineal

  1. Igualdad de varianza: En el primer gráfico evaluar si los datos están distribuido más o menos igual. y que no un más variación un un área del gráfico que otro.

  2. Normalidad de los datos: En el segundo gráfico evaluar la figura de qqplot.

  3. Evaluar si hay datos sesgados (atípicos) que influencia los resultados

  4. En el cuarto gráfico evaluar si hay valores que tienen mucho peso si se incluyen o no en el análisis. Estos van a ser identificado

En el presente gráfico hay tres valores que hay que evaluar (1,42,169), estos valores se tiene que asegurar que son correctos. Siempre es bueno remover los valores sesgados y rehacer el análisis para observar cuan diferente son los resultas.

plot(model1)  # Evaluar los supuestos, 1. Igualdad de varianza, 2. Normalidad, 4. Datos sesgados (Cook's Distance)


Cook’s Distance

Evaluar si hay unos valores que sesgan los resultados.

FROM WIKIPEDIA Definition. Data points with large residuals (outliers) and/or high leverage may distort the outcome and accuracy of a regression. Cook’s distance measures the effect of deleting a given observation. Points with a large Cook’s distance are considered to merit closer examination in the analysis.

Detecting highly influential observations There are different opinions regarding what cut-off values to use for spotting highly influential points. A simple operational guideline of D_i>1 has been suggested.

Others have indicated that D_i>4/n, where n is the number of observations, might be used.

Aqui enseño como añadir 1. los valores de “cook.distance” a su archivo 2. Añadir una columna de “secuencia” de los datos 3. Crear una gráfica de las distancia de Cook.
4. Which values have Di larger than 1, which values have Cook’s distances larger 4/n?

4/length(Album_Sales_1_new$adverts)
## [1] 0.02
Album_Sales_1_new$cooks.distance<-cooks.distance(model1)
head(Album_Sales_1_new)
advertssalescooks.distance
10.3     3300.0572  
986       1200.011   
1.45e+033600.0178  
1.19e+032700.000662
575       2200.000548
569       1700.000207
Album_Sales_1_new$sequence=c(1:200)
head(Album_Sales_1_new)
advertssalescooks.distancesequence
10.3     3300.0572  1
986       1200.011   2
1.45e+033600.0178  3
1.19e+032700.0006624
575       2200.0005485
569       1700.0002076
tail(Album_Sales_1_new)
advertssalescooks.distancesequence
701      2500.00141 195
911      1900.000804196
889      2400.000321197
801      2500.00101 198
1.5e+032300.00608 199
786      1100.00649 200
ggplot(Album_Sales_1_new, aes(sequence, cooks.distance))+
  geom_point()+
  geom_hline(aes(yintercept=4/length(Album_Sales_1_new$adverts), colour="red"))


Repaso de los PASOS

Paso 1

  1. Primer paso, mirar los coeficientes.
  • El resultado: El coeficiente (intercepto) y la pendiente del model1

Paso 2

  1. Evaluar si los coeficientes son diferentes de cero.
  1. La primera hipótesis, Determinar si el intercepto es igual a cero. Mira el valor de p, Pr(>|t|), determinar si el valor es menor de 0.05, si lo es se rechaza la Ho y por consecuencia tenemos confianza que el intercepto no incluye cero.

  2. La segunda hipótesis nula Deteminar si la pendiente es igual a cero. Mira el valor de p, Pr(>|t|), como el valor es menor de p=0.05, se rechaza la Ho y por consecuencia tenemos confianza que la pendiente no incluye cero.

Paso 3

  1. Evaluar si los datos cumple con los supuestos.

  2. Igualdad de varianza, usa la gráfica de residuales

  3. Normalidad, qqplot

  4. Valores sesgados, la prueba de Cook’s

advertizingReg <- lm(sales~adverts, data = Album_Sales_1_new)
summary(advertizingReg)
## 
## Call:
## lm(formula = sales ~ adverts, data = Album_Sales_1_new)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -152.949  -43.796   -0.393   37.040  211.866 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.341e+02  7.537e+00  17.799   <2e-16 ***
## adverts     9.612e-02  9.632e-03   9.979   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 65.99 on 198 degrees of freedom
## Multiple R-squared:  0.3346, Adjusted R-squared:  0.3313 
## F-statistic: 99.59 on 1 and 198 DF,  p-value: < 2.2e-16
tab_model(
  advertizingReg,show.df = TRUE,
  CSS = list(
    css.depvarhead = 'color: red;',
    css.centeralign = 'text-align: left;', 
    css.firsttablecol = 'font-weight: bold;', 
    css.summary = 'color: blue;'
  )
)
  sales
Predictors Estimates CI p df
(Intercept) 134.14 119.28 – 149.00 <0.001 198.00
adverts 0.10 0.08 – 0.12 <0.001 198.00
Observations 200
R2 / R2 adjusted 0.335 / 0.331

Alternativa para Gráficar los residuales

Plot the residuals and the fitted values

Evaluating residuals:

If we observe that the residuals are randomly distributed (no pattern, a straigth line), then we can assume that there is equality of variance along the “x” axis.

Definition of Residuals

A quantity remaining after other things have been subtracted or allowed for. FROM Wikipedia

A residual (or fitting deviation), on the other hand, is an observable estimate of the unobservable statistical error.

# "advertizingReg", nota que este no es un data frame pero un modelo
# The main figure
ggplot(Album_Sales_1_new,aes(x=adverts, y=sales))+
  geom_smooth(method=lm, se = TRUE)+
  geom_point()

# Graficando los residuales, metodo 1
ggplot(model1, aes(x=.fitted, y=.resid))+
  geom_point()+
  geom_smooth(method=lm)

Porque podemos usar la alternativas anterior

Cuando R corre un analisis, hace muchos calculos y guarda parate de estos calculos como lista de datos.

Para ver algunos de los valores calculado uno pone el nombre del modelo y la añade el $ y deberia aparecer la lista de los valores calculados. o otra alternativa es usar la función ls() = list() para ver el nombre de las variables

ls(model1)
##  [1] "assign"        "call"          "coefficients"  "df.residual"  
##  [5] "effects"       "fitted.values" "model"         "qr"           
##  [9] "rank"          "residuals"     "terms"         "xlevels"
# model1$residuals  # Por ejemplo para ver los residuales del modelo

Ejercicio #1

Usa los enviado por MSTeam llamado “Alturas_Humanos

  • En este archivo vemos tres columnas. Genero, Altura de la persona (Altura_cm) y Peso de la persona (Peso_kg)

  • Evalua cada paso para determinar si hay una relación entre la altura y el peso de las personas.

  • Haz un gráfico de la relación entre la altura y el peso

  • Evalúa el modelo de regresión entre las dos variables

  • Evalúa los supuestos

library(readr)
Alturas_Humanos <- read_csv("Data_files_csv/Alturas_Humanos.csv")

head(Alturas_Humanos)
GeneroAltura_cmPeso_kg
Hombres17496
Hombres18987
Mujer185110
Mujer195104
Hombres14961
Hombres189104

Ejercicio #2

Datos sobre la orquidea Dipodium

Los datos se encuentra en el paquete ggversa debajo el nombre de dipodium.

  • Evalua todos los pasos para ver si la altura de inflorescencia es un predictor de la cantidad de flores que tiene una planta.