Modelo probit.
Orlando Moscote Flórez
2022-04-07
1. El modelo probit.
El modelo probit trata de estimar la probabilidad de que un elemento con determinadas variables se encuentre en uno de dos grupos.
Suponer que un elemento pertenece a uno u otro grupo y depende de un índice de conveniencia no observable \(I_{i}\) que se determina por una o más variables explicatorias. Este índice se puede expresar como:
\(I_{i}=\beta X\)
Sea \(Y=\) si el elemento pertenece al grupo \(1\) y \(Y=0\) si pertenece al grupo \(2\).
Se supone, además, que existe un valor de umbral del índice o valor crítico \(I_{i}^{*}\) tal que si \(I_{i}\) excede a \(I_{i} ^{*}\), el elemento pertenece al grupo \(1\), de lo contrario, pertenece al grupo \(2\).
\(I_{i}\) e \(I_{i}^{*}\) son no observables pero se supone que están distribuidos normalmente con la misma media y varianza. Según esto, la probabilidad de que \(I_{i}^{*}\) sea menor o igual que \(I_{i}\) se puede obtener a partir de la distribución normal:
\(P_{i}=P(Y_{i}=1|X_{i})=P( I_{i}^{*}\leq I_{i})=F(I_{i})\)
\(P(Y_{i}=1|X_{i})=\dfrac{1}{\sqrt{2\pi}}\displaystyle\int_{-\infty}^{\beta X}e^{-\dfrac{y^{2}}{2}}dy\)
\(P(Y_{i}=1|X_{i})=\Phi(\beta X)\)
donde \(Y\) es una variable normal estándar.
Para obtener información sobre \(I_{i}\), el índice de conveniencia, e igual que para \(\beta\), se toma la inversa de la función de distribución acumulada:
\(I_{i}=\Phi^{-1}(P_{i})\)
\(I_{i}=\beta X\)
Donde \(\Phi^{-1}\) es la inversa de la función de densidad normal acumulada.
2. Estimación.
Los parámetros (\(\beta\)) pueden estimarse a partir del método de máxima verosimilitud, y a partir del modelo estimado pueden obtenerse las probabilidades \(P_{i}\).
curve(pnorm(x,0,1),xlim=c(-4,4),lwd=2)
3. Ejemplo 1.
Ejemplo tomado del libro Introducción al análisis de regresión lineal de Montgomery-Peck-Vining. Página 432.
Datos de 20 familias.
x: ingreso.
y: propiedad de la vivienda (1= sí, 0= no)
x<-c(38000,51200,39600,43400,47700,53000,41500,40800,45400,52400,38700,40100,
49500,38000,42000,54000,51700,39400,40900,52800)
y<-c(0,1,0,1,0,0,1,0,1,1,1,0,1,0,1,1,1,0,0,1)3.1. Gráfico de los datos.
plot(x,y,pch=16)
3.2. Modelo.
modelo1=glm(y~x,family=binomial(link="probit"))
summary(modelo1)##
## Call:
## glm(formula = y ~ x, family = binomial(link = "probit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0155 -0.8765 0.5090 0.8175 1.6038
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.234e+00 2.532e+00 -2.067 0.0388 *
## x 1.199e-04 5.677e-05 2.112 0.0347 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27.526 on 19 degrees of freedom
## Residual deviance: 22.479 on 18 degrees of freedom
## AIC: 26.479
##
## Number of Fisher Scoring iterations: 43.3. Coeficientes del modelo.
coe<-coefficients(modelo1);coe## (Intercept) x
## -5.2336234340 0.0001198948c1<-coe[1]
c2<-coe[2]
exp(coefficients(modelo1))## (Intercept) x
## 0.005334162 1.0001199023.4. Intervalo de confianza de los parámetros.
confint(modelo1,level=0.95)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) -1.047116e+01 -0.4993096581
## x 1.468827e-05 0.00023755943.5. Gráfico del modelo.
curve(pnorm(c1+c2*x),xlim=c(30000,60000),lwd=2)
3.6. Predichos.
pred<-predict(modelo1)
pred## 1 2 3 4 5 6
## -0.67762222 0.90498872 -0.48579059 -0.03019047 0.48535703 1.12079931
## 7 8 9 10 11 12
## -0.25799053 -0.34191687 0.20959907 1.04886245 -0.59369588 -0.42584321
## 13 14 15 16 17 18
## 0.70116762 -0.67762222 -0.19804315 1.24069408 0.96493611 -0.50976955
## 19 20
## -0.32992739 1.09682035exp(pred)## 1 2 3 4 5 6 7 8
## 0.5078230 2.4719041 0.6152106 0.9702607 1.6247550 3.0673049 0.7726025 0.7104073
## 9 10 11 12 13 14 15 16
## 1.2331835 2.8544022 0.5522823 0.6532188 2.0161054 0.5078230 0.8203345 3.4580128
## 17 18 19 20
## 2.6246200 0.6006340 0.7189759 2.9946290pre1<-predict(modelo1,type="response");pre1## 1 2 3 4 5 6 7 8
## 0.2490056 0.8172643 0.3135578 0.4879576 0.6862884 0.8688134 0.3982071 0.3662067
## 9 10 11 12 13 14 15 16
## 0.5830097 0.8528793 0.2763578 0.3351111 0.7584008 0.2490056 0.4215057 0.8926406
## 17 18 19 20
## 0.8327116 0.3051065 0.3707274 0.86364003.7. Predicción para ingreso de 38000.
oda<- predict(modelo1,data.frame(x=38000))
oda## 1
## -0.6776222pnorm(oda)## 1
## 0.24900563.8. Curva ROC.
roc(y,pre1,plot = TRUE, legacy.axes = TRUE,
percent = TRUE, xlab = "% Falsos positivos",
ylab = "% verdaderos postivios", col = "red", lwd = 2,
print.auc = TRUE)## Setting levels: control = 0, case = 1## Setting direction: controls < cases
##
## Call:
## roc.default(response = y, predictor = pre1, percent = TRUE, plot = TRUE, legacy.axes = TRUE, xlab = "% Falsos positivos", ylab = "% verdaderos postivios", col = "red", lwd = 2, print.auc = TRUE)
##
## Data: pre1 in 9 controls (y 0) < 11 cases (y 1).
## Area under the curve: 79.8%3.9. Matriz de confusión.
tabla1<-table(true=y,pred=round(fitted(modelo1)))
sum(diag(tabla1))/sum(tabla1)## [1] 0.7matriz1<-confusionMatrix(tabla1);matriz1## Confusion Matrix and Statistics
##
## pred
## true 0 1
## 0 7 2
## 1 4 7
##
## Accuracy : 0.7
## 95% CI : (0.4572, 0.8811)
## No Information Rate : 0.55
## P-Value [Acc > NIR] : 0.1299
##
## Kappa : 0.4059
##
## Mcnemar's Test P-Value : 0.6831
##
## Sensitivity : 0.6364
## Specificity : 0.7778
## Pos Pred Value : 0.7778
## Neg Pred Value : 0.6364
## Prevalence : 0.5500
## Detection Rate : 0.3500
## Detection Prevalence : 0.4500
## Balanced Accuracy : 0.7071
##
## 'Positive' Class : 0
## 4. Ejemplo 2.
Ejemplo tomado del libro Introducción al análisis de regresión lineal de Montgomery-Peck-Vining. Página 434.
Datos de 20 familias.
x1: ingreso.
x2: edad (años) del coche mas viejo.
y: si seis meses después habian comprado otro coche (1= sí, 0= no)
x1<-c(45000,40000,60000,50000,55000,50000,35000,65000,53000,48000,37000,31000,
40000,75000,43000,49000,37500,71000,34000,27000)
x2<-c(2,4,3,2,2,5,7,2,2,1,5,7,4,2,9,2,4,1,5,6)
y2<-c(0,0,1,1,0,1,1,1,0,0,1,1,1,0,1,0,1,0,0,0)4.1. Modelo.
modelo2=glm(y2~x1+x2,family=binomial(link="probit"))
summary(modelo2)##
## Call:
## glm(formula = y2 ~ x1 + x2, family = binomial(link = "probit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5640 -0.8053 -0.1545 0.9542 1.8009
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.350e+00 2.722e+00 -1.598 0.1101
## x1 4.559e-05 3.787e-05 1.204 0.2286
## x2 6.099e-01 2.997e-01 2.035 0.0418 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27.726 on 19 degrees of freedom
## Residual deviance: 20.954 on 17 degrees of freedom
## AIC: 26.954
##
## Number of Fisher Scoring iterations: 64.2. Coeficientes.
com<-coefficients(modelo2)
d1<-com[1]
d2<-com[2]
d3<-com[3]
exp(coefficients(modelo2))## (Intercept) x1 x2
## 0.01291076 1.00004559 1.840292304.3. Intervalo de confianza de los parámetros.
confint(modelo2,level=0.95)## Waiting for profiling to be done...## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## 2.5 % 97.5 %
## (Intercept) -1.028031e+01 0.5559561291
## x1 -2.547525e-05 0.0001236394
## x2 1.027560e-01 1.31612239384.4. Predichos.
pred3<-predict(modelo2)
pred3## 1 2 3 4 5 6
## -1.07822889 -0.08633748 0.21556779 -0.85027147 -0.62231405 0.97950178
## 7 8 9 10 11 12
## 1.51547835 -0.16639921 -0.71349702 -1.55137886 0.38681249 1.33311241
## 13 14 15 16 17 18
## -0.08633748 0.28951563 3.10005906 -0.89586296 -0.20031619 -0.50277472
## 19 20
## 0.25003803 0.54082206exp(pred3)## 1 2 3 4 5 6 7
## 0.3401975 0.9172846 1.2405661 0.4272989 0.5367010 2.6631291 4.5515979
## 8 9 10 11 12 13 14
## 0.8467082 0.4899279 0.2119555 1.4722804 3.7928299 0.9172846 1.3357803
## 15 16 17 18 19 20
## 22.1992623 0.4082551 0.8184719 0.6048500 1.2840743 1.7174181pre1<-predict(modelo2,type="response");pre1## 1 2 3 4 5 6 7
## 0.14046581 0.46559907 0.58533767 0.19758709 0.26686769 0.83633395 0.93517435
## 8 9 10 11 12 13 14
## 0.43392140 0.23776913 0.06040546 0.65055248 0.90875254 0.46559907 0.61390659
## 15 16 17 18 19 20
## 0.99903259 0.18516298 0.42061665 0.30756133 0.59872103 0.705684884.5. Estimación para x1=45000,x2=2.
oma<- predict(modelo2,data.frame(x1=45000,x2=2))
oma## 1
## -1.078229pnorm(oma)## 1
## 0.14046584.6. Curva ROC.
roc(y2,pred3,plot = TRUE, legacy.axes = TRUE,
percent = TRUE, xlab = "% Falsos positivos",
ylab = "% verdaderos postivios", col = "red", lwd = 2,
print.auc = TRUE)## Setting levels: control = 0, case = 1## Setting direction: controls < cases
##
## Call:
## roc.default(response = y2, predictor = pred3, percent = TRUE, plot = TRUE, legacy.axes = TRUE, xlab = "% Falsos positivos", ylab = "% verdaderos postivios", col = "red", lwd = 2, print.auc = TRUE)
##
## Data: pred3 in 10 controls (y2 0) < 10 cases (y2 1).
## Area under the curve: 77.5%4.7. Matriz de confusión.
tabla2<-table(true=y2,pred=round(fitted(modelo2)))
sum(diag(tabla2))/sum(tabla2)## [1] 0.65matriz2<-confusionMatrix(tabla2);matriz2## Confusion Matrix and Statistics
##
## pred
## true 0 1
## 0 7 3
## 1 4 6
##
## Accuracy : 0.65
## 95% CI : (0.4078, 0.8461)
## No Information Rate : 0.55
## P-Value [Acc > NIR] : 0.252
##
## Kappa : 0.3
##
## Mcnemar's Test P-Value : 1.000
##
## Sensitivity : 0.6364
## Specificity : 0.6667
## Pos Pred Value : 0.7000
## Neg Pred Value : 0.6000
## Prevalence : 0.5500
## Detection Rate : 0.3500
## Detection Prevalence : 0.5000
## Balanced Accuracy : 0.6515
##
## 'Positive' Class : 0
## |—|
O.M.F.