Apuntes y ejercitation
CLASE 1
El curso de estadistica modelo lineal. Porque tiene análisis de la varianza, regresión lineal y regresión logística. Todos estos tienen el mismo principio matemático, por ende entran todos dentro del modelo lineal.
El método de análisis de la varianza:
Cuando uno tiene un conjunto de datos y ver la información que hay allí uno debe aplicar estadística. Estos métodos estadísticos se distinguen entre los descriptivos y lo inferenciales.
Los métodos inferenciales busca extrapolar los resultados obtenidos sobre una muestra hacia la población. Un ejemplo de esto sería el estudio sobre medicamentos, en dónde tomando una muestra se busca poder establecer su eficacia en la población. Al no poder analizar una población completa debemos restringirnos a una muestra.
El método básico de la inferencia estadística es la prueba de hipótesis. Ésta tiene la particularidad de evaluar la veracidad de una hipótesis planteada sobre una población de referencia, dicha información se obtiene a partir de nuestra muestra. El instrumento utilizado para ello es la prueba de hipótesis.
Los test más comúnes de prueba de hipótesis son aquellos que toman como referencia la curva normal comparando dos medias. (Buscar explicación online de test para comparación de dos medias normales)
En los casos en que la distribución de la variable observada no se distribuya normalmente tendrémos que realizar purebas de hipótesis no paramétricas(Test de Matt Wittney).
La razón por la cual uno debe realizar estos supuestos para hacer inferencia estadística es porque necesitamos expandir las conclusiónes de las muestras a toda la población. Siempre que uno trabaje con una prueba de hipótesis se va a tener un estadístico de prueba.
El test de t se utiliza para la comparación de dos medias normales. Sin embargo, cuando uno tiene más de poblaciones para comparar podemos utilizar el analisis de la varianza (ANOVA)
ANOVA
El ANOVA no es una prueba de hipótesis, sino que un conjunto o una metodología a partir de la cual se pueden derivar distintas pruebas de hipótesis.Las pruebas de hipótesis son
Analsis de la varianza de un factor
Análisis de varianza de varios factores
Diseños factoriales
Análisis de la varianza de medida repetidas
La razon por la que usamos ANOVA y no T ya que si bien podríamos comparar varios grupos usando T de a pares, pero se aumentaria el error por cada grupo sometido a observación.
Para evitar eso se debe hacer primero un test global. Ver si hay diferencias entre los grupos. Luego determinan entre que pares de grupos se da.
ANOVA es un test global. Nos muestra si hay o no diferencias entre los grupos sin decirnos entre que grupos esas diferencias se manifiestan. Compara la varianza que hay entre todas las unidades dentro cada grupo con la que hay en los promedios de los grupos.
Los supuestos del ANOVA son:
Que tenga distribución normal.
Que las poblaciones que proveen las muestras tienen varianzas iguales.
Que las muestras sean independientes.
Con ANOVA uno busca tener k poblaciones y k muestras. Uno supone que la muestra extraida proviene de una variable con distribucion normal. Las hipotesis que queremos contrastar son que H0 “las medias son todas iguales” y H1 “Al menos dos de las medias no son iguales”.
Si no hubiera diferencia entre los grupos, las medias deberian ser todas iguales entre sí e iguales a la media total (media de las medias de cada grupo).
SST=SSentre + SSintra
Las hipotesis se plantean en terminos de u(mu, la media de la poblacion) pero el analisis de hace en base a las observaciones de la muestra.
La ANOVA da cuenta de la variabilidad existente entre los grupos y al interior de los grupos.
PRACTICA
EJERCICIO 1
El gabinete psicológico desea determinar si existen diferencias en la cantidad media de horas dormidas por los pacientes con insomnio durante la primera noche en los tres grupos de edad a) Realice un gráfico boxplot que permita comparar las horas dormidas durante la primera noche para los 3 grupos etarios considerados. ¿Qué puede concluir? b) Realice la verificación de los supuestos para realizar un ANOVA.( Normalidad y Homocedasticidad de varianzas) c) ¿Puede afirmar que existen diferencias? Utilizar = 0,05 d) En caso de existir diferencias verificar para que grupos etarios las medias son distintas de los otros grupos (Tukey y comparar con el BoxPlot) e) Repetir los ítems anteriores pero considerando la tercera noche
La primera parte del ejercicio nos pide hacer un gráfico de boxplot para comparar los tres grupos de edad solicitados. Eso nos va a permitir tener una primera aproximación al fenómeno para entender la dispersión en torno a la media de cada uno de los grupos.
Lo primero que realizamos es importar la base “insomnio” para trabajar con los datos que nos pide el ejercicio.
#importación de la base
Insomnio_2 <- read_sav("Insomnio_2.sav")
View(Insomnio_2)Siempre se debe mirar los datos primero desde un ángulo descriptivo. Para ello agrupamos los segmentos de edad y le pedimos que nos resuma la información principal proveniente de los datos. Para ello le pediremos media y desvío usando tidyverse.
#descriptivo por edad
Insomnio_2 %>%
group_by(Edad) %>%
get_summary_stats(PrimerNoche, type = "mean_sd")## # A tibble: 3 x 5
## Edad variable n mean sd
## <dbl+lbl> <chr> <dbl> <dbl> <dbl>
## 1 1 [Menor a 20 años] PrimerNoche 33 5.70 1.69
## 2 2 [Entre 20 y 25 años] PrimerNoche 33 7.74 2.37
## 3 3 [Mas de 25 años] PrimerNoche 33 6.88 1.91Luego pasamos a visualizar mediante el uso de boxplot. Eso nos permite observar la distribución de las observaciones por cada uno de los grupos.
Debajo se puede observar visualmente que no hay observaciones atipicas y que tampoco hay mucha simetria. Las distribuciones no son simetricas y por lo tanto se entiende que no estamos hablando de una distribución normal. Es importante ver si no hay observaciones atipicas ya que las mismas afectan gravemente al modelo.
En caso de tener observaciones atipicas es mejor sacarlas o arreglarlas. En caso de tener grandes muestras deberíamos sacar los casos atípicos los cuales afectan al analisis.
#grafico de boxplot por grupo de edad
ggboxplot(Insomnio_2, x = "Edad", y = "PrimerNoche")
Para realizar el analisis de anova usamos la funcion aov(analysis of variance). El primer argumento que necesitamos es poner donde estan las variables que vamos a usar (data insomnio) y luego determinar cuales son las variables que seran incluidas en el analisis.
El rol que indica la funcion que las variables van a complir en el analisis es ~ (el rulo de muinoz). A la izquierda del rulo se encuentra la variable dependiente, la cual nosotros queremos analizar. A la derecha ponemos la o las variables predictoras, independientes o de corte. En nuestro caso ponemos la variable edad.
En nuestro ejercicio usamos la variable edad factorizada (factor) ya que es el tipo de variables que necesitamos para correr el modelo. Esto se puede solucionar ya sea con la funcion factor o creando una nueva variable que sea factor.
# Anova de un factor
model <- aov( PrimerNoche ~ factor(Edad), data=Insomnio_2)Con la funcion summmary podemos ver el contenido el modelo. La tabla anova muestra el estadistico F y el P value. La tabla solamente no entrega una probabilidad (P value).
Entonces lo que debemos hacer es comparar el P value con el alfa con que nosotros queremos trabajar. Por ejemplo, si queremos trabajar con un alfa de 0,05 y nuestra probabilidad es menor que ella, en este caso rechazamos la H0 porque el P value es 0.000371 ***. Por ende, vamos a rechazar la hipotesis nula da cuenta de que hay diferencia entre los grupos de edad.
summary(model)## Df Sum Sq Mean Sq F value Pr(>F)
## factor(Edad) 2 69.3 34.66 8.588 0.000371 ***
## Residuals 96 387.5 4.04
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Todo esto se basa en la normalidad de la variancia. La prueba de hipotesis que se usan para calcular la normalidad son la de Shapiro y la otra es la de Kolmodoro - Smirnoff. Se parte de que la hipotesis nula de esos test es que la muestra proviene de una muestra normal.
Debajo se ejecuta el test de Shapiro para cada uno de los grupos. Si asumimos mediante la hipotesis nula que las muestras provienen de una distribucion normal, mirando el valor P entonces aceptamos la hipotesis nula. Esto nos dice que las muestras tienen una distribucion normal ya que P Value > 0,05.
# testeo de normalidad de las variables
Insomnio_2 %>%
group_by(Edad) %>%
shapiro_test(PrimerNoche)## # A tibble: 3 x 4
## Edad variable statistic p
## <dbl+lbl> <chr> <dbl> <dbl>
## 1 1 [Menor a 20 años] PrimerNoche 0.973 0.563
## 2 2 [Entre 20 y 25 años] PrimerNoche 0.969 0.458
## 3 3 [Mas de 25 años] PrimerNoche 0.973 0.574Podemos visualizarlo con la siguiente funcion qqplot. Lo que se observa en esta visualizacion es la comparacion de los cuantiles teoricos de una distribucion normal con los cuantiles de una muestra. Si la muestra corresponde a una distribucion normal, entonces los cuantiles teoricos y muestral deberian coincidir.
ggqqplot(Insomnio_2, "PrimerNoche", facet.by = "Edad")
Una vez comparada la normalidad debemos comparar las varianzas. Si las varianzas coinciden entonces decimos que hay homocedasticidad. En caso contrario, decimos que hay homocedasaticidad. Para ello usamos el test de Levene. H0 nos dice que las varianzas son iguales, H1 las varianzas no son iguales. En el ejemplo el P value hace que la hipotesis nula no se rechace. Podemos suponer que las varianzas son iguales.
# testeo de homocedasticidad de varianzas
testeo_homocedasticidad<- Insomnio_2 %>% levene_test(PrimerNoche ~ as.factor(Edad))
testeo_homocedasticidad## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 2 96 1.18 0.313EJERCICIO 2
Un psicólogo clínico está interesado en evaluar la efectividad de las siguientes tres técnicas para tratar la depresión leve: reestructuración cognitiva, capacitación asertiva y un programa de ejercicio y nutrición. Cuarenta estudiantes que sufren de depresión leve son seleccionados al azar de la lista de espera del centro de asesoría de una Universidad y se asignaron aleatoriamente a cada una de las tres técnicas mencionadas anteriormente. Los diez restantes a un grupo control placebo. Se lleva a cabo el tratamiento durante 10 semanas, tras lo cual la depresión se mide mediante el inventario de depresión de Beck (que es un cuestionario auto administrado que consta de 21 preguntas de respuesta múltiple con la finalidad de medir la severidad de una depresión). A continuación se muestran las puntuaciones de depresión obtenidas después del tratamiento. Las puntuaciones más altas indican una mayor depresión.
Realice un gráfico boxplot que permita comparar las puntuaciones medias para los tratamientos considerados ¿Qué puede concluir?
Realice la verificación de los supuestos para realizar un ANOVA. (Normalidad y Homogeneidad de varianzas)
¿Puede afirmar que existen diferencias? Utilizar = 0,05
De existir diferencias realizar los test post hoc y extraer conclusiones.
# Creamos el dataset.
placebo <- c(27,16,18,26,18,28,25,20,24,26)
reestruc_cogni <- c(10,8,14,16,18,8,12,14,9,7)
capac_asertiva <- c(16,18,12,15,9,13,17,20,21,19)
ejercicio_nutric <- c(26,24,17,23,25,22,16,15,18,23)- Realice un gráfico boxplot que permita comparar las puntuaciones medias para los tratamientos considerados ¿Qué puede concluir?
Para poder correr el boxplot es necesario poder tener dos ejes. Por ende, lo que debemos hacer es no solamente tener un vector con los valores, sino que tambien otro vector al que le asignemos la categoria que representa cada uno de los valores. Para asignar las categorias, usaremos la funcion rep y asignaremos el nombre al vector.
#creamos un vector numerico con todos los valores de los tratamientos.
depresion<- c(placebo,reestruc_cogni,capac_asertiva,ejercicio_nutric)
#creamos un vector asignando el nombre a los valores, para ello usamos la funcion rep
tratamiento <- c(rep("placebo",10), rep("reestructuracion",10), rep("capacidad",10), rep("ejercicio",10))
#Ahora si podemos crear el data frame para luego hacer el boxplot. Asi tendriamos dos columnas/vectores/ejes para cruzar entre si en el boxplot
matriz <- data.frame(tratamiento,depresion)
View(matriz)#Boxplot
ggboxplot(matriz, "tratamiento", "depresion")
A partir del grafico de boxplot se puede observar que el tratamiento de reestructuracion arroja mejores resultados para combatir la depresion. Por otro lado, el tratamiento de ejercicios y placebos no responden con la misma eficacia y a su vez presentan una mayor dispersion.
- Realice la verificación de los supuestos para realizar un ANOVA. (Normalidad y Homogeneidad de varianzas)
#Testeo de normalidad
matriz %>% group_by(tratamiento) %>% shapiro_test(depresion)## # A tibble: 4 x 4
## tratamiento variable statistic p
## <chr> <chr> <dbl> <dbl>
## 1 capacidad depresion 0.967 0.860
## 2 ejercicio depresion 0.901 0.225
## 3 placebo depresion 0.887 0.157
## 4 reestructuracion depresion 0.929 0.440#Testeo de homocedasticidad
matriz %>% levene_test(depresion ~ as.factor(tratamiento))## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 3 36 0.110 0.954- ¿Puede afirmar que existen diferencias? Utilizar = 0,05
Para determinar si existen diferencias entre los grupos debemos proceder a realizar ANOVA.
anova_trata <- aov (depresion ~ as.factor(tratamiento), data = matriz )
summary(anova_trata)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(tratamiento) 3 762.9 254.29 15.92 9.4e-07 ***
## Residuals 36 574.9 15.97
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Se puede afirmar que existe diferencias significativas entre los grupos.
- De existir diferencias realizar los test post hoc y extraer conclusiones.
Como se cumplen los supuestos de normalidad y homocedasticidad usamos Tukeys o REWQ.
anova_trata <- aov (depresion ~ as.factor(tratamiento), data = matriz ) %>% tukey_hsd()
anova_trata## # A tibble: 6 x 9
## term group1 group2 null.value estimate conf.low conf.high p.adj
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 as.factor(t~ capaci~ ejercicio 0 4.9 0.0868 9.71 4.47e-2
## 2 as.factor(t~ capaci~ placebo 0 6.8 1.99 11.6 2.86e-3
## 3 as.factor(t~ capaci~ reestruc~ 0 -4.4 -9.21 0.413 8.35e-2
## 4 as.factor(t~ ejerci~ placebo 0 1.90 -2.91 6.71 7.14e-1
## 5 as.factor(t~ ejerci~ reestruc~ 0 -9.3 -14.1 -4.49 4.62e-5
## 6 as.factor(t~ placebo reestruc~ 0 -11.2 -16.0 -6.39 1.79e-6
## # ... with 1 more variable: p.adj.signif <chr>Viendo los resultados del test de Turkey podemos decir que hay diferencia significativas entre capacidad contra ejercicio y placebo, asi como el tratamiento de reestructuracion es significativamente diferente contra el ejercicio y el placebo.
EJERCICIO 3
Un profesor quiere determinar la mejor forma de presentar un importante tema al grupo que tiene a su cargo. Este profesor puede escoger entre una de las tres opciones siguientes: 1) dar clase, 2) dar clase y asignar una lectura complementaria o 3) mostrar una película y asignar una lectura complementaria. Para ello, decide hacer un experimento para evaluar las tres opciones. Solicita a 27 voluntarios del grupo al que imparte clases y los divide de manera aleatoria en 3 grupos de 9 individuos cada uno, asignando a cada grupo una de las tres condiciones. Bajo la condición 1, él da la clase. Bajo la condición 2, él da la clase y asigna una lectura complementaria. Bajo la condición 3, los alumnos presencian una película acerca del tema y reciben la misma lectura complementaria de los estudiantes que están bajo la condición 2. Después se aplica a todos los estudiantes un examen sobre el material. Se obtuvieron los siguientes resultados (porcentaje de respuestas correctas)
- Cuál es la hipótesis nula global?
- Utilice un ANOVA para determinar la conclusión del experimento
- Cuales son los supuestos que debe considerar para resolver estadísticamente el problema?
- Utilice Test post hoc para determinar los grupos con diferencias
La hipotesis nula global dice que no hay diferencias entre las distintas modalidad de impartir un tema al grupo.
b)Utilice un ANOVA para determinar la conclusión del experimento
- Cuales son los supuestos que debe considerar para resolver estadísticamente el problema?
Los supuestos que se deben considerar son los de normalidad y homocedasticidad.
- Utilice Test post hoc para determinar los grupos con diferencias
#creamos el dataset
condicion1 <- c(92,86,87,76,80,87,92,83,84)
condicion2 <- c(86,93,97,81,94,89,98,90,91)
condicion3 <- c(81,80,72,82,83,89,76,88,83)
tablacondi<- data.frame(condicion1,condicion2,condicion3)
View(tablacondi)- Utilice un ANOVA para determinar la conclusión del experimento
#Primero vemos los estadisticos descripttivos
get_summary_stats(tablacondi, type="mean_sd")## # A tibble: 3 x 4
## variable n mean sd
## <chr> <dbl> <dbl> <dbl>
## 1 condicion1 9 85.2 5.22
## 2 condicion2 9 91 5.34
## 3 condicion3 9 81.6 5.32#lo visualizamos. Para ello debemos crear dos vectores
Eje1b <- c(condicion1, condicion2, condicion3)
Eje2b <- c(rep("condicion1",9),rep("condicion2",9),rep("condicion3",9))
Ejes<- data.frame(Eje1b,Eje2b)
View(Ejes)
ggboxplot(Ejes, "Eje2b", "Eje1b")
#Testeo de normalidad * shapiro
normalidadcondiciones<- Ejes %>% group_by(Eje2b) %>% shapiro_test(Eje1b)
View(normalidadcondiciones)
#Visualizamos la normalidad
ggqqplot(data = Ejes, "Eje1b", facet.by = as.factor("Eje2b") )
#Testeo de varianza levene
varianzacondiciones<- Ejes %>% levene_test(Eje1b~ as.factor(Eje2b))
View(varianzacondiciones)
#Grafico de medias por grupo.
#Testeo de anova
anovadecondic<- aov( Eje1b ~ as.factor(Eje2b), Ejes)
summary(anovadecondic)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Eje2b) 2 408.1 204.04 7.289 0.00336 **
## Residuals 24 671.8 27.99
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Cuales son los supuestos que debe considerar para resolver estadísticamente el problema?
Los supuestos que se deben considerar son los de normalidad y homoecedasticidad los cuales se cumplen y por ende, podemos proceder a ejecutar el test de Tukey para ver la diferencia entre los grupos.
- Utilice Test post hoc para determinar los grupos con diferencias
condiposthoc<- aov( Eje1b ~ as.factor(Eje2b), Ejes) %>% tukey_hsd()
condiposthoc## # A tibble: 3 x 9
## term group1 group2 null.value estimate conf.low conf.high p.adj p.adj.signif
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 as.f~ condi~ condi~ 0 5.78 -0.451 12.0 0.0726 ns
## 2 as.f~ condi~ condi~ 0 -3.67 -9.89 2.56 0.323 ns
## 3 as.f~ condi~ condi~ 0 -9.44 -15.7 -3.22 0.0025 **La diferencia entre grupos es significativa solamente entre la condicion 2 y la condicion 3. Por ende, el metodo mas efectivo para el aprendizaje de los alumnos es presentar el tema y dar una clase, en comparacion con la eficacia que provee ver una pelicula.
EJERCICIO 4
Para verificar si la memoria cambia con la edad un investigador realiza un experimento en el cual hay cuatro grupos de seis sujetos cada uno. Los grupos difieren en cuanto a la edad de los sujetos siendo ésta 30, 40, 50 y 60 años respectivamente. Suponga que todos los sujetos gozan de buena salud y que además concuerdan en otras variables importantes como la escolaridad, IQ, género, motivación, etc. Se muestra a cada sujeto una serie de sílabas sin sentido (una combinación sin significado de tres letras como DAF o FUM) a razón de una sílaba cada 4 segundos. La serie se muestra dos veces, después de los cual se pide a los sujetos que escriban el mayor número posible de sílabas que puedan recordar. La cantidad de sílabas recordadas por cada sujeto es la siguiente:
memoria <- data.frame(grupo1,grupo2, grupo3,grupo4)
knitr::kable(memoria,caption = "Cantidad de silabas recordada por grupo")| grupo1 | grupo2 | grupo3 | grupo4 |
|---|---|---|---|
| 14 | 12 | 17 | 13 |
| 13 | 15 | 14 | 10 |
| 15 | 16 | 14 | 7 |
| 17 | 11 | 9 | 8 |
| 12 | 12 | 13 | 6 |
| 20 | 18 | 15 | 9 |
- Realizar un gráfico de medias de los grupos
- Calcular media y desvío de cada grupo
- Utilice el ANOVA para determinar si la edad influye sobre la memoria
- Utilizar test post hoc para determinar cuáles grupos difieren entre sí, en caso de haber diferencias globales
ejememo1<- c(grupo1,grupo2,grupo3, grupo4)
ejememo2 <- c(rep("grupo1",6), rep("grupo2",6),rep("grupo3",6),rep("grupo4",6))
graficmedia<- data.frame(ejememo1,ejememo2)grafdata<- ddply(graficmedia, c(as.factor("ejememo2")), summarise,
N = length(ejememo1),
media = mean(ejememo1),
sd = sd(ejememo1),
error = sd/sqrt(N))| ejememo2 | N | media | sd | error |
|---|---|---|---|---|
| grupo1 | 6 | 15.166667 | 2.926887 | 1.194897 |
| grupo2 | 6 | 14.000000 | 2.756810 | 1.125463 |
| grupo3 | 6 | 13.666667 | 2.658320 | 1.085255 |
| grupo4 | 6 | 8.833333 | 2.483277 | 1.013794 |
#Grafico de medias
ggplot(grafdata,aes(x=ejememo2, y=media)) +
geom_errorbar(aes(ymin=media-error, ymax=media+error),width=.1) +
geom_line()+
geom_point()## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
- Calcular media y desvío de cada grupo
- Utilice el ANOVA para determinar si la edad influye sobre la memoria
anovamemoria<- aov( ejememo1 ~ as.factor(ejememo2), memoria )
summary(anovamemoria)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(ejememo2) 3 140.8 46.94 6.387 0.00327 **
## Residuals 20 147.0 7.35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Efectivamente hay una influencia de la edad sobre la memoria.
- Utilizar test post hoc para determinar cuáles grupos difieren entre sí, en caso de haber diferencias globales
postanovamemoria<- aov( ejememo1 ~ as.factor(ejememo2), memoria ) %>% tukey_hsd()
postanovamemoria## # A tibble: 6 x 9
## term group1 group2 null.value estimate conf.low conf.high p.adj
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 as.factor(ejememo2) grupo1 grupo2 0 -1.17 -5.55 3.21 0.878
## 2 as.factor(ejememo2) grupo1 grupo3 0 -1.50 -5.88 2.88 0.774
## 3 as.factor(ejememo2) grupo1 grupo4 0 -6.33 -10.7 -1.95 0.00325
## 4 as.factor(ejememo2) grupo2 grupo3 0 -0.333 -4.71 4.05 0.996
## 5 as.factor(ejememo2) grupo2 grupo4 0 -5.17 -9.55 -0.786 0.0173
## 6 as.factor(ejememo2) grupo3 grupo4 0 -4.83 -9.21 -0.452 0.0273
## # ... with 1 more variable: p.adj.signif <chr>Soluciones al ejercicio 2 propuesta por el docente.
placebo <- c(27,16,18,26,18,28,25,20,24,26)
rc <- c(10,8,14,16,18,8,12,14,9,7)
ca <- c(16,18,12,15,9,13,17,20,21,19)
eyn <- c(26,24,17,23,25,22,16,15,18,23)depresion <- c(placebo,rc,ca,eyn)
trat <- c(rep("placebo",10),rep("Reestruccog",10),rep("Capacasert",10),rep("Ejercynut",10))
ftrat <- factor(trat, labels=c("Placebo","Reestruccog","Capacasert","Ejercynut"))
eje2 <- data.frame(trat, depresion,ftrat)
eje2## trat depresion ftrat
## 1 placebo 27 Capacasert
## 2 placebo 16 Capacasert
## 3 placebo 18 Capacasert
## 4 placebo 26 Capacasert
## 5 placebo 18 Capacasert
## 6 placebo 28 Capacasert
## 7 placebo 25 Capacasert
## 8 placebo 20 Capacasert
## 9 placebo 24 Capacasert
## 10 placebo 26 Capacasert
## 11 Reestruccog 10 Ejercynut
## 12 Reestruccog 8 Ejercynut
## 13 Reestruccog 14 Ejercynut
## 14 Reestruccog 16 Ejercynut
## 15 Reestruccog 18 Ejercynut
## 16 Reestruccog 8 Ejercynut
## 17 Reestruccog 12 Ejercynut
## 18 Reestruccog 14 Ejercynut
## 19 Reestruccog 9 Ejercynut
## 20 Reestruccog 7 Ejercynut
## 21 Capacasert 16 Placebo
## 22 Capacasert 18 Placebo
## 23 Capacasert 12 Placebo
## 24 Capacasert 15 Placebo
## 25 Capacasert 9 Placebo
## 26 Capacasert 13 Placebo
## 27 Capacasert 17 Placebo
## 28 Capacasert 20 Placebo
## 29 Capacasert 21 Placebo
## 30 Capacasert 19 Placebo
## 31 Ejercynut 26 Reestruccog
## 32 Ejercynut 24 Reestruccog
## 33 Ejercynut 17 Reestruccog
## 34 Ejercynut 23 Reestruccog
## 35 Ejercynut 25 Reestruccog
## 36 Ejercynut 22 Reestruccog
## 37 Ejercynut 16 Reestruccog
## 38 Ejercynut 15 Reestruccog
## 39 Ejercynut 18 Reestruccog
## 40 Ejercynut 23 ReestruccogView(eje2)one.way <- aov(depresion ~ ftrat, data = eje2)summary(one.way)## Df Sum Sq Mean Sq F value Pr(>F)
## ftrat 3 762.9 254.29 15.92 9.4e-07 ***
## Residuals 36 574.9 15.97
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1######################################################
#descriptivo por edad
eje2 %>%
group_by(trat) %>%
get_summary_stats(depresion, type = "mean_sd")## # A tibble: 4 x 5
## trat variable n mean sd
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 Capacasert depresion 10 16 3.80
## 2 Ejercynut depresion 10 20.9 4.01
## 3 placebo depresion 10 22.8 4.37
## 4 Reestruccog depresion 10 11.6 3.78#grafico de boxplot por grupo de edad
ggboxplot(eje2, x = "trat", y = "depresion")
# Anova de un factor
model <- lm(depresion ~ trat, data = eje2)
# testeo de normalidad de las variables
eje2 %>%
group_by(trat) %>%
shapiro_test(depresion)## # A tibble: 4 x 4
## trat variable statistic p
## <chr> <chr> <dbl> <dbl>
## 1 Capacasert depresion 0.967 0.860
## 2 Ejercynut depresion 0.901 0.225
## 3 placebo depresion 0.887 0.157
## 4 Reestruccog depresion 0.929 0.440ggqqplot(eje2, "depresion", facet.by = "trat")
# testeo de homocedasticidad de varianzas
eje2 %>% levene_test(depresion ~ as.factor(trat))## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 3 36 0.110 0.954##################################
#graficos de medias POR GRUPO #
##################################
cdata <- ddply(eje2, c("trat"), summarise,
N = length(depresion),
mean = mean(depresion),
sd = sd(depresion),
se = sd / sqrt(N))
cdata## trat N mean sd se
## 1 Capacasert 10 16.0 3.800585 1.201850
## 2 Ejercynut 10 20.9 4.012481 1.268858
## 3 placebo 10 22.8 4.366539 1.380821
## 4 Reestruccog 10 11.6 3.777124 1.194432ggplot(cdata, aes(x=trat, y=mean)) +
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.1) +
geom_line() +
geom_point()## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggplot(cdata, aes(x=trat, y=mean)) +
geom_point() +
geom_line() ## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
Sobre la clase pasada de ANOVA: Hay otra forma de mostrar graficamente el analisis de la varianza.
anova_apa Esta funcion nos presenta un test F de Snedor para el interceptor (variable dependientes) y para la variable independiente. Incluye el eta cuadraro parcial. Este modo es el que uno escribe en los papers, no pone la tabla de ANOVA simplemente presenta la linea que se desprende de este analisis. La primer fila testa si mu es igual a cero o si es distinta a cero (variable dependiente), la segunda fila nos da la informacion del analisis de la varianza. En lugar de darnos una tabla nos da sus caractristicas en cuanto a los grados de libertad, el valor del estadistico, su nivel de significacion y el tamanio de su efecto(leve, moderado o alto).
El eta cuadrado se puede solo obtener con la funcion de anova_test y no con la de aov. (ver parcial eta square)
CLASE 2
Test post hoc
Aplicacion de un test de analisis de la varianza. Una vez que el anova me dio positivo queremos saber entre que grupos hay diferencias para ello utilizamos el test post hoc.
El ANOVA no tiene la posibilidad de darnos toda la informacion que necesitmaos. Solo nos dice si hay o no diferencia pero no entre que grupos esa diferencia se da o se manifiesta. Por ende, es en si mismo de poca utilidad si no seguimos indagando. Por ende lo que queremos hacer es hacer una comparacion entre los grupos entre si, o sea comparaciones multiples.
Problemas de las comparaciones multiples: Es un problema tecnico estadistico. Si uno hace multiples test del mismo tipo incrementa la probabilidad de encontrar diferencias significativas simplemente por azar o por repeticion (error del tipo 1),
La pregunta es como hacer comparaciones multiples resolviendo el problema que surge del incremento de error del tipo 1. La solucion a este problema son los test post hoc.
Posibilidades:
Viieja escuela: Minima diferencia significativa (LSD): Es el clasido uso de t-test. Recomendado para hasta tres grupos
Bonferroni y Sidak test: Trabajan sobre el test de significacion.
Bonferroni utiliza un alpha prima dividido la cantidad comparaciones que uno necesite. Es decir uno fija un alpha general de comparacion y utilizo para comparar entre los pares de medias un alpha mucho mas chico. Esto me asegura que la probabilidad final de error tipo 1 se mantenga en los niveles que yo quiero. El problema de Bonferroni es que se torna muy exigente, ya que si tengo muchaas comparaciones el alpha sera mucho mas pequenio.
Sidak propone algo similar a Bonferroni, solo con la diferencia de que el alpha utilizado en la comparacion no sea tan exigente.
Otra opcion es trabajar sobre la distribucion del estadistico para hacer las comparaciones multiples sin afectar el alpha. Para ello se utiliza el rango estudentizado. Aparece entonces el test de John Tukey.
Rando Studentizado. La logica es que yo tengo los promedios de cada grupo y se que agarrando la diferencia entre el par de grupos tiene una t de Studente. La diferencia extrema se va a dar cuando yo comparo el grupo del promedio mas grando contra el minimo. Lo cual permite usar la distribucion de la diferencia entre el promedio maximo y minimo lo cual introduce el concepto de orden entre los promedios. Entonces, las distribuciones estadisticas cambian radicalmente. Eso se llama distribucion del rango studentizado. Habiando fijado un alpha, un desvio y un tamanio de muestra, lo que definimos es la diferencia honestamente sifnificativa.
Cual usar?
Si se cumplen los supuestos de normalidad y homoeslaticidad usamos Tukeys o REWQ.
Si tenemos N desiguales usamos Gabril o Hochberg
Si las varianzas son desiguales usamos Games-Howell
Si se compara contra un grupo de control usamos Dunnet.
Cuando no se cumplen los supuestos, el test de anova pierde eficacia. Sin embargo, hay una version robusta. Hay un test de anova de Wells que es utilizado cuando no se tiene igualdad de varianzas.
Practica.
como armar una variable factor en tidyverse!!
Insomnio_2 <- Insomnio_2%>%
mutate(fedad=factor(Edad, labels=(c("menor20","entre 20 y 25","mayor25"))))# estadÃsticos descriptivos de la variable dependiente, en cada grupo
Insomnio_2 %>%
group_by(fedad) %>%
get_summary_stats(PrimerNoche, type = "mean_sd")## # A tibble: 3 x 5
## fedad variable n mean sd
## <fct> <chr> <dbl> <dbl> <dbl>
## 1 menor20 PrimerNoche 33 5.70 1.69
## 2 entre 20 y 25 PrimerNoche 33 7.74 2.37
## 3 mayor25 PrimerNoche 33 6.88 1.91gráficos de la dependiente en cada grupo
ggboxplot(Insomnio_2, x = "fedad", y = "PrimerNoche")
test de normalidad en cada grupo
Insomnio_2 %>%
group_by(fedad) %>%
shapiro_test(PrimerNoche)## # A tibble: 3 x 4
## fedad variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 menor20 PrimerNoche 0.973 0.563
## 2 entre 20 y 25 PrimerNoche 0.969 0.458
## 3 mayor25 PrimerNoche 0.973 0.574qqplot en cada grupo
ggqqplot(Insomnio_2, "PrimerNoche", facet.by = "fedad")
# test de homocedasticidad de varianzas
Insomnio_2 %>% levene_test(PrimerNoche ~ fedad)## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 2 96 1.18 0.313test de anova
original
summary(aov(PrimerNoche~fedad, data=Insomnio_2))## Df Sum Sq Mean Sq F value Pr(>F)
## fedad 2 69.3 34.66 8.588 0.000371 ***
## Residuals 96 387.5 4.04
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1otra en tidy
Insomnio_2 %>% anova_test(PrimerNoche~fedad)## Coefficient covariances computed by hccm()## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 fedad 2 96 8.588 0.000371 * 0.152post hoc
Insomnio_2 %>% t_test(PrimerNoche ~ fedad, p.adjust.method = "bonferroni")## # A tibble: 3 x 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 PrimerNoche menor~ entre~ 33 33 -4.03 57.8 1.63e-4 4.89e-4 ***
## 2 PrimerNoche menor~ mayor~ 33 33 -2.66 63.0 1 e-2 3 e-2 *
## 3 PrimerNoche entre~ mayor~ 33 33 1.63 61.3 1.09e-1 3.27e-1 nsInsomnio_2 %>% t_test(PrimerNoche ~ fedad, p.adjust.method = "holm")## # A tibble: 3 x 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 PrimerNoche menor~ entre~ 33 33 -4.03 57.8 1.63e-4 4.89e-4 ***
## 2 PrimerNoche menor~ mayor~ 33 33 -2.66 63.0 1 e-2 2 e-2 *
## 3 PrimerNoche entre~ mayor~ 33 33 1.63 61.3 1.09e-1 1.09e-1 nsaov(PrimerNoche ~ fedad, data = Insomnio_2) %>% tukey_hsd()## # A tibble: 3 x 9
## term group1 group2 null.value estimate conf.low conf.high p.adj
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 fedad menor20 entre 20 y~ 0 2.04 0.864 3.22 2.29e-4
## 2 fedad menor20 mayor25 0 1.18 0.00139 2.36 4.97e-2
## 3 fedad entre 20 y 25 mayor25 0 -0.863 -2.04 0.315 1.94e-1
## # ... with 1 more variable: p.adj.signif <chr>tabla anova con formato apa
ANOVA ROBUSTO
Insomnio_2 %>% welch_anova_test(PrimerNoche ~ fedad)## # A tibble: 1 x 7
## .y. n statistic DFn DFd p method
## * <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 PrimerNoche 99 8.76 2 62.9 0.000439 Welch ANOVAgames_howell_test(Insomnio_2, PrimerNoche ~ fedad, conf.level = 0.95, detailed = F)## # A tibble: 3 x 8
## .y. group1 group2 estimate conf.low conf.high p.adj p.adj.signif
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 PrimerNoche menor20 entre~ 2.04 0.824 3.26 4.72e-4 ***
## 2 PrimerNoche menor20 mayor~ 1.18 0.113 2.24 2.7 e-2 *
## 3 PrimerNoche entre 20 y 25 mayor~ -0.863 -2.14 0.410 2.42e-1 nsgrafico de medias
m <- Insomnio_2%>%
group_by(fedad)%>%
get_summary_stats(PrimerNoche)
ggplot(m, aes(fedad, mean)) +
geom_point()+
geom_line(line=2)## Warning: Ignoring unknown parameters: line## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
CLASE 3
Anova de dos factores.
El anova de dos factores permite resolver si la determinacion es debido a una variable independiente o a otra, o si la interaccion entre estas dos variables independientes tienen un factor potencial o de atenuacion.
La definicion de los grupos a comparar viene dada por la modalidad de las dos variables. En el sentido tal de que el anova de dos factores va a tomar la combinacion de las dos variables independientes al momento de crear los grupos a comparar. Por ejemplo, si tenemos edad y sexo, un grupo seria hombres jovenes, otro seria mujeres jovenes, otro hombres adultos, mujeres adultas, etc. Asi con cada grupo que se forme a partir de las variables independientes.
Hay tres pruebas de hipotesis que se llevan a cabo en el anova de dos factores.
- Las medias son iguales en todas las categorias del factor uno controlando por factor dos
- Las medias son iguales en todas las categorias del factor dos controlado por factor uno.
- No hay interaccion entre el factor uno y dos. La diferencia entre las medias de un factor, es igual en las categorias de otro factor.
Obtendremos 3 valores P para cada prueba de hipotesis.
Los supuestos para realizar anova de dos factores son:
Normalidad Muestreo probabilistico Varianza
Practica:
Anova de dos factores sin interaccion
View(Insomnio_2)
dd0 <- aov(PrimerNoche ~ Edad+Sexo, data=Insomnio_2)
#el signo mas suma otra variable indep. considerando ambas autonomassummary(dd0)## Df Sum Sq Mean Sq F value Pr(>F)
## Edad 1 22.9 22.931 5.092 0.0263 *
## Sexo 1 1.6 1.636 0.363 0.5481
## Residuals 96 432.3 4.503
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova de dos factores con interaccion
dd <- aov(PrimerNoche ~ Edad*Sexo, data=Insomnio_2)summary(dd)## Df Sum Sq Mean Sq F value Pr(>F)
## Edad 1 22.9 22.931 5.052 0.0269 *
## Sexo 1 1.6 1.636 0.360 0.5497
## Edad:Sexo 1 1.1 1.095 0.241 0.6245
## Residuals 95 431.2 4.539
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Cuando uno dice que los factores interactuan se refiere a que uno de los factores se comporta de manera diferente sobre los factores. La manifestacion del efecto esta condicionada por la presencia del otro factor cuando hay interaccion. Ver graficos de dos ejes compartidos en la clase.
GRAFICO DE INTERACCION
# gráfico de interacción
with(Insomnio_2, interaction.plot(Edad, Sexo, PrimerNoche, fun = mean,
main = "Interaction Plot"))
Hacer el grafico de medias ayuda a entender si nuestro anova de dos factores da cuenta de interaccion o no entre grupos.
Tamanio del efecto
Como determinamos si la relacion tiene un efecto significativo. Uno lo que quiere es tener una medida del tamanio del efecto de un factor dentro de un determinado problema. Una cosa es la significacion estadistica (P value) y otra cosa es la signficacion conceptual que se ve a partir del tamanio del efecto.
Es la importancia de un factor sobre una variable dependiente en el analisis de una varianza. El coeficiciente utilizado mas comunmmente es el Eta cuadrado.
El eta cuadrado es la suma de cuadrado del efecto dividido la suma de cuadrado total. Cuanto explica la variabilidad de los grupos sobre la variabilidad total. Se considera debil cuando es menor a 0.04 y es grande cuando es mayor a 0.34.
El eta cuadrado parcial descarta la suma de cuadrados de los efectos que no van a ser considerados. Se descarta la interaccion del factor que no tomemos en cuenta. Ver ppt. Es debil cuando es menor a 0.02 y grande cuando .
Tablas con resumenes estadisticos
cuadrodosfactores<- Insomnio_2 %>%
group_by(Sexo,Edad) %>%
get_summary_stats(PrimerNoche, type = "mean_sd")# testeo de normalidad de las variables
cuadrodosfactores2<- Insomnio_2 %>%
group_by(Sexo,Edad) %>%
shapiro_test(PrimerNoche)ggqqplot(Insomnio_2, "PrimerNoche", facet.by = "Edad*Sexo")
En el grafico se ve que uno de los grupos no permite aceptar el supuesto de normalidad. Hay valores por fuera de la zona de aceptancia de normalidad (grupo de edad 2 para mujeres)
PRACTICA
EJERCICIO 1
Un gabinete de psicología clínica pretende estudiar la eficacia de cuatro terapias (Psicoanalítica, conductista, cognitivista y gestáltica) en el tratamiento de los trastornos del sueño. Para ello asigna aleatoriamente a un grupo de 24 pacientes a cada terapia y mide las horas que duermen transcurrido un mes después de la terapia. Además se tiene como información si los pacientes son o no fumadores. Los resultados obtenidos son los que se muestran en el archivo terapia.sav
Existe diferencia por sexo en las horas dormidas? Justifique e indique el método estadístico utilizado para responder.
Existe diferencia debida al tipo de terapia recibida en las horas dormidas? Justifique e indique el método estadístico utilizado para responder. Si existen diferencias en horas dormidas, indique entre qué grupos se manifiestan dichas diferencias.
Comparar las medias de las horas dormidas utilizando como factores la terapia y el sexo en conjunto. A su criterio, existe algún vínculo entre ambas variables que afecte a las horas de sueño?
- Analizar la significancia de la interacción. ¿Qué significa? b. Analizar la influencia del sexo sobre las horas dormidas.
- Analizar la influencia de la terapia sobre las horas d. Analizar los correspondientes gráficos de medias 5) Comparar las medias de las horas dormidas tomando como factores el hecho de fumar y el sexo. Utilizar = 0,05 a. Analizar la significancia de la interacción. ¿Qué significa esto? b. Analizar la influencia del sexo sobre las horas dormidas. c. Analizar la influencia de ser fumador en las horas dormidas. d. ¿Alcanza con la tabla de ANOVA para responder b y c? Justifique
terapia<- read_sav(file = "Terapia_2020.sav")View(terapia)terapia$Sexo<- factor(terapia$Sexo, levels=c(1,2), labels= c("Hombre", "Mujer"))
is.factor(terapia$Sexo)## [1] TRUEterapia<- as.data.frame(terapia)
terapia %>% group_by(Sexo) %>% summarise(N = length(Horas),
media= mean(Horas),
desvio = sd(Horas),
error = desvio/sqrt(N))## N media desvio error
## 1 24 5.916667 1.863066 0.3802967resumenterapia<- terapia %>% group_by(Sexo) resumenterapia## # A tibble: 24 x 4
## # Groups: Sexo [2]
## Terapia Sexo Fuma Horas
## <dbl+lbl> <fct> <dbl+lbl> <dbl>
## 1 1 [Psicoanalítica] Hombre 0 [No] 4
## 2 1 [Psicoanalítica] Mujer 1 [Si] 3
## 3 1 [Psicoanalítica] Hombre 0 [No] 5
## 4 1 [Psicoanalítica] Mujer 1 [Si] 4
## 5 1 [Psicoanalítica] Mujer 1 [Si] 4
## 6 1 [Psicoanalítica] Hombre 0 [No] 5
## 7 2 [Conductista] Mujer 1 [Si] 5
## 8 2 [Conductista] Mujer 0 [No] 6
## 9 2 [Conductista] Hombre 1 [Si] 7
## 10 2 [Conductista] Hombre 1 [Si] 7
## # ... with 14 more rows- Existe diferencia por sexo en las horas dormidas? Justifique e indique el método estadístico utilizado para responder.
Para poder justificar el metodo a utilizar debemos primero comprobar si se cumplen las condiciones de normalidad y homoecedasticidad de las muestras hombre y mujer. Para ello procedemos a realizar el test de shapiro y de levene respectivamente…
#Test de shapiro para determinar la normalidad
terapia %>% group_by(Sexo) %>% shapiro_test(Horas)## # A tibble: 2 x 4
## Sexo variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Hombre Horas 0.876 0.0776
## 2 Mujer Horas 0.937 0.466#Para determinar la normalidad de manera visual usamos ggqqplot
ggqqplot(terapia, x = "Horas", facet.by = "Sexo")
Podemos verificar entonces que mediante el test de shapiro o la forma grafica, la normalidad se cumple. Entonces, debemos pasar a comprobar si la varianza entre las muestras cumple la condicion de homoecedasticidad o no. Para ello utilizaremos el test de levene y su forma grafica.
#Test de levene para determinar la homoecedasticidad
terapia %>% levene_test( Horas ~ as.factor(Sexo) )## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 1 22 0.186 0.670#Visualmente como
ggboxplot(terapia, "Sexo", "Horas")
A partir del test de levene podemos comprobar que la varianza de las muestras cumplen la condicion de homoecedasticidad.
Al observar que las condiciones de normalidad y de igualdad de varianzas se cumplen, podemos proceder a hacer el analisis de varianza y aplicar el test post hoc de Tukey.
anovaterapiasexo<- aov(Horas ~ as.factor(Sexo), terapia)
summary(anovaterapiasexo)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Sexo) 1 13.50 13.500 4.477 0.0459 *
## Residuals 22 66.33 3.015
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El test de anova nos arroja que la diferencia de las horas dormidas entre hombres y mujeres es significativa estadisticamente. Debido a que nuestro testo solo tiene dos muestras, cuando corremos el test de Tukey obtenemos los mismos resultados. En caso de tener mas de dos muestras nuestro testeo post hoc de Tukey nos diria entre que grupos la diferencia es significativa.
anovaterapiasexotukey<- aov(Horas ~ as.factor(Sexo), terapia) %>% tukey_hsd()summary(anovaterapiasexo)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Sexo) 1 13.50 13.500 4.477 0.0459 *
## Residuals 22 66.33 3.015
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Existe diferencia debida al tipo de terapia recibida en las horas dormidas? Justifique e indique el método estadístico utilizado para responder. Si existen diferencias en horas dormidas, indique entre qué grupos se manifiestan dichas diferencias.
terapiahorasanova<- aov( Horas ~ as.factor(Terapia), terapia)
summary(terapiahorasanova)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Terapia) 3 54.83 18.28 14.62 2.84e-05 ***
## Residuals 20 25.00 1.25
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hay diferencias significativas entre las variables. Para determinar entre que grupos se muestran esas diferencias debemos comprobar la normalidad y la homocedasticidad para determinar el testo post hoc mas adecuado.
#Testeo de normalidad
terapia %>% group_by(as.factor(Terapia)) %>% shapiro_test(Horas)## # A tibble: 4 x 4
## `as.factor(Terapia)` variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 1 Horas 0.866 0.212
## 2 2 Horas 0.915 0.473
## 3 3 Horas 0.831 0.110
## 4 4 Horas 0.921 0.514#Testeo de normalidad de modo visual
ggqqplot(terapia, "Horas" , facet.by = "Terapia")
#Testeo de homoecedasticidad
terapia %>% levene_test(Horas ~ as.factor(Terapia))## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 3 20 0.422 0.739#Test post hoc de Tukey debido a la normalidad y homocedasticidad
aov(Horas ~ as.factor(Terapia), terapia) %>% tukey_hsd()## # A tibble: 6 x 9
## term group1 group2 null.value estimate conf.low conf.high p.adj
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 as.factor(Terapia) 1 2 0 2.50 0.693 4.31 4.82e-3
## 2 as.factor(Terapia) 1 3 0 3.83 2.03 5.64 4.58e-5
## 3 as.factor(Terapia) 1 4 0 0.667 -1.14 2.47 7.33e-1
## 4 as.factor(Terapia) 2 3 0 1.33 -0.473 3.14 1.98e-1
## 5 as.factor(Terapia) 2 4 0 -1.83 -3.64 -0.0266 4.59e-2
## 6 as.factor(Terapia) 3 4 0 -3.17 -4.97 -1.36 4.59e-4
## # ... with 1 more variable: p.adj.signif <chr>ggboxplot(terapia, "Terapia", "Horas")
Las diferencias entre las horas dormidas y el tipo de tratamiento se observan entre la terapia 1 vs 2 y 3, 2 vs 4 y 3 vs 4.
Comparar las medias de las horas dormidas utilizando como factores la terapia y el sexo en conjunto. A su criterio, existe algún vínculo entre ambas variables que afecte a las horas de sueño?
#Primero debemos determinar si la interaccion de los dos factores se da de manera conjunta o no. Para ello aplicamos el metodo grafico
with(terapia, interaction.plot(as.factor(Sexo), as.factor(Terapia), Horas, fun = mean, main= "Interaction Plot" ) )
El grafico nos muestra que no hay interaccion entre los grupos por ende debemos correr las funciones para anova de dos factores sin interaccion.
sexoterapiahoras<- aov( Horas ~ as.factor(Sexo)+as.factor(Terapia), terapia)
summary(sexoterapiahoras)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Sexo) 1 13.50 13.500 22.3 0.000148 ***
## as.factor(Terapia) 3 54.83 18.278 30.2 1.95e-07 ***
## Residuals 19 11.50 0.605
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ggqqplot(terapia, "Horas", facet.by = "Terapia+Sexo")
#Grafico correcto de medias para la interaccion
with(terapia, interaction.plot(as.factor(Terapia), as.factor(Sexo), Horas, fun=mean, main= "Interaction plot"))
EJERCICIO 2
El Departamento de nutrición de la Ciudad desea analizar los efectos del tipo de harina y del porcentaje de endulzamiento sobre los atributos físicos de un determinado tipo de torta, a la venta actualmente. Para ello, se diseñó un experimento utilizando dos tipos de harina (multipropósito y harina para tortas), y distintos porcentajes de endulzamiento, en cuatro niveles diferentes. Los siguientes datos corresponden a la información de la densidad específica de muestras de tortas preparadas con las distintas combinaciones posibles. Se prepararon 3 tortas con cada una de las combinaciones posibles
Los resultados obtenidos son los que se muestran en el archivo terapia.sav 1) Existen diferencias estadísticamente significativas debidas al tipo de harina utilizada?
harinagen <- c(0.90, 0.87, 0.90, 0.86,0.89,0.91,0.93,0.88,0.87,0.79,0.82,0.80)
harinator <- c(0.91,0.90,0.80,0.88,0.82,0.83,0.86,0.85,0.80,0.86,0.85,0.85)
concazu<- c(rep(0,3), rep(50,3), rep(75,3), rep(100,3))
ejerc2 <- data.frame(concazu,harinagen,harinator)
ejerc2## concazu harinagen harinator
## 1 0 0.90 0.91
## 2 0 0.87 0.90
## 3 0 0.90 0.80
## 4 50 0.86 0.88
## 5 50 0.89 0.82
## 6 50 0.91 0.83
## 7 75 0.93 0.86
## 8 75 0.88 0.85
## 9 75 0.87 0.80
## 10 100 0.79 0.86
## 11 100 0.82 0.85
## 12 100 0.80 0.85resp1 <- c(harinagen,harinator)
resp1fac <- c(rep("general",12), rep( "tortas",12))
ejerc1<- data.frame(resp1fac,resp1)
#Cumple la normalidad segun el test de shapiro
ejerc1 %>% group_by(factor(resp1fac)) %>%
shapiro_test(resp1)## # A tibble: 2 x 4
## `factor(resp1fac)` variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 general resp1 0.924 0.319
## 2 tortas resp1 0.950 0.643#Se cumple la condicion de homoecedasticidad
ejerc1 %>% levene_test(resp1 ~ as.factor(resp1fac))## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 1 22 0.528 0.475#Pasamos a hacer el analisis de varianza
anova1<- aov(resp1 ~ as.factor(resp1fac), ejerc1)
summary(anova1)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(resp1fac) 1 0.00184 0.001837 1.16 0.293
## Residuals 22 0.03486 0.001584ggqqplot(ejerc1, "resp1", facet.by="resp1fac")
- Existen diferencias estadísticamente significativas debidas al nivel de endulzamiento utilizado?
#Usamos la funcion melt para cambiar el formato de la tabal y poder analizar los resultados
tablaendu<- melt(ejerc2, id.vars = "concazu" ) %>% as.data.frame()
colnames(tablaendu)[3] <- "densidad"
colnames(tablaendu)[2] <- "tipoharina"
tablaendu## concazu tipoharina densidad
## 1 0 harinagen 0.90
## 2 0 harinagen 0.87
## 3 0 harinagen 0.90
## 4 50 harinagen 0.86
## 5 50 harinagen 0.89
## 6 50 harinagen 0.91
## 7 75 harinagen 0.93
## 8 75 harinagen 0.88
## 9 75 harinagen 0.87
## 10 100 harinagen 0.79
## 11 100 harinagen 0.82
## 12 100 harinagen 0.80
## 13 0 harinator 0.91
## 14 0 harinator 0.90
## 15 0 harinator 0.80
## 16 50 harinator 0.88
## 17 50 harinator 0.82
## 18 50 harinator 0.83
## 19 75 harinator 0.86
## 20 75 harinator 0.85
## 21 75 harinator 0.80
## 22 100 harinator 0.86
## 23 100 harinator 0.85
## 24 100 harinator 0.85#determinamos la normalidad
tablaendu %>% group_by(as.factor(concazu)) %>% shapiro_test(densidad)## # A tibble: 4 x 4
## `as.factor(concazu)` variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 0 densidad 0.728 0.0122
## 2 50 densidad 0.951 0.748
## 3 75 densidad 0.964 0.851
## 4 100 densidad 0.889 0.310#Visualizacion de normalidad
ggqqplot(tablaendu, "densidad", facet.by = "concazu" )
#determinamos la homoecedasticidad
tablaendu %>% levene_test(densidad ~ as.factor(concazu))## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 3 20 0.0530 0.983#realizamos el analisis de varianza
aov(densidad ~ as.factor(concazu), tablaendu)## Call:
## aov(formula = densidad ~ as.factor(concazu), data = tablaendu)
##
## Terms:
## as.factor(concazu) Residuals
## Sum of Squares 0.00871250 0.02798333
## Deg. of Freedom 3 20
##
## Residual standard error: 0.03740544
## Estimated effects may be unbalanced- Existe interacción estadísticamente significativa entre los factores considerados?
View(tablaendu)with(tablaendu, interaction.plot(concazu,tipoharina, densidad, fun = mean , main="Interacion Plot"))
###anova de dos factores sin interaccion
anovados<- aov(densidad ~ as.factor(tipoharina)+as.factor(concazu), tablaendu)
summary(anovados)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(tipoharina) 1 0.001838 0.001837 1.335 0.262
## as.factor(concazu) 3 0.008713 0.002904 2.110 0.133
## Residuals 19 0.026146 0.001376###anova de dos factores con interaccion. Observando el grafico, estos son los factores que debemos mirar
anovados2 <- aov(densidad ~ as.factor(tipoharina)+as.factor(concazu), tablaendu)
summary(anovados2)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(tipoharina) 1 0.001838 0.001837 1.335 0.262
## as.factor(concazu) 3 0.008713 0.002904 2.110 0.133
## Residuals 19 0.026146 0.001376- Realizar los gráficos de medias e interpretar
EJERCICIO 3.
Un investigador lleva a cabo un experimento para determinar si un medicamento hipnótico llamado Drowson, contra el insomnio ayuda a promover el sueño. Además está interesado en determinar si la tolerancia a la droga se desarrolla con el uso crónico. El diseño del experimento se desarrolló de la siguiente manera: Una de las variables es la concentración de Drowson. Hay dos niveles:
• concentración cero (Placebo)
• dosis recomendada por el fabricante
La otra variable se refiere al uso previo de Drowson, también clasificada en dos niveles:
• sujetos sin uso previo
• consumidores crónicos.
Dieciséis individuos con dificultades para conciliar el sueño y que no han tenido ningún uso previo de Drowson son asignados al azar a las dos condiciones de concentración, tales que hay ocho sujetos en cada condición. Dieciséis usuarios crónicos de Drowson también se asignan al azar a las dos condiciones, ocho sujetos por condición. Todos los sujetos toman el medicamente prescripto en 3 noches consecutivas, y se registra el tiempo en minutos para conciliar el sueño. Los resultados que se muestran en la siguiente tabla son los tiempos medios en cuestión de minutos para conciliar el sueño para cada sujeto, en las 3 noches:
¿Qué puede concluir?
CLASE 4
Correlacion y Regresion lineal
El analisis de la asociacion entre variables depende del nivel de medicion de las mismas
El analsis estaditico de las asociaciones puede darse de :
Manera grafica (Graficos de dispersion)
Metodo estadisticos.
La covarianza permite evaluar asociaciones del tipo lineal, puede tomar cualquier valor y depende de las unidades de medidas de las variables. En estadistica uno siempre realiza primero el analisis de linealidad. Lo mas probable es que uno siempre intente resolver el problema ajustando una recta. Si uno tiene que ajustar algo que no es una recta, entonces uno habla de un polinomio.
Cuando queremos analizar si hay asociacion o no entre las variables debemos analizar la covarianza:
Covarianza mayor a cero es una asociacion lineal directa
Si es menor a cero es una lineal inversa
Si es igual a cero no estan asociadas linealmente
El mayor defecto de la covarianza es que puede tomar cualquier valor por lo cual no permite decir si una asociacion es fuerte. A su vez, el valor de la covarianza depende de las unidades de medidas de las variables. Ej: Si para medir peso y altura usas kgs y cms, tenes un tipo de covarianza. Pero si en lugar de usar kgs y cms, usamos libras y pies nos da una asociacion distinta. Por lo tanto, no es un indicador adecuado para medir asociacion. Es entonces que se introduce el R de pearson, lo cual es una covarianza estandarizada. Esta medida es independiente de las unidades de medida de las variables.
Cuando tenemos un R con signo menos hablamos de una relacion lineal inversa. Cuando el signo es mas, estamos hablando de una asociacion positiva.
Los valores de R entre 1 y - 1, nos da cuenta de que tan alejados los puntos estan de la recta de regresion. Un R igual a 1 da cuenta de una asociacion perfecta ya que estan todos los puntos sobre la pendiente de la linea de regresion.
Uno considera que una correlacion es alta dependiendo del la calidad de las variables que uno esta asociando. En ciencias naturales las variables definen mejor la medicion y por su propia construccion dan cuenta del fenomeno de manera mas fehaciente. En cambio, las ciencias sociales utilizan variables que no definen su medicion con la misma precision y por lo tanto, se les exigue menos.
Hay distintos tipos de asociacion, los cuales pueden ser clasificados en linealidad o en relacion curvilinea.
Analisis de regresion lineal.
Si uno quiere ver si dos variables estan asociadas mira el coeficiente de correlacion, pero si uno busca predecir utiliza un modelo de regresion lineal. Uno busca vincular las variables en funcion de poder predecir o explicar una variable dependiente en terminos de una variable predictora o independiente.
Los modelos de regresion lineal simple estan absados en una ecuacion que vincula dos variables x e y. El vinculo que se propone es de asociacion lineal que viene dado por alpha y beta. Alpha es lo que se conoce como la ordenada al origien y beta es el efecto de la variable independiente x sobre la variable dependiente y. Estos parametros son poblacionales y deben ser estimados utilizando los datos disponibles, es por ello que se agrega un termino de error ya que nuestra suposicion se basa en la incertidumbre. El error es el error del modelo, que se supone aleatorio y de media cero para todo i. Los parametros de alpha y beta deben estimarse a partir de los datos.
yi=a+bxi+ei
La diferencia entre y e yi es el residuo y describe el error en el ajuste del modelo de cada punto.
Los estimadores de a y b de los parametros de la regresion alpha y beta, pueden hallarse mediante el metodo de minimos cuadrados. Se trata de encontrar los estimadores que reduzcan los valores de a y b.
Lo primero que uno debe ver es si los estimadores son significativos. Para testear a y b se utiliza el test de prueba T.
Para hacerlo en R utilizamos la funcion lm.
El estimate intercept es el a de la formula y la linea que le sigue es el b. El t test lo miramos en la columna t value y vemos si es significativo a partir de Prob<0,5. Cuando una muestra es chica, aun con un R alto la prueba de hipotesis podria dar no significativo.
El R2 se conoce como coeficiente de determinacion, es el R al cuadrado y su interpretacion es simple. Da cuenta de la proporcion de varianza de la variable dependiente que es explicada por el modelo.
ansiedad <- c(28,41,35,39,31,42,50,46,45,37)
nota <- c(82,58,63,89,92,64,55,70,51,72)
examen <- data.frame(ansiedad,nota)
plot(ansiedad,nota)
correl <- cor(ansiedad,nota)
correl## [1] -0.6907746regre <- lm(nota~ansiedad,data=examen)
summary(regre)##
## Call:
## lm(formula = nota ~ ansiedad, data = examen)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.885 -7.957 -1.457 7.507 18.829
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 125.8830 21.1109 5.963 0.000337 ***
## ansiedad -1.4285 0.5287 -2.702 0.026986 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.87 on 8 degrees of freedom
## Multiple R-squared: 0.4772, Adjusted R-squared: 0.4118
## F-statistic: 7.301 on 1 and 8 DF, p-value: 0.02699abline(lm(nota ~ ansiedad))
coefficients(regre)## (Intercept) ansiedad
## 125.883049 -1.428504anova(regre)## Analysis of Variance Table
##
## Response: nota
## Df Sum Sq Mean Sq F value Pr(>F)
## ansiedad 1 861.96 861.96 7.3013 0.02699 *
## Residuals 8 944.44 118.06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fitted.values(regre)## 1 2 3 4 5 6 7 8
## 85.88494 67.31439 75.88542 70.17140 81.59943 65.88589 54.45786 60.17187
## 9 10
## 61.60038 73.02841residuals(regre)## 1 2 3 4 5 6
## -3.8849432 -9.3143939 -12.8854167 18.8285985 10.4005682 -1.8858902
## 7 8 9 10
## 0.5421402 9.8281250 -10.6003788 -1.0284091CLASE 5
Regresion lineal multiple
En la mayor parte de los problemas de investigacion en que se aplica el analisis de regresion, se requiere mas de una variable como predictor. La relacion se restrigue a la linealidad porque es mas simple de comprender, aunque relaciones lineales pueden transformarse en relaciones curvilineas.
La complejidad de la mayoria de los mecanismos cientificos es tal que se necesita un modelo de regresion lineal multiple.
- Como hacemos para considerar mas de un predictor?
Desde el punto de vista teorico conceptual es bastante simple. Las variables predictoras pueden ser incluidos como terminos sumatorios en la ecuacion de regresion. Dadas x1…xk, k variables independientes, el modelo de regresion lineal simple es un modelo sumatorio:
Utilizamos el modelo de R2 ajustado (coeficiente de determinacion ajustado) para regresiones lineales multiples. Cada vez que uno agrega una variable al modelo el R2 sera mas grande, por el solo hecho de agregar otro elemento. Ese crecimiento puede que no sea significativo, sin embargo el r2 crecera. Para evitar considerar un modelo como bueno por el simple hecho de tener muchas variables, se utiliza el r2. Este coeficiente considera la cantidad de variables que uno agrega. Ese resultado suele ser menor.
Lo primero que uno hace antes de correr una regresion multiple es correlacionarlas o calcular los graficos de dispersion entre las variables. Eso permite entender como se asocian las variables independientes con respecto a la dependiente.
Debajo se muestra una de las salidas que provee R para funcion lineal.
En la primera linea se observa el R o coeficiente de variacion el cual incluye las seis variables independientes del modelo vs la dependiente. Es el coeficiente de variacion multiple.
El R2 es el coeficiente de determinacion del modelo.
El R2 ajustado es el que muestra el ajuste por la cantidad de variables que uno posee.
El test ANOVA para la regresion lineal simple nos dice la suma de cuadrados de la regresion, del residual y el total. Aqui el test ANOVA contrasta dos hipotesis, H0 vs H1 (La regresion es significativa). En la tabla con el grado de significancia de 0.000 decimos que la regresion es significativa.
Los parametros estimados nos da el valor de la estimacion, el T que permite determinar si ese valor es distitno de cero o no y la signigicacion para ese test. Uno lo que tiene que observar aqui es ver si la variables introducidas en el modelos son significativas. En el ejemplo se ve que las variables X2 y X4 no son predictores de la variable dependiente en presencia de todas las otras. Esas son las variables que uno podria quitar del modelo ya que no son significativas.
Como hacemos para elegir el mejor conjunto de predictores para un modelo?
No siempre el mejor modelo es el que tiene mas predictores. El mejor modelo de regresion depende de los objetivos que se tengan, Muchas veces no existe el mejor modelo. Por eso, existen varios criterios de comparacion entre modelos y es tan importante la seleccion de predictores. El modo en que nuestro modelo se desarrollara depende los objetivos.
Es importante recordar que cualquier variable que se ingrese en la regresion lineal incrementara el R2 asi como tambien el error. Por ende, el criterio que uno tiene para agregar los predictores dependera de dos estrategias:
Metodo de seleccion automatica de a pasos (Stepwise methods)
Metodo de ajuste a todos los posibles modelos (all possible subsets)
Para comparar los modelos podemos usar algunos criterios, esto se conoce como metodos secuendiales de seleccion de modelos(paquete osllr).:
Forward: El paso uno calcula todas las posibles regresiones con un unico predictor. Por ejemplo si uno tiene 6 predictores, calcula todas con cada uno. Paso 2: Habiendo elegido la mejor variable, una inserta la que da el mayor incremento de la regresion en presencia de X1( la variable predictora ya ingresada en el modelo con mayor aporte). Se continua este proceso hasta que ninguna de las variables produce un incremento significativo en la regresion.
Backward: Este modelo es el opuesto. Paso 1, se selecciona y se elimina del modelo aquella variable que da el valor mas pequeño del R2.
Stepwise: Es una modeificacion del forward, que incorpora en cada etapa un chequeo de la eficacia de las variables ingresadas en el modelo en las etapas anteriores. Este chequeo puede derivar en la eliminacion de alguna de estas variables. El ingreso de variables predictoras puede variar la importancia de variables que ya estaban, por ello el stepwise lo que hace es chequear que no quede en el modelo variables no utiles cuando uno ingresa una variable mas.
Estos metodos automaticos no son la panacea pero permite ver mejor la estrategia de modelacion sin caer en el subjetivismo absoluto pero sin olvidarnos tampoco de no caer en el objetivismo absoluto.
Que criterios son los mejores para tomar un modelo?
R2 ajustado:
Cp de Mallows: Es un indicador del sesgo que se introduce al predecir la variable dependiente con un modelo mal especificado. Un modelo adecuado tiene un Cp igual o menor a la cantidad de parametros del modelo.
Indicadores de Akaike: Indicadores de teoria de la informacion. Provee una medida de seleccion de modelos. AIC estima la calidad de un modelo de acuerdo a dos criterios: El nivel de precision de sus estimadores y la cantidad de variables que contiene. SBC o BIC es una alternativa al AIC de Akaike, Se diferencia del AIC por su criterio de penalizacion de los modelos.
Multicolinealidad.
Tener mas de una variable predictora trae algun problema al cual debamos prestar atencion?
Practica:
library(olsrr)## Warning: package 'olsrr' was built under R version 4.1.1##
## Attaching package: 'olsrr'## The following object is masked from 'package:datasets':
##
## riverslibrary(dplyr)
library(ggplot2)
library(gridExtra)## Warning: package 'gridExtra' was built under R version 4.1.1##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinelibrary(purrr)
library(tibble)
library(nortest)## Warning: package 'nortest' was built under R version 4.1.1library(goftest)## Warning: package 'goftest' was built under R version 4.1.1##
## Attaching package: 'goftest'## The following objects are masked from 'package:nortest':
##
## ad.test, cvm.testlibrary(haven)
library(plyr)un_ejemplo_de_regresion <- read_sav("un ejemplo de regresion.sav")
View(un_ejemplo_de_regresion)Debajo se muestra la funcion que utilizamos para regresion lineal simple. La unica diferencia es que se le suman todas las variables predictoras.
# REGRESION MULTIPLE CON lm() y con ols_regress()
model <- lm(y ~ x1 + x2 + x3 + x4 + x5 + x6, data = un_ejemplo_de_regresion)Con esta funcion observamos que las variables X1 X3, X5 y X6. Debajo de la tabla con el estimados, el error estandar, el t valor y la significancia, esta funcion nos muestra e; R2 ajustado para el modelo.
summary(model)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6, data = un_ejemplo_de_regresion)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.4699 -0.8499 0.0141 0.9996 5.3806
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.022e+02 1.245e+01 8.210 1.99e-08 ***
## x1 -2.199e-01 9.959e-02 -2.208 0.03703 *
## x2 -7.238e-02 5.467e-02 -1.324 0.19802
## x3 -2.681e+00 3.749e-01 -7.150 2.17e-07 ***
## x4 -8.442e-04 5.863e-02 -0.014 0.98863
## x5 -3.732e-01 1.207e-01 -3.092 0.00498 **
## x6 3.047e-01 1.372e-01 2.221 0.03606 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.322 on 24 degrees of freedom
## Multiple R-squared: 0.848, Adjusted R-squared: 0.81
## F-statistic: 22.32 on 6 and 24 DF, p-value: 1.023e-08Debajo se muestra la misma regresion solo que en lugar de utilizar la funcion lm, se utiliza ols_regress.
model1 <- ols_regress(y ~ x1 + x2 + x3 + x4 + x5 + x6, data = un_ejemplo_de_regresion)En la practica la funcion ols utiliza lm para calcularse, solo que formatea la salida de otras manera.
Por ejemplo, si queremos ver la correlacion del modelo utilizando las funciones del paquete ols, entonces debemos utilizar el modelo creado a partir de la funcion lm. Abajo se observa como dentro de ols correlation usamos el modelo creado en lm.
Zero order es la correlacion de cada variable independiente con la variable dependiente.
# CORRELACION DE CADA PREDICTOR CON LA DEPENDIENTE
ols_correlations(model)## Correlations
## -------------------------------------------
## Variable Zero Order Partial Part
## -------------------------------------------
## x1 -0.305 -0.411 -0.176
## x2 -0.163 -0.261 -0.105
## x3 -0.862 -0.825 -0.569
## x4 -0.346 -0.003 -0.001
## x5 -0.398 -0.534 -0.246
## x6 -0.237 0.413 0.177
## -------------------------------------------# DIAGNOSTICOS DE COLINEALIDAD
ols_vif_tol(model)## Variables Tolerance VIF
## 1 x1 0.6672145 1.498768
## 2 x2 0.8668311 1.153627
## 3 x3 0.6643874 1.505146
## 4 x4 0.7599631 1.315853
## 5 x5 0.1174186 8.516537
## 6 x6 0.1136534 8.798680Debajo se muestra las estrategias para calcular todos los posibles modelos. Se corre este modelo usando lm. La cantidad de posibles modelos va a dependender de la cantidad de variables predictoras.
En la tabla de los predictores debajo se muestra los posibles modelos.
# TODOS LOS MODELOS POSIBLES
k <- ols_step_all_possible(model)
plot(k)## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
k## Index N Predictors R-Square Adj. R-Square Mallow's Cp
## 3 1 1 x3 0.74338010 0.7345311403 13.519765
## 5 2 1 x5 0.15838344 0.1293621761 105.889559
## 4 3 1 x4 0.11999670 0.0896517619 111.950747
## 1 4 1 x1 0.09277653 0.0614929635 116.248757
## 6 5 1 x6 0.05604592 0.0234957795 122.048447
## 2 6 1 x2 0.02648849 -0.0070808735 126.715506
## 8 7 2 x1 x3 0.76424693 0.7474074256 12.224935
## 17 8 2 x3 x5 0.76142381 0.7443826564 12.670699
## 18 9 2 x3 x6 0.74522106 0.7270225683 15.229081
## 12 10 2 x2 x3 0.74493479 0.7267158412 15.274283
## 16 11 2 x3 x4 0.74338145 0.7250515578 15.519551
## 10 12 2 x1 x5 0.37599543 0.3314236752 73.529064
## 21 13 2 x5 x6 0.28941948 0.2386637282 87.199232
## 11 14 2 x1 x6 0.25998174 0.2071232889 91.847392
## 9 15 2 x1 x4 0.24760963 0.1938674585 93.800923
## 19 16 2 x4 x5 0.21215907 0.1558847212 99.398495
## 14 17 2 x2 x5 0.16685536 0.1073450310 106.551859
## 7 18 2 x1 x2 0.15063534 0.0899664306 109.112969
## 13 19 2 x2 x4 0.14399971 0.0828568276 110.160721
## 20 20 2 x4 x6 0.14331054 0.0821184329 110.269540
## 15 21 2 x2 x6 0.06751590 0.0009098958 122.237360
## 27 22 3 x1 x3 x5 0.81109446 0.7901049580 6.827804
## 40 23 3 x3 x5 x6 0.80998844 0.7888760478 7.002442
## 28 24 3 x1 x3 x6 0.78173017 0.7574779672 11.464366
## 22 25 3 x1 x2 x3 0.77083060 0.7453673364 13.185386
## 26 26 3 x1 x3 x4 0.76562507 0.7395834084 14.007329
## 38 27 3 x3 x4 x5 0.76227983 0.7358664753 14.535536
## 33 28 3 x2 x3 x5 0.76182904 0.7353655984 14.606715
## 34 29 3 x2 x3 x6 0.74615485 0.7179498380 17.081637
## 39 30 3 x3 x4 x6 0.74526986 0.7169665163 17.221375
## 32 31 3 x2 x3 x4 0.74494195 0.7166021675 17.273152
## 29 32 3 x1 x4 x5 0.43845650 0.3760627779 65.666587
## 31 33 3 x1 x5 x6 0.42227346 0.3580816224 68.221856
## 24 34 3 x1 x2 x5 0.40912553 0.3434728125 70.297888
## 30 35 3 x1 x4 x6 0.35682183 0.2853575841 78.556538
## 41 36 3 x4 x5 x6 0.32686542 0.2520726885 83.286594
## 37 37 3 x2 x5 x6 0.32077932 0.2453103606 84.247576
## 23 38 3 x1 x2 x4 0.30753079 0.2305897655 86.339493
## 25 39 3 x1 x2 x6 0.29021246 0.2113471762 89.074022
## 35 40 3 x2 x4 x5 0.22232446 0.1359160617 99.793401
## 36 41 3 x2 x4 x6 0.15778804 0.0642089332 109.983571
## 50 42 4 x1 x3 x5 x6 0.83681815 0.8117132468 4.766086
## 43 43 4 x1 x2 x3 x5 0.81649255 0.7882606383 7.975456
## 54 44 4 x2 x3 x5 x6 0.81584902 0.7875181054 8.077068
## 56 45 4 x3 x4 x5 x6 0.81162125 0.7826399053 8.744625
## 48 46 4 x1 x3 x4 x5 0.81115901 0.7821065501 8.817612
## 44 47 4 x1 x2 x3 x6 0.78622430 0.7533357276 12.754753
## 49 48 4 x1 x3 x4 x6 0.78244195 0.7489714787 13.351978
## 42 49 4 x1 x2 x3 x4 0.77297910 0.7380528027 14.846143
## 52 50 4 x2 x3 x4 x5 0.76260533 0.7260830738 16.484140
## 53 51 4 x2 x3 x4 x6 0.74617451 0.7071244322 19.078534
## 45 52 4 x1 x2 x4 x5 0.47593322 0.3953075564 61.749089
## 51 53 4 x1 x4 x5 x6 0.47235647 0.3911805443 62.313850
## 47 54 4 x1 x2 x5 x6 0.47171966 0.3904457653 62.414401
## 46 55 4 x1 x2 x4 x6 0.39278110 0.2993628100 74.878640
## 55 56 4 x2 x4 x5 x6 0.35917422 0.2605856441 80.185098
## 59 57 5 x1 x2 x3 x5 x6 0.84800181 0.8176021762 5.000207
## 61 58 5 x1 x3 x4 x5 x6 0.83690359 0.8042843054 6.752595
## 62 59 5 x2 x3 x4 x5 x6 0.81712216 0.7805465875 9.876043
## 57 60 5 x1 x2 x3 x4 x5 0.81677090 0.7801250814 9.931505
## 58 61 5 x1 x2 x3 x4 x6 0.78744812 0.7449377489 14.561513
## 60 62 5 x1 x2 x4 x5 x6 0.52421007 0.4290520785 56.126272
## 63 63 6 x1 x2 x3 x4 x5 x6 0.84800313 0.8100039082 7.000000Debajo se utiliza el metodo de los mejores modelos. Se observa que solo debemos comparar seis modelos, los cuales se comparan debajo.
# LOS MEJORES MODELOS DE CADA TAMAÃO
ols_step_best_subset(model)## Best Subsets Regression
## --------------------------------
## Model Index Predictors
## --------------------------------
## 1 x3
## 2 x1 x3
## 3 x1 x3 x5
## 4 x1 x3 x5 x6
## 5 x1 x2 x3 x5 x6
## 6 x1 x2 x3 x4 x5 x6
## --------------------------------
##
## Subsets Regression Summary
## ---------------------------------------------------------------------------------------------------------------------------------
## Adj. Pred
## Model R-Square R-Square R-Square C(p) AIC SBIC SBC MSEP FPE HSP APC
## ---------------------------------------------------------------------------------------------------------------------------------
## 1 0.7434 0.7345 0.7053 13.5198 154.5083 65.4840 158.8103 233.5838 8.0199 0.2691 0.2920
## 2 0.7642 0.7474 0.7054 12.2249 153.8792 64.8457 159.6151 222.5380 7.8621 0.2655 0.2863
## 3 0.8111 0.7901 0.7591 6.8278 149.0115 61.3486 156.1814 185.1748 6.7253 0.2291 0.2449
## 4 0.8368 0.8117 0.7785 4.7661 146.4737 60.4483 155.0776 166.3575 6.2053 0.2137 0.2259
## 5 0.8480 0.8176 0.7867 5.0002 146.2728 61.6283 156.3107 161.4127 6.1782 0.2157 0.2250
## 6 0.8480 0.8100 0.7628 7.0000 148.2725 64.2115 159.7444 168.4292 6.6095 0.2344 0.2407
## ---------------------------------------------------------------------------------------------------------------------------------
## AIC: Akaike Information Criteria
## SBIC: Sawa's Bayesian Information Criteria
## SBC: Schwarz Bayesian Criteria
## MSEP: Estimated error of prediction, assuming multivariate normality
## FPE: Final Prediction Error
## HSP: Hocking's Sp
## APC: Amemiya Prediction CriteriaEl metodo de setpwise nos arroja cual es el mejor modelo, Basicamente lo que nos dice es que nos quedemos con el modelo que utiliza solo una variable X3 para explicar la variacion sobre la variable dependiente.
# stepwise regression
model <- lm(y ~ x1 + x2 + x3 + x4 + x5 + x6 , data = un_ejemplo_de_regresion)
kkk <- ols_step_both_p(model)
plot(kkk)## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
kkk##
## Stepwise Selection Summary
## -------------------------------------------------------------------------------------
## Added/ Adj.
## Step Variable Removed R-Square R-Square C(p) AIC RMSE
## -------------------------------------------------------------------------------------
## 1 x3 addition 0.743 0.735 13.5200 154.5083 2.7448
## -------------------------------------------------------------------------------------El modelo nos dice entonces que yo predigo todo lo bien que puedo con las variables que tengo. En nuestro caso del set de variables, solo con X3 predigo todo lo bien que podria.
#ALGUNOS GRAFICOS DE RESIDUOS
ols_plot_resid_stud_fit(model)
ols_plot_resid_qq(model)
#test de homocedasticidad
ols_test_breusch_pagan(model)##
## Breusch Pagan Test for Heteroskedasticity
## -----------------------------------------
## Ho: the variance is constant
## Ha: the variance is not constant
##
## Data
## -----------------------------
## Response : y
## Variables: fitted values of y
##
## Test Summary
## ----------------------------
## DF = 1
## Chi2 = 0.6056411
## Prob > Chi2 = 0.4364337# graficos de la dependiente
ols_plot_response(model)
CLASE 6
Multicolinealidad
Cuando las variables predictoras presentan una alta correlacion entre si ademas de con la variable dependiente, entonces lo que explican es casi lo mismo. Cuando las variables independientes estan altamente correlacionadas se habla del fenomeno de multicolinealidad.
Efecto de la colinealidad:
Limita el tamanio de R2 porque los predictores correlacionados explican la misma variacion Dificulta la determinacion de la importancia de los predictores Incrementa la varianza de los estimadores de los parametros.
Diagnosticos de multicolinealidad.
VIF o Tolerancia son los indicadores que se calculan para cada una de las variables predictoras que uno incluye en el modelo. Se calcula una regresion de esa variable respecto de los otros predictores incluidos en el modelo, eso da un R2 que se pone en el VIF. Si ese valor es mayor de 10 estamos hablando de colinealidad. La tolerancia es la inversa del VIF.
Debajo se observa el testeo de multicolinealidad a partir de los datos de la clase 5. Se observa que si bien X5 y x6 tienen un VIF alto o tolerancia baja, asi y todo se siguen encontrando dentro de los parametros de aceptabilidad aunque con una alta correlacion.
# DIAGNOSTICOS DE COLINEALIDAD
ols_vif_tol(model)## Variables Tolerance VIF
## 1 x1 0.6672145 1.498768
## 2 x2 0.8668311 1.153627
## 3 x3 0.6643874 1.505146
## 4 x4 0.7599631 1.315853
## 5 x5 0.1174186 8.516537
## 6 x6 0.1136534 8.798680Valores atipicos y puntos de influencias
Una vez que uno tiene un modelo que le resulta eficiente, uno debe analizar una serie de aspectos del modelo para determinar si los supuestos basicos para una regresion multiple se cumplen.
Estos modelos tienen una particularidad, son susceptibles a los valores atipicos.
Los outliers o valores atipicos son valores que tienen una probabilidad muy baja de ocurrir. A veces estos valores responden a errores en el volcano del dato o en su recopilaion, aunque en otros casos puede estar dando cuenta de una situacion anomala o poco frecuente que debe ser analizada en forma puntual.
Los valores atipicos van a tener un impacto fuerte sobre sumas y promedios. Las medidas de tendencia central son facilmente desviadas por los valores atipicos.
En las regresiones, los valores atipicos pueden estar en las variables dependientes tanto como en las independientes. En las dependientes uno las evalua mediante los residuos, mientras que las independientes mediante los puntos de influencia.
Hay dos corrientes sobre como proceder con los outliers. Por un lado, podemos detectarlos y eliminarlis del analisis o por otro lado, uno puede hacer lo que se conoce commo regresiones robustas. En la practica, la regresion robusta no es el metodo mas utilizado aunque es el mas recomendado.
Los valores atipicos en las variables independientes llamados puntos de influencia o leverage points son aquellos datos que ejercen una considerable influencia en el ajuste del modelo y en general, se asocian a observaciones atipicas en los predictores. En estos casos, la presencia de estos valores cambie drasticamente la pendiente de la recta de regresion.
Supuesto del modelo de regresion lineal
La relacion es lineal Varianza constante Los errores tienen distribucion normal con media cero y varianza constante Muestras independientes
Estos supuestos pueden evaluarse analizando los residuos en donde tenemos estrategias graficas y matematicas para ello.
En la forma matematica tenemos los residuos, residuos studentizado y residuo R student.
En la forma grafica se pueden observar de la siguiente forma
(chequear regresion multinivel para mi proyecto de indicadores economicos)
Deteccion de valores atipicos.
La distancia de Cook indica la influencia conjunta de un caso, sea outlieer o de los predictores. Un valor de C mayor que 1 es considerado grande o bien 4/n.
Los DFBETA compara el estimador del parametro incluyendo el caso y excluyendo el caso atipico, si un valor es influyente el DFBETA sera grande permitiendo detectar este tipo de problemas. Este calculo lo realiza sobre los outliers. Los DFFITS calcular lo mismo que el DFBETA pero para la variable predictora. Por ende, se observa sobre los puntos de influencia. Estos diagnosticos se realizan sobre cada uno de los puntos de la muesta..
Debajo se muestran tres funciones para graficar los distintos tipos de residuos. La ultima funcion es la mas recomendable.
#gráfico de los residuos
ols_plot_resid_fit(model, print_plot = TRUE)
ols_plot_resid_stand(model, print_plot = TRUE)
ols_plot_resid_stud_fit(model, print_plot = TRUE)
En el grafico de arriba se puede observar el residuo r de student en donde se muestra la franja entre 2 y 2 en donde hay dos valores que podrian ser valores atipicos. Se utiliza 2 y -2 porque entre ellos esta el 95% de la distribucion de la curva normal. Es por ello que algun caso por fuera de 2 y -2 que puedo llegar a tener algun caso, en el ejemplo son el caso 17 y el caso 15.
# gráfico que detecta outliers y puntos influyentes
ols_plot_resid_lev(model, print_plot = TRUE)
Debajo se muestra el leverage de los puntos de influencia que es el H sub i (Hi) de los valores de la diagonal del modelo de regresion..
# punto de influencia
ols_leverage(model)## [1] 0.13665170 0.20359032 0.24257045 0.26178893 0.25534189 0.08329544
## [7] 0.24663828 0.19496137 0.20202175 0.48791299 0.21940923 0.40318418
## [13] 0.19289603 0.14504530 0.12893677 0.25748479 0.07904722 0.06942112
## [19] 0.25609953 0.43434498 0.21070284 0.25729841 0.10889814 0.22649735
## [25] 0.15426106 0.35063786 0.24059674 0.32532958 0.26284389 0.26056844
## [31] 0.10172346Debajo se muestra el grafico de Cook que nos seniala que observaciones se encuentran en los residuos. Se observa el caso 10 como caso extremo.
# gráficos de la d de Cook}
ols_plot_cooksd_bar(model, print_plot = TRUE)
ols_plot_cooksd_chart(model, print_plot = TRUE)
Debajo se muestra el dffits, el cual nos dice que el caso 10 varia mucho la aceptabilidad del modelo.
# dfits
ols_plot_dffits(model, print_plot = TRUE)
Debajo el dfbeta nos dice cuanto varian los valores del parametro incluyendo o excluyendo cada uno de los casos. Se muestras 4 graficos porque tenemos 4 parametros que son la ordenada al origen, el parametro de x3, el parametro de x1 y el parametro de x5. Con respecto al intercept se ve que el caso 10 tiene fuerte influencia, en x3 y x5 el caso 10 se mantiene dentro de los parametros pero no asi con x1. Por ende, la conclusion que obtenemos es que el valor que toma la variable X1 para el caso 10 impacta en la estimacion de nuestro modelo.
#dfbetas
ols_plot_dfbetas(model, print_plot = TRUE)
Entonces, a partir de este ejemplo puntual uno lo que haria seria testear el modelo con y sin el caso 10. De esa manera podemos comprobar si hay grandes cambios provocados gracias al ajuste. De no proceder de este modo, podemos buscar hacer una regresion robusta. Uno podria sacar ese caso y ver como funciona el modelo pero vale la pena tener en cuenta que uno lo que busca es que el modelo se comporte de modo correcto para cualquier muestra. Si bien uno puede sacar de su muestra este caso extremo, la idea es que nuestro modelo se ajuste a la realidad de manera tal que se compruebe cierta regularidad independientemente de las muestras.
CLASE 7
Regresion logistica
Cuando tenemos una variable dependiente que es cualitativa utilizamos regresion logistica. En la regresion logistica lo que uno busca ver es la presencia o ausencia de la variable dependiente. Si uno hiciera un grafico de dispersion nos seria imposible hacer una linea de regresion, los grupos nos quedarian agrupados sobre dos lineas del eje.
Por lo tanto, lo que podemos hacer es agrupar los datos por intervalos de clase y obtener la media de la proporcion de aparicion del evento como se ve debajo y a partir de ello graficar.
Algunas conclusiones que se derivan de esto son que:
La relacion debe buscarse entre la probabilidad de la variale dependiente y la de los predictores que uno quiera considerar.
La relacion no puede ser lineal, sino mas bien curvilinea. La S alargada que se ve en el grafico es un comportamiento tipico de las funciones de distribucion de probabilidad. Lo que se debe es plantear una relacion entre la probabilidad del evento y los predictores a traves de una funcion que modele adecuadamente este comportamiento.
Por ende debemos recurrir a dos posibles soluciones.:
Modelos logit: Utilizan la funcion logista.
Modelos probit: Utilizan la inversa de la distribucion normal.
En la practica no hay diferencias entre usar un modelo o el otro. Nosotros abordaremos el modelo de regresion losgistica que permitiria linealizar la funcion como se ve en la ultima linea:
Riesgo relativo y odds ratio:
El riesgo o probabilidad se define como la cantidad de ocurrencias del evento con respecto al total de casos, es una probabilidad simple. El riesgo relativo es el cosciente de la probabilidad de x entre el grupo a y b. El riesgo relativo de un suceso se define como la razón entre la probabilidad de que dicho suceso ocurra (p) y la probabilidad de que no ocurra (1-p).
Si queremos comparar la probabilidad del efecto de la variable independiente sobre la dependiente usamos odds ratio. En el numerado de la ecuacion se ubica la probabilidad de la presencia de ocurrencia del factor que uno desea observar y el denominador la ausencia del mismo. Un odds ratio igual uno significa que la presencia de ese evento no cambia segun el factor, en nuestro ejemplo seria que en presencia o ausencia de obesidad no cambia el riesgo de tener hipertension. La lectura del cuadro de abajo seria que ser obeso aumenta la probabilidad de tener hipertension. La probabilidad de tener hipertension siendo obeso se cuadruplica.
Para sacar conclusiones a partir de una prueba de hipotesis a partir de odds ratio diriamos que:
H0 es igual a 1.
H1 es distinto a 1.
Para rechazar H0 utilizando un intervalo de confianza el valor no deberia estar incluido dentro del mismo.
En el ejemplo debajo calculamos la regresion logistica con un predictor dicotomico utilizando la funcion glm.
Cuando vemos la ultima fila del output, observamos que ci low esta por fuera de 1, lo cual nos dice que rechazamos la hipotesis nula y aceptamos la hipotesis alternativa. En este caso diriamos que fumar (smoke) tiene un efecto de riesgo de un doble de aquellos que no fuman.
Ejemplo:
# EJEMPLO USADO EN LA EXPLICACION #
# lectura de la base
ejemplo1hosmer <- read_sav("ejemplo1hosmer.sav")
View(ejemplo1hosmer)# grafico de puntos y recta de regresión lineal
z <- lm(CHD~AGE,data=ejemplo1hosmer)
plot(ejemplo1hosmer$AGE,ejemplo1hosmer$CHD)
abline(z)
En la funcion debajo, le damos la indicacion “binomial” en family para ajustar el modelo a una regresion logistica.
# ajuste del modelo de regresión logistica
eh <- glm(CHD~AGE,data=ejemplo1hosmer,family="binomial")
summary(eh)##
## Call:
## glm(formula = CHD ~ AGE, family = "binomial", data = ejemplo1hosmer)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9718 -0.8456 -0.4576 0.8253 2.2859
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.30945 1.13365 -4.683 2.82e-06 ***
## AGE 0.11092 0.02406 4.610 4.02e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 136.66 on 99 degrees of freedom
## Residual deviance: 107.35 on 98 degrees of freedom
## AIC: 111.35
##
## Number of Fisher Scoring iterations: 4CLASE 8
Mucha de la terminologia de la regresion logistica se relaciona con las ciencias medicas debido a su alto uso en esta disciplina. La regresion logistica puede hacerse con un predictor cualitativo con varios niveles y no necesariamente solo con variables dicotomicas. Lo que se hace es reemplazar la variable cualitativa por una cantidad de variables dicotomicas, las cuales en conjunto puedan simular el comportamiento de las variables cualitativas. Esto es lo que se conoce como variables dummy.
Para hacer esto en R debemos tomar la variable cualitativa y colocarla como factor, asi como de ve debajo con RACE. De esa manera, la regresion logistica se calcula creando las variables dummy para la variable cualitativa original. En el siguiente ejemplo se observa en el output como la variable race que tiene tres niveles “blanco”, “negro” y “otro”, se compara a modo de factor. Debajo del estimador del intercept se encuentran los estimadores correspondientes a las dos variables race. La primera race2 compara a las madres de raza negra vs las blancas y factor race 3 compara las madres de otras raza contra las de raza blanca.
A su vez, se observa el oddsratio en el output. En esta caso nos diria que las madres de raza negra tendrian tres veces mas probabilidades de ocurrencia que las de raza blanca. El p valor nos indica si la variable es predictora o no, en caso de no ser significativa no podriamos tomarla.
En el caso de una regresion logistica con un predictor cuantitativo, podemos proceder de la misma manera simplemente incluyendo la variable cuantitativa . Para ello:
Predictor cuantitativo
Deben tener relación lineal con la distribución logística
Se recomienda centrar o estandarizar los predictores cuantitativos
En el ejemplo debajo la variable cuantitativa “LWT” actua como factor “protector” por su signo negativo.
En la regresion logistica uno no puede utilizar el metodo de los minimos cuadrados para poder realizar estimaciones. Por lo tanto, el metodo que utilizamos es el metodo de maxima verosimilitud, el cual maxima los parametros dandole cierta similitud al metodo de los minimos cuadrados. El criterio es maximizar la probabilidad conjunta de la muestra obtenida. Si uno obtuvo ciertos valores de una variable mediante una muestra, se asume que la muestra obtenida es la que tenia la mayor probabilidad de ser seleccionada u obtenida. La muestra que se obtiene es la mas probable de ser obtenida, entonces los valores obtenidos deben ser los mas probables. Este razonamiento lleva a la funcion de verosimilitud que funciona para varios modelos estadisticos.
Tomando el logaritmo de la funcion de verosimilitud multiplicado por dos obtenemos la deviance. Es un valor tipico de cada modelo que permite la comparacion de deviancias entre modelos a partir del test de coeciente de maxima verosimilitud.
A partir de la deviancia, podemos tomar decisiones similares a las que obtendriamos con los metodos de seleccion de modelos de la regresion lineal multiple. La deviancia nos puede decir que modelos son los mas convenientes para explicar nuestra variable dependiente. Esto permite clasificar a los individuos en funcion de los estimadores del modelo. Es posible calcular la probabilidad de cada sujeto, se puede predecir un valor de Y en terminos de esa probabilidad, y calcular los aciertos. Para esto es necesario decidir a partir de que probabilidad se considera al individuo como que tuvo el evento. Esto se obtiene a partir del modelo logit permitiendo predecir cuando la dependiente vale 1 o vale 0. Segun los modelos, la dependiente vale 1 cuando la probabilidad estimada para ese sujeto vale mas de 0,5.
La tabla de confusion calcula la cantidad de veces que el modelo predice exito cuando la variable es exito (vale 1) y cuando predice fracaso (vale 0). Si el modelo de regresion logistica fuera muy bueno en nuestro ajuste, C y B (como se ve debajo) tendrian que ser 0 (fracaso y fracaso). La sensitividad y especificidad de un modelo son dos conceptos importantes y se espera que ambos sean altos mientras el error de clasificacion del modelo sea bajo.
La curva ROC grafica en conjunto la sensitividad y la especifividad. Lo ideal de la curva ROC seria que el area representada sea igual a uno, cuanto mas cerca de la diagonal este peor es el modelo. Por ende se grafica de la siguiente manera y se toman los criterios de calidad de la clasificacion como los que siguen:
Al igual que el modelo de regresion lineal multiple, tenemos metodos de seleccion automatica de predictores. Como se menciono mas arriba estaran basados en la deviance y en el AIC y sus variantes.
En este tipo de modelos cuando uno tiene una variable ordinal queda a criterio del investigador si la considerara como una variable cualitativa o cuantitativa. Si tiene mas de 5 categorias, las ordinales tienen un comportamiento similar a las variables cuantitativas. El problema es que con estos casos es mejor plantear escenarios posibles y tomar decisiones a partir de ello. Este es un inconveniente clasico en los estudios psicologicos en donde hay variables Likert. En este tipo de escalas hay dos estrategias o agruparlas como si fuera cualitativas o trabajarla como cuantitativas. Lo mejor seria correr una distribucion de frecuencia y a partir de alli proceder.
Practica