Econometric2
Le Anh
9/9/2021
- I. LOGIT MODEL
- 1. Regress a logistic regresion th “attorney” as the y, marital status and loss as x.
- 2. Read the result
- 3. Predict the values of P(y =1) for each observation
- 4.Predict the value of y at marital = 1, loss = 8 using R commands
- 5. Retrieve McFadden R square (pseudo_R=1-ln(Lfull)/ln(Lnull))
- 6. Get predicted p for a new dataset:
- 7.Create the Upper and Lower limits for Confidence interval for the p
- 8.Plot for datanew
- 9.Test for Goodness of fit: Log Likelihood ratio Test
- 9. Margin effects
I. LOGIT MODEL
library(readr)
AUTOBIsim <- read_csv("D:/Kinhteluong_actuary/Actuary_Econometric/AUTOBIsim.csv")## Warning: Missing column names filled in: 'X1' [1]##
## -- Column specification --------------------------------------------------------
## cols(
## X1 = col_double(),
## ATTORNEY = col_double(),
## CLMSEX = col_double(),
## MARITAL = col_double(),
## CLMINSUR = col_double(),
## SEATBELT = col_double(),
## CLMAGE = col_double(),
## LOSS = col_double()
## )autobi=AUTOBIsim
names(autobi)## [1] "X1" "ATTORNEY" "CLMSEX" "MARITAL" "CLMINSUR" "SEATBELT" "CLMAGE"
## [8] "LOSS"summary(autobi)## X1 ATTORNEY CLMSEX MARITAL
## Min. : 1.0 Min. :1.00 Min. :1.000 Min. :1.000
## 1st Qu.: 335.8 1st Qu.:1.00 1st Qu.:1.000 1st Qu.:1.000
## Median : 670.5 Median :1.00 Median :2.000 Median :2.000
## Mean : 670.5 Mean :1.49 Mean :1.581 Mean :1.569
## 3rd Qu.:1005.2 3rd Qu.:2.00 3rd Qu.:2.000 3rd Qu.:2.000
## Max. :1340.0 Max. :2.00 Max. :2.000 Max. :4.000
## CLMINSUR SEATBELT CLMAGE LOSS
## Min. :1.000 Min. :1.000 Min. : 0.100 Min. : 0.010
## 1st Qu.:2.000 1st Qu.:1.000 1st Qu.: 1.800 1st Qu.: 0.630
## Median :2.000 Median :1.000 Median : 3.350 Median : 2.000
## Mean :1.916 Mean :1.016 Mean : 5.801 Mean : 5.682
## 3rd Qu.:2.000 3rd Qu.:1.000 3rd Qu.: 7.000 3rd Qu.: 5.438
## Max. :2.000 Max. :2.000 Max. :161.100 Max. :341.720Review Graphical Analysis:
i.For Scatter plot: Visualise the linear relationship between the predictor and response.
ii.For Box plot: To spot any outlier observations in the variable. Having outliers in your predictor can drastically affect the predictions as they can affect the direction/slope of the line of best fit.
3i.For Density plot: To see the distribution of the predictor variable. Ideally, a close to normal distribution (a bell shaped curve), without being skewed to the left or right is preferred.
library(e1071)
par(mfrow=c(1, 2)) # divide graph area in 2 columns
plot(density(autobi$ATTORNEY), main="Density Plot: ATTORNEY", ylab="Frequency", sub=paste("Skewness:", round(e1071::skewness(autobi$ATTORNEY), 2))) # density plot for 'ATTORNEY'
polygon(density(autobi$ATTORNEY), col="red")
plot(density(autobi$MARITAL), main="Density Plot: MARITAL", ylab="Frequency", sub=paste("Skewness:", round(e1071::skewness(autobi$MARITAL), 2))) # density plot for 'MARITAL'
polygon(density(autobi$MARITAL), col="red")
1. Regress a logistic regresion th “attorney” as the y, marital status and loss as x.
#Create new variable
autobi$ATTORNEY2=ifelse(autobi$ATTORNEY==1,0,1)
table(autobi$ATTORNEY2)##
## 0 1
## 683 657model1=glm(ATTORNEY2~MARITAL+LOSS, data = autobi,family = binomial)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(model1)##
## Call:
## glm(formula = ATTORNEY2 ~ MARITAL + LOSS, family = binomial,
## data = autobi)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3484 -1.1757 -0.1461 1.0558 4.2923
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.50698 0.15866 3.195 0.0014 **
## MARITAL -0.07809 0.08973 -0.870 0.3841
## LOSS -0.10147 0.01282 -7.914 2.49e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1857.1 on 1339 degrees of freedom
## Residual deviance: 1745.3 on 1337 degrees of freedom
## AIC: 1751.3
##
## Number of Fisher Scoring iterations: 62. Read the result
+Very small P value => the β1 is statistically significant
+Very small P value => loss has an impact on the P(hiring an attorney)
+(-) sign: the impact is negative.
+If loss increases by 1 unit, the log(odd) of hiring an attorney will decreases by 0.1013 Or the odd increases by e^0.1013 times
Residual Deviance = -2 LL(model)
Null deviance = -2 LL(null model). null model: the model without any X
3.Are marital status and loss statistically significant?=>marital status and loss aren’t statistically significant. 4.If the signs are as expected?=>The signs of coeficients are positive as We expected.
3. Predict the values of P(y =1) for each observation
model1_probs = data.frame(probs = predict(model1, type="response"))
head(model1_probs)## # A tibble: 6 x 1
## probs
## <dbl>
## 1 0.426
## 2 0.582
## 3 0.540
## 4 0.603
## 5 0.412
## 6 0.583predictedY = predict(model1, autobi, type="response")
head(predictedY)## 1 2 3 4 5 6
## 0.4263938 0.5823654 0.5399234 0.6029389 0.4123116 0.5831157autobi$phat = predict(model1, data=autobi, type = "response")
autobi[1:5,c("phat", "ATTORNEY","LOSS", "MARITAL")]## # A tibble: 5 x 4
## phat ATTORNEY LOSS MARITAL
## <dbl> <dbl> <dbl> <dbl>
## 1 0.426 1 6.38 2
## 2 0.582 2 0.95 1
## 3 0.540 1 1.88 2
## 4 0.603 2 0.11 1
## 5 0.412 1 6.95 24.Predict the value of y at marital = 1, loss = 8 using R commands
predictedY = predict(model1, newdata= data.frame(MARITAL=1,LOSS=8), type="response")
predictedY## 1
## 0.40543465. Retrieve McFadden R square (pseudo_R=1-ln(Lfull)/ln(Lnull))
Pseudo R-squared Mc Fadden = 1−LL(model)/ LL(null) => 0 < Pseudo R-squared Mc Fadden <, the larger the better, other things being the same
McF = 1 - model1$deviance/model1$null.deviance
McF## [1] 0.06024277Note: the larger the Pseudo R-squared the better in condtion other things being equal
6. Get predicted p for a new dataset:
Create a new dataset, called autobin
Note: The option se.fit = TRUE will create standard deviation for the prediction of p
autobinew= with(autobi, data.frame(LOSS=c(1,2,3, 4,5),MARITAL=c(1,2,2,1,1) ))
p = predict.glm(model1, newdata=autobinew, type="response", se.fit = TRUE)
p## $fit
## 1 2 3 4 5
## 0.5811309 0.5368973 0.5115949 0.5057518 0.4803947
##
## $se.fit
## 1 2 3 4 5
## 0.02013285 0.01763597 0.01721615 0.01912253 0.01962356
##
## $residual.scale
## [1] 17.Create the Upper and Lower limits for Confidence interval for the p
UL = p$fit+1.96*p$se.fit
LL = p$fit-1.96*p$se.fit
mean = p$fit
df=data.frame(t=1:1340, LL, UL, mean)## Warning in data.frame(t = 1:1340, LL, UL, mean): row names were found from a
## short variable and have been discardeddf## # A tibble: 1,340 x 4
## t LL UL mean
## <int> <dbl> <dbl> <dbl>
## 1 1 0.542 0.621 0.581
## 2 2 0.502 0.571 0.537
## 3 3 0.478 0.545 0.512
## 4 4 0.468 0.543 0.506
## 5 5 0.442 0.519 0.480
## 6 6 0.542 0.621 0.581
## 7 7 0.502 0.571 0.537
## 8 8 0.478 0.545 0.512
## 9 9 0.468 0.543 0.506
## 10 10 0.442 0.519 0.480
## # ... with 1,330 more rows8.Plot for datanew
plot(df$t, df$mean, ylim=c(0,0.7), xlim=c(1,100) ,type="l")
lines(df$t, df$mean, lwd = 1, col="green")
lines(df$t, df$LL, lwd = 1, col="red")
lines(df$t, df$UL, lwd = 1, col="blue")
9.Test for Goodness of fit: Log Likelihood ratio Test
LRT = 2(LL(bmodel)-LL(null))
Note: The model is statistically significant as P-value is very small
lrt = -(model1$deviance-model1$null.deviance)
pchisq(lrt, df=3, lower.tail=FALSE)## [1] 4.325182e-249. Margin effects
Q: If X increases by 1 unit, what happens to the p?
library(margins)
margins(model1)## # A tibble: 1,340 x 18
## X1 ATTORNEY CLMSEX MARITAL CLMINSUR SEATBELT CLMAGE LOSS ATTORNEY2 phat
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 2 2 2 1 1.8 6.38 0 0.426
## 2 2 2 2 1 2 1 5.6 0.95 1 0.582
## 3 3 1 1 2 2 1 1.4 1.88 0 0.540
## 4 4 2 1 1 2 1 5.4 0.11 1 0.603
## 5 5 1 2 2 2 1 0.9 6.95 0 0.412
## 6 6 2 2 2 2 1 2.3 0.15 1 0.583
## 7 7 1 2 1 2 1 0.6 6.84 0 0.434
## 8 8 1 2 2 2 2 9.3 15.2 0 0.232
## 9 9 2 1 2 2 1 8.2 0.1 1 0.584
## 10 10 2 2 1 2 1 3 0.94 1 0.583
## # ... with 1,330 more rows, and 8 more variables: fitted <dbl>,
## # se.fitted <dbl>, dydx_MARITAL <mrgnlffc>, dydx_LOSS <mrgnlffc>,
## # Var_dydx_MARITAL <dbl>, Var_dydx_LOSS <dbl>, _weights <dbl>,
## # _at_number <int>