¡Descarga ejemplo regresion logistica y más Guías, Proyectos, Investigaciones en PDF de Econometría solo en Docsity!
Ejemplo simple de regresión logística
CONTENIDO:
Regresión logística con variables cuantitativas y cualitativas
OBJETIVO:
1. Mostrar los modos en que se puede introducir la información para el análisis
2. Familiarizar con los términos y funciones usuales del análisis de regresión logística
3. Interpretar los resultados del análisis
Regresión logística
1 .Datos: Formas en que se pueden presentar
1.1.Datos en modo frecuencias relativas
1.2.Datos en modo matriz de éxitos y fracasos
1.3.Datos en modo vector de éxitos y fracasos
2. Presentación de la ecuación del modelo en términos de logits, de odds y de
probabilidad
3. Interpretación de los coeficientes del modelo
4. Otros resultados del análisis del modelo ajustado sin interacción
5. Comparación de modelos mediante anova
6. Valores predichos del modelo sin interacción (Predicción o valores ajustados)
7. Predicción de valores nuevos con el modelo ajustado
8. Representación gráfica del modelo
9. Diagnóstico de residuos
9.1Gráfico de valores ajustados frente a residuos:
10. Resumen de funciones básica usadas en el análisis de regresión logística
CODIGO
Práctica
Regresión logística
1.Datos: Formas en que se pueden presentar
Los datos pueden darse de varios modos, según se presente la información relativa a los
exitos y fracasos de la variable dependiente.
1) Un vector de valores que representan proporciones de éxitos. Nº de éxitos yi
entre el total ( ni=éxitos + fracasos ). En este caso los totales deben introducirse
como el argumento weights.
2) Un vector de 0’s y 1’s (fracasos y éxitos, respectivamente). En este caso no hay
que especificar el argumento weights.
3) Un vector con valores que representan más de dos niveles. En este caso se trata
como en el caso 2), anterior, asumiendo que el nivel más bajo representa el cero
o fracaso y los otros el 1(éxito).
4) Una matriz formada por dos columnas que representan los éxitos y fracasos. En
este caso se asume que la primera columna contiene los éxitos ( yi ) y la segunda
los fracasos ( ni-yi ). Tampoco es necesario el argumento weights.
Ejemplo:
Datos:
Varios grupos de trabajadores (hombres), han estado expuestos a diversas dosis de un
compuesto químico. Los datos se muestran en la tabla siguiente:
Introduzca estos datos
> datos=edit(data.frame())
> datos
dosis numexpos numintoxi 1 1 20 2 2 2 30 6 3 3 30 9 4 4 25 10 5 5 20 12
d=datos
d$frec=d$numintoxi/d$numexpos #proporción de intoxicados (éxitos)
d$frac= d$numexpos-d$numintoxi #número de no intoxicados (fracasos)
> d
dosis numexpos numintoxi frec frac 1 1 20 2 0.1 18 2 2 30 6 0.2 24 3 3 30 9 0.3 21 4 4 25 10 0.4 15 5 5 20 12 0.6 8
Representación gráfica:
plot(d$dosis, d$frec,type="l",col=2,ylab="proporción de intoxicados",xlab="dosis")
title("Relación entre dosis e intoxicación",col.main="red")
Null Deviance: 14. Residual Deviance: 0.2037 AIC: 20.
1.3.Datos en modo vector de éxitos y fracasos
Preparación de los datos en formato vector de éxitos y fracasos
>y=c(rep(1,sum(d$numintoxi)) , rep(0,sum(d$frac)) )
>dosis2=c(rep(d$dosis,d$numintoxi),rep(d$dosis,d$frac))
y [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 [38] 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [112] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 dosis [1] 1 1 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 [38] 5 5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 [75] 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 [112] 4 4 4 4 4 4 5 5 5 5 5 5 5 5
>glm(y~dosis2,family=binomial)
Call: glm(formula = y ~ dosis2, family = binomial)
Coefficients: (Intercept) dosis -2.6814 0.
Degrees of Freedom: 124 Total (i.e. Null); 123 Residual Null Deviance: 155. Residual Deviance: 140.5 AIC: 144.
Observe que sí han cambiado, no obstante, los grados de libertad, el valor de la deviance
y el valor del estadístico AIC. Antes los datos se estructuraban en una tabla de 5 filas.
Ahora en una de 125.
b) Suponga que a iguales niveles de dosis y para iguales tamaños de expuestos
(numexpos) se han obtenido los siguientes intoxicadas para subgrupos de mujeres.
dosis numexpos numintoxi 1 1 20 4 2 2 30 10 3 3 30 16 4 4 25 20 5 5 20 19
Genere un nuevo archivo que contenga la información completa diferenciando por sexo.
> #Data frame con nuevos datos para mujeres añadidos
> d2=rbind(datos,datos)
> d
dosis numexpos numintoxi 1 1 20 2 2 2 30 6 3 3 30 9 4 4 25 10 5 5 20 12 11 1 20 2 21 2 30 6 31 3 30 9 41 4 25 10 51 5 20 12
> d2$numintoxi=c(datos$numintoxi,4,10,16,20,19) #sustituye los valores para mujeres
> d2$sexo=c(rep('h',5),rep('m',5))
> d
dosis numexpos numintoxi sexo 1 1 20 2 h 2 2 30 6 h 3 3 30 9 h
4 4 25 10 h 5 5 20 12 h 11 1 20 4 m 21 2 30 10 m 31 3 30 16 m 41 4 25 20 m 51 5 20 19 m
>d2$sexo=factor(d2$sexo)
> glm(d2$numintoxi/d2$numexpos~d2$dosis+d2$sexo, weights=d2$numexpos, family=binomial)
>#o la orden equivalente
> glm(numintoxi/ numexpos~ dosis+ sexo, weights= numexpos, family=binomial,data=d2)
Call: glm(formula = numintoxi/numexpos ~ dosis + sexo, family = binomial, data = d2, weights = numexpos)
Coefficients: (Intercept) dosis sexom -3.3385 0.7972 1.
Degrees of Freedom: 9 Total (i.e. Null); 7 Residual Null Deviance: 69. Residual Deviance: 3.803 AIC: 42.
> summary( glm(numintoxi/ numexpos~ dosis+ sexo+dosis*sexo, weights= numexpos, family =binomial,
data=d2))
Call: glm(formula = numintoxi/numexpos ~ dosis + sexo + dosis * sexo, family = binomial, data = d2, weights = numexpos)
Deviance Residuals: Min 1Q Median 3Q Max -0.56947 -0.13528 0.06966 0.19441 0.
Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.68138 0.59066 -4.540 5.63e-06 *** dosis 0.60134 0.16726 3.595 0.000324 *** sexom -0.01132 0.82765 -0.014 0. dosis:sexom 0.41042 0.25600 1.603 0.
Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 69.0618 on 9 degrees of freedom Residual deviance: 1.1771 on 6 degrees of freedom AIC: 41.
Number of Fisher Scoring iterations: 4
Nota: el mismo resultado da si se sustituye la fórmula anterior (más extensa) por la más
simple: numintoxi/numexpos~dosis*sexo, dado que se asume que los efectos de orden
inferior al de interacción se incluyen por defecto (modelos jerárquicos).
En la tabla anterior, eliminamos el término menos significativo (p-valor más grande)
que es dosis:sexom (recuerde que aunque sexom presenta un p-valor mayor, no puede
salir del modelo antes que otro término de orden superior del que forma parte (principio
de jerarquía). La eliminación del término de interacción da lugar a un nuevo contexto,
dentro del cual vemos que la variable sexo recupera importancia.
> summary( glm(numintoxi/ numexpos~ dosis+ sexo, weights= numexpos, family=binomial,data=d2))
Call: glm(formula = numintoxi/numexpos ~ dosis + sexo, family = binomial, data = d2, weights = numexpos)
Deviance Residuals: Min 1Q Median 3Q Max
Degrees of Freedom: 9 Total (i.e. Null); 7 Residual Null Deviance: 69. Residual Deviance: 3.803 AIC: 42.
Coeficiente de dosis:
1. 0,7972 es el cambio esperado en el logit al aumentar una unidad la dosis, supuestas
estables el resto de las variables en el modelo.
2. exp(0,7972) es la razón de Odds al aumentar una unidad la dosis, supuestas estables
el resto de las variables en el modelo.
exp(0.7972)= 2.
Es decir, al aumentar una unidad la dosis, la Odds o ventaja de intoxicación se duplica;
es 2.22 (podríamos decir que el riesgo de intoxicación frente a no intoxicación se
multiplica por 2.22)
Coeficiente de sexo:
1.2467 es el cambio esperado en el logit al pasar de un hombre a una mujer.
La razón de odds que compara mujeres con hombres es igual a exp(1.2467)= 3.478844. Es
decir, es 3,5 veces superior la odds o ventaja de intoxicación en la mujer que en el hombre.
4. Otros resultados del análisis del modelo ajustado sin interacción
> logit1=glm(numintoxi/numexpos~dosis+sexo, weights=numexpos, family=binomial,data=d2)
> residuals(logit1)
> fitted.values(logit1)
> coef(logit1)
(Intercept) d2$dosis d2$sexom -3.3385342 0.7971514 1.
> deviance(logit1)
[1] 3.
> logit1$null.deviance
[1] 69.
> logit1$y
> logit1$xlevels
$d2$sexo [1] "h" "m"
> logit1$residuals
5. Comparación de modelos mediante anova
La función anova permite comparar modelos anidados. Cuando se usa un solo modelo
se determina la significatividad de cada término añadido.
> anova(logit1,test="Chisq")
Analysis of Deviance Table
Model: binomial, link: logit
Response: numintoxi/numexpos
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev P(>|Chi|) NULL 9 69. dosis 1 46.997 8 22.065 7.109e- sexo 1 18.262 7 3.803 1.925e-
Comprobaremos de otro modo que el término interacción no es significativo:
> logit2=update(logit1,.~.+sexo*dosis)
> logit
Call: glm(formula = numintoxi/numexpos ~ dosis + sexo + dosis:sexo, family = binomial, data = d2, weights = numexpos)
Coefficients: (Intercept) dosis sexom dosis:sexom -2.68138 0.60134 -0.01132 0.
Degrees of Freedom: 9 Total (i.e. Null); 6 Residual Null Deviance: 69. Residual Deviance: 1.177 AIC: 41.
> anova(logit2,test="Chisq")
Analysis of Deviance Table
Model: binomial, link: logit
Response: numintoxi/numexpos
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev P(>|Chi|) NULL 9 69. dosis 1 46.997 8 22.065 7.109e- sexo 1 18.262 7 3.803 1.925e- dosis:sexo 1 2.626 6 1.177 0.
La última, fila de la tabla anterior, correspondiente al término de interacción presenta un
p-valor mayor que 0.05. No se puede rechazar la hipótesis de nulidad o no importancia
del término de interacción para explicar la probabilidad de éxito (intoxicación)
La reducción de la deviance en 2.626 con 1 g.l. no es importante (p-valor=0.105 > 0.05)
Otro modo de comparar los modelos es mediante
> anova(logit1,logit2,test = "Chisq")
Analysis of Deviance Table
Model 1: numintoxi/numexpos ~ dosis + sexo Model 2: numintoxi/numexpos ~ dosis + sexo + dosis:sexo Resid. Df Resid. Dev Df Deviance P(>|Chi|)
h m
Probabilidades de intoxicación por sexo
sexo
probabilidad
#predicción de valores nuevos para hombres
ndosis=c(0.5,7)
sexo=c("h","h")
ndatosh=data.frame(ndosis,sexo) #datos con valores en las variables independientes
ndatosh
ph=predict(logit1, newdata = ndatosh,type = "response")
#predicción de valores nuevos para mujeres
ndosis=c(0.5,7)
sexo=c("m","m")
ndatosm=data.frame(ndosis,sexo)
ndatosm
pm=predict(logit1, newdata = ndatosm,type = "response")
Representación gráfica de la curva
>#Valores ajustados
>plot(probabilidad~dosis, data=d2,subset=sexo=="h",type="l",lty=1,col=2,ylim=c(0,1),xlim=c(0,7))
>lines(probabilidad~dosis, data=d2,subset=sexo=="m",type="l",lty=2,col=4,ylim=c(0,1),xlim=c(0,7))
>#Valores nuevos ajustados
>points(ndosis,pm,col=4,pch=1)
>points(ndosis,ph,col=2,pch=2)
title("Probabilidades de intoxicación por sexo",col.main=3)
legend(0,1,legend=c("hombre","mujer"),col=c(2,4),pch=2:1)
0 1 2 3 4 5 6 7
dosis
probabilidad
Probabilidades de intoxicación por sexo
hombre mujer
9. Diagnóstico de residuos
La representación gráfica de los valores ajustados y los residuos tipo “deviance”,
permite ver si existen o no anomalías en el ajuste.
p=predict(logit1)#valores de los logits ajustados
r=residuals(logit1)
plot(p,r,xlab="valores ajustados",ylab="residuos",ylim=c(-1,1)*max(abs(r)),pch=levels(d2$sexo))
abline(0,0,col="red")
title("Valores ajustados frente a residuos",col.main="red")
Valores ajustados frente a residuos
-2 -1 0 1 2
-1.
-0.
valores ajustados
residuos
Los residuos no caen fuera de la banda (-2, 2). Prácticamente están en la banda (-1,1),
lo que es buen indicio de ajuste; sin embargo, parece existir cierta pauta de descenso
logit1=glm(numintoxi/numexpos~dosis+sexo, weights=numexpos, family=binomial,data=d2)
residuals(logit1)
fitted.values(logit1)
coef(logit1)
deviance(logit1)
logit1$null.deviance
logit1$y
logit1$residuals
anova(logit1,test="Chisq")
logit2=update(logit1,.~.+sexo*dosis)
predict(logit1) #predice los valores de los logits (valores de los logits ajustados)
predict(logit1,type="response") #predice las probabilidades de intoxicación
predict(logit1, newdata = ndatos1,type = "response")
plot(probabilidad~dosis, data=d2,subset=sexo=="h",type="l",lty=1,col=2,ylim=c(0,1),xlim=c(0,7))
lines(probabilidad~dosis, data=d2,subset=sexo=="m",type="l",lty=2,col=4,ylim=c(0,1),xlim=c(0,7))
points(ndosis,pm,col=4,pch=1)
Prácticas:
1) Con la opción copiar y pegar seleccione y copie en el editor bloc de notas los
datos del data.frame d en un archivo de datos denominado tóxico. Léalo
desde R.
2) Con la función de R:( write.table) guarde los datos del data data.frame d en
un archivo denominado toxico2. Léalo luego desde R.
3) Efectúe una predicción para una mujer expuesta a una dosis de 9
4) Razón de odds que compara hombres con mujeres
5) Incremento esperado en el logit al aumentar 2 unidades la dosis de
exposición
6) Probabilidad de intoxicación de una mujer expuesta a una dosis de 15
