Modelo probit.

Orlando Moscote Flórez

2022-04-07

 1. El modelo probit.
 2. Estimación.
 3. Ejemplo 1.
o 3.1. Gráfico de los datos.
o 3.2. Modelo.
o 3.3. Coeficientes del modelo.
o 3.4. Intervalo de confianza de los parámetros.
o 3.5. Gráfico del modelo.
o 3.6. Predichos.
o 3.7. Predicción para ingreso de 38000.
o 3.8. Curva ROC.
o 3.9. Matriz de confusión.
 4. Ejemplo 2.
o 4.1. Modelo.
o 4.2. Coeficientes.
o 4.3. Intervalo de confianza de los parámetros.
o 4.4. Predichos.
o 4.5. Estimación para x1=45000,x2=2.
o 4.6. Curva ROC.
o 4.7. Matriz de confusión.

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 IiIi que se determina por una o más variables explicatorias. Este
índice se puede expresar como:
Ii=βXIi=βX
Sea Y=Y= si el elemento pertenece al grupo 11 y Y=0Y=0 si pertenece al grupo 22.
Se supone, además, que existe un valor de umbral del índice o valor crítico I∗iIi∗ tal que
si IiIi excede a I∗iIi∗, el elemento pertenece al grupo 11, de lo contrario, pertenece al
grupo 22.
IiIi e I∗iIi∗ 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∗iIi∗ sea menor o igual
que IiIi se puede obtener a partir de la distribución normal:
Pi=P(Yi=1|Xi)=P(I∗i≤Ii)=F(Ii)Pi=P(Yi=1|Xi)=P(Ii∗≤Ii)=F(Ii)
P(Yi=1|Xi)=12π−−√∫βX−∞e−y22dyP(Yi=1|Xi)=12π∫−∞βXe−y22dy
P(Yi=1|Xi)=Φ(βX)P(Yi=1|Xi)=Φ(βX)
donde YY es una variable normal estándar.
Para obtener información sobre IiIi, el índice de conveniencia, e igual que para ββ, se toma la
inversa de la función de distribución acumulada:
Ii=Φ−1(Pi)Ii=Φ−1(Pi)
Ii=βXIi=βX
Donde Φ−1Φ−1 es la inversa de la función de densidad normal acumulada.

2. Estimación.
Los parámetros (ββ) pueden estimarse a partir del método de máxima verosimilitud, y a partir
del modelo estimado pueden obtenerse las probabilidades PiPi.
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: 4

3.3. Coeficientes del modelo.

coe<-coefficients(modelo1);coe
## (Intercept) x
## -5.2336234340 0.0001198948
c1<-coe[1]
c2<-coe[2]
exp(coefficients(modelo1))
## (Intercept) x
## 0.005334162 1.000119902
3.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.0002375594

3.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.09682035
exp(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.9946290
pre1<-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.8636400

3.7. Predicción para ingreso de 38000.

oda<- predict(modelo1,data.frame(x=38000))
oda
## 1
## -0.6776222
pnorm(oda)
## 1
## 0.2490056

3.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.7
matriz1<-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: 6

4.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.84029230

4.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.3161223938

4.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.54082206
exp(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.7174181
pre1<-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.70568488

4.5. Estimación para x1=45000,x2=2.

oma<- predict(modelo2,data.frame(x1=45000,x2=2))
oma
## 1
## -1.078229
pnorm(oma)
## 1
## 0.1404658

4.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.65
matriz2<-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
##