library(radiant.data)
library(texreg)
data(titanic)
head(titanic)
## pclass survived sex age sibsp parch fare name cabin embarked
## 1 1st Yes female 29.0000 0 0 211.3375 Allen, Miss. Elisabeth Walton B5 Southampton
## 2 1st Yes male 0.9167 1 2 151.5500 Allison, Master. Hudson Trevor C22 C26 Southampton
## 3 1st No female 2.0000 1 2 151.5500 Allison, Miss. Helen Loraine C22 C26 Southampton
## 4 1st No male 30.0000 1 2 151.5500 Allison, Mr. Hudson Joshua Crei C22 C26 Southampton
## 5 1st No female 25.0000 1 2 151.5500 Allison, Mrs. Hudson J C (Bessi C22 C26 Southampton
## 6 1st Yes male 48.0000 0 0 26.5500 Anderson, Mr. Harry E12 Southampton
lm0 <- lm(fare ~ age + sex + pclass, data = titanic)
summary(lm0)
##
## Call:
## lm(formula = fare ~ age + sex + pclass, data = titanic)
##
## Residuals:
## Min 1Q Median 3Q Max
## -87.12 -12.66 -2.52 4.02 424.66
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 109.2997 4.8989 22.311 < 0.0000000000000002 ***
## age -0.2727 0.1043 -2.614 0.00909 **
## sexmale -11.8135 2.8513 -4.143 0.000037 ***
## pclass2nd -72.2467 3.8881 -18.581 < 0.0000000000000002 ***
## pclass3rd -81.4502 3.6302 -22.437 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 43.6 on 1038 degrees of freedom
## Multiple R-squared: 0.3908, Adjusted R-squared: 0.3884
## F-statistic: 166.4 on 4 and 1038 DF, p-value: < 0.00000000000000022
lm1 <- lm(fare ~ age + sex*pclass, data = titanic)
summary(lm1)
##
## Call:
## lm(formula = fare ~ age + sex * pclass, data = titanic)
##
## Residuals:
## Min 1Q Median 3Q Max
## -90.95 -8.94 -3.94 4.05 436.17
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 122.2950 5.3331 22.931 < 0.0000000000000002 ***
## age -0.2663 0.1028 -2.590 0.00974 **
## sexmale -36.5516 5.1482 -7.100 0.00000000000232 ***
## pclass2nd -91.7052 5.7392 -15.979 < 0.0000000000000002 ***
## pclass3rd -101.7318 5.3394 -19.053 < 0.0000000000000002 ***
## sexmale:pclass2nd 35.1016 7.4790 4.693 0.00000304768544 ***
## sexmale:pclass3rd 34.9785 6.6161 5.287 0.00000015168935 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 42.97 on 1036 degrees of freedom
## Multiple R-squared: 0.4095, Adjusted R-squared: 0.4061
## F-statistic: 119.7 on 6 and 1036 DF, p-value: < 0.00000000000000022
htmlreg(list(lm0, lm1))
##
## <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">
## <table cellspacing="0" align="center" style="border: none;">
## <caption align="bottom" style="margin-top:0.3em;">Statistical models</caption>
## <tr>
## <th style="text-align: left; border-top: 2px solid black; border-bottom: 1px solid black; padding-right: 12px;"><b></b></th>
## <th style="text-align: left; border-top: 2px solid black; border-bottom: 1px solid black; padding-right: 12px;"><b>Model 1</b></th>
## <th style="text-align: left; border-top: 2px solid black; border-bottom: 1px solid black; padding-right: 12px;"><b>Model 2</b></th>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">(Intercept)</td>
## <td style="padding-right: 12px; border: none;">109.30<sup style="vertical-align: 0px;">***</sup></td>
## <td style="padding-right: 12px; border: none;">122.29<sup style="vertical-align: 0px;">***</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(4.90)</td>
## <td style="padding-right: 12px; border: none;">(5.33)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">age</td>
## <td style="padding-right: 12px; border: none;">-0.27<sup style="vertical-align: 0px;">**</sup></td>
## <td style="padding-right: 12px; border: none;">-0.27<sup style="vertical-align: 0px;">**</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.10)</td>
## <td style="padding-right: 12px; border: none;">(0.10)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">sexmale</td>
## <td style="padding-right: 12px; border: none;">-11.81<sup style="vertical-align: 0px;">***</sup></td>
## <td style="padding-right: 12px; border: none;">-36.55<sup style="vertical-align: 0px;">***</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(2.85)</td>
## <td style="padding-right: 12px; border: none;">(5.15)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">pclass2nd</td>
## <td style="padding-right: 12px; border: none;">-72.25<sup style="vertical-align: 0px;">***</sup></td>
## <td style="padding-right: 12px; border: none;">-91.71<sup style="vertical-align: 0px;">***</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(3.89)</td>
## <td style="padding-right: 12px; border: none;">(5.74)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">pclass3rd</td>
## <td style="padding-right: 12px; border: none;">-81.45<sup style="vertical-align: 0px;">***</sup></td>
## <td style="padding-right: 12px; border: none;">-101.73<sup style="vertical-align: 0px;">***</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(3.63)</td>
## <td style="padding-right: 12px; border: none;">(5.34)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">sexmale:pclass2nd</td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">35.10<sup style="vertical-align: 0px;">***</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(7.48)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">sexmale:pclass3rd</td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">34.98<sup style="vertical-align: 0px;">***</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(6.62)</td>
## </tr>
## <tr>
## <td style="border-top: 1px solid black;">R<sup style="vertical-align: 0px;">2</sup></td>
## <td style="border-top: 1px solid black;">0.39</td>
## <td style="border-top: 1px solid black;">0.41</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">Adj. R<sup style="vertical-align: 0px;">2</sup></td>
## <td style="padding-right: 12px; border: none;">0.39</td>
## <td style="padding-right: 12px; border: none;">0.41</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">Num. obs.</td>
## <td style="padding-right: 12px; border: none;">1043</td>
## <td style="padding-right: 12px; border: none;">1043</td>
## </tr>
## <tr>
## <td style="border-bottom: 2px solid black;">RMSE</td>
## <td style="border-bottom: 2px solid black;">43.60</td>
## <td style="border-bottom: 2px solid black;">42.97</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;" colspan="4"><span style="font-size:0.8em"><sup style="vertical-align: 0px;">***</sup>p < 0.001, <sup style="vertical-align: 0px;">**</sup>p < 0.01, <sup style="vertical-align: 0px;">*</sup>p < 0.05</span></td>
## </tr>
## </table>
titanic <- titanic %>%
mutate(surv = as.integer(survived))
titanic<- titanic %>%
mutate(survival = sjmisc::rec(surv, rec = "2=0; 1=1"))%>%
select(pclass, survived, survival, everything())
head(titanic)
## pclass survived survival sex age sibsp parch fare name cabin embarked surv
## 1 1st Yes 1 female 29.0000 0 0 211.3375 Allen, Miss. Elisabeth Walton B5 Southampton 1
## 2 1st Yes 1 male 0.9167 1 2 151.5500 Allison, Master. Hudson Trevor C22 C26 Southampton 1
## 3 1st No 0 female 2.0000 1 2 151.5500 Allison, Miss. Helen Loraine C22 C26 Southampton 2
## 4 1st No 0 male 30.0000 1 2 151.5500 Allison, Mr. Hudson Joshua Crei C22 C26 Southampton 2
## 5 1st No 0 female 25.0000 1 2 151.5500 Allison, Mrs. Hudson J C (Bessi C22 C26 Southampton 2
## 6 1st Yes 1 male 48.0000 0 0 26.5500 Anderson, Mr. Harry E12 Southampton 1
nm0 <- glm(survival ~ age + sex + pclass + fare, family = binomial, data = titanic)
summary(nm0)
##
## Call:
## glm(formula = survival ~ age + sex + pclass + fare, family = binomial,
## data = titanic)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6431 -0.6982 -0.4384 0.6706 2.3949
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.4784425 0.3759530 9.252 < 0.0000000000000002 ***
## age -0.0343766 0.0063680 -5.398 0.0000000673 ***
## sexmale -2.4878826 0.1671013 -14.888 < 0.0000000000000002 ***
## pclass2nd -1.2510727 0.2552347 -4.902 0.0000009503 ***
## pclass3rd -2.2555724 0.2655107 -8.495 < 0.0000000000000002 ***
## fare 0.0003464 0.0017796 0.195 0.846
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1409.99 on 1042 degrees of freedom
## Residual deviance: 981.69 on 1037 degrees of freedom
## AIC: 993.69
##
## Number of Fisher Scoring iterations: 4
nm1 <- glm(survival ~ age + sex*pclass + fare, family = binomial, data = titanic)
summary(nm1)
##
## Call:
## glm(formula = survival ~ age + sex * pclass + fare, family = binomial,
## data = titanic)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0634 -0.6641 -0.4943 0.4336 2.4941
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.8978215 0.6131092 7.988 0.00000000000000137 ***
## age -0.0387245 0.0068044 -5.691 0.00000001262563680 ***
## sexmale -3.8996177 0.5015659 -7.775 0.00000000000000755 ***
## pclass2nd -1.5923247 0.6024844 -2.643 0.00822 **
## pclass3rd -4.1382715 0.5601819 -7.387 0.00000000000014976 ***
## fare -0.0009058 0.0020509 -0.442 0.65874
## sexmale:pclass2nd -0.0600076 0.6372949 -0.094 0.92498
## sexmale:pclass3rd 2.5019110 0.5479901 4.566 0.00000498035051247 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1409.99 on 1042 degrees of freedom
## Residual deviance: 931.45 on 1035 degrees of freedom
## AIC: 947.45
##
## Number of Fisher Scoring iterations: 5
nm2 <- glm(survival ~ age + sex*pclass + fare + I(fare^2), family = binomial, data = titanic)
summary(nm2)
##
## Call:
## glm(formula = survival ~ age + sex * pclass + fare + I(fare^2),
## family = binomial, data = titanic)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.9910 -0.6698 -0.4810 0.4311 2.4911
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.82312425 0.72093499 8.077 0.000000000000000663 ***
## age -0.04174592 0.00699881 -5.965 0.000000002450641030 ***
## sexmale -4.15467555 0.52059154 -7.981 0.000000000000001455 ***
## pclass2nd -2.14286950 0.64569176 -3.319 0.000904 ***
## pclass3rd -4.82128298 0.62929970 -7.661 0.000000000000018399 ***
## fare -0.01384608 0.00543044 -2.550 0.010781 *
## I(fare^2) 0.00003571 0.00001611 2.217 0.026615 *
## sexmale:pclass2nd 0.14735788 0.64832374 0.227 0.820198
## sexmale:pclass3rd 2.73660984 0.56404801 4.852 0.000001223878466529 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1409.99 on 1042 degrees of freedom
## Residual deviance: 923.48 on 1034 degrees of freedom
## AIC: 941.48
##
## Number of Fisher Scoring iterations: 6
library(texreg)
library(pander)
htmlreg(list(nm0, nm1, nm2))
##
## <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">
## <table cellspacing="0" align="center" style="border: none;">
## <caption align="bottom" style="margin-top:0.3em;">Statistical models</caption>
## <tr>
## <th style="text-align: left; border-top: 2px solid black; border-bottom: 1px solid black; padding-right: 12px;"><b></b></th>
## <th style="text-align: left; border-top: 2px solid black; border-bottom: 1px solid black; padding-right: 12px;"><b>Model 1</b></th>
## <th style="text-align: left; border-top: 2px solid black; border-bottom: 1px solid black; padding-right: 12px;"><b>Model 2</b></th>
## <th style="text-align: left; border-top: 2px solid black; border-bottom: 1px solid black; padding-right: 12px;"><b>Model 3</b></th>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">(Intercept)</td>
## <td style="padding-right: 12px; border: none;">3.48<sup style="vertical-align: 0px;">***</sup></td>
## <td style="padding-right: 12px; border: none;">4.90<sup style="vertical-align: 0px;">***</sup></td>
## <td style="padding-right: 12px; border: none;">5.82<sup style="vertical-align: 0px;">***</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.38)</td>
## <td style="padding-right: 12px; border: none;">(0.61)</td>
## <td style="padding-right: 12px; border: none;">(0.72)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">age</td>
## <td style="padding-right: 12px; border: none;">-0.03<sup style="vertical-align: 0px;">***</sup></td>
## <td style="padding-right: 12px; border: none;">-0.04<sup style="vertical-align: 0px;">***</sup></td>
## <td style="padding-right: 12px; border: none;">-0.04<sup style="vertical-align: 0px;">***</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.01)</td>
## <td style="padding-right: 12px; border: none;">(0.01)</td>
## <td style="padding-right: 12px; border: none;">(0.01)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">sexmale</td>
## <td style="padding-right: 12px; border: none;">-2.49<sup style="vertical-align: 0px;">***</sup></td>
## <td style="padding-right: 12px; border: none;">-3.90<sup style="vertical-align: 0px;">***</sup></td>
## <td style="padding-right: 12px; border: none;">-4.15<sup style="vertical-align: 0px;">***</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.17)</td>
## <td style="padding-right: 12px; border: none;">(0.50)</td>
## <td style="padding-right: 12px; border: none;">(0.52)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">pclass2nd</td>
## <td style="padding-right: 12px; border: none;">-1.25<sup style="vertical-align: 0px;">***</sup></td>
## <td style="padding-right: 12px; border: none;">-1.59<sup style="vertical-align: 0px;">**</sup></td>
## <td style="padding-right: 12px; border: none;">-2.14<sup style="vertical-align: 0px;">***</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.26)</td>
## <td style="padding-right: 12px; border: none;">(0.60)</td>
## <td style="padding-right: 12px; border: none;">(0.65)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">pclass3rd</td>
## <td style="padding-right: 12px; border: none;">-2.26<sup style="vertical-align: 0px;">***</sup></td>
## <td style="padding-right: 12px; border: none;">-4.14<sup style="vertical-align: 0px;">***</sup></td>
## <td style="padding-right: 12px; border: none;">-4.82<sup style="vertical-align: 0px;">***</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.27)</td>
## <td style="padding-right: 12px; border: none;">(0.56)</td>
## <td style="padding-right: 12px; border: none;">(0.63)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">fare</td>
## <td style="padding-right: 12px; border: none;">0.00</td>
## <td style="padding-right: 12px; border: none;">-0.00</td>
## <td style="padding-right: 12px; border: none;">-0.01<sup style="vertical-align: 0px;">*</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.00)</td>
## <td style="padding-right: 12px; border: none;">(0.00)</td>
## <td style="padding-right: 12px; border: none;">(0.01)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">sexmale:pclass2nd</td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">-0.06</td>
## <td style="padding-right: 12px; border: none;">0.15</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.64)</td>
## <td style="padding-right: 12px; border: none;">(0.65)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">sexmale:pclass3rd</td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">2.50<sup style="vertical-align: 0px;">***</sup></td>
## <td style="padding-right: 12px; border: none;">2.74<sup style="vertical-align: 0px;">***</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.55)</td>
## <td style="padding-right: 12px; border: none;">(0.56)</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">I(fare^2)</td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">0.00<sup style="vertical-align: 0px;">*</sup></td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;"></td>
## <td style="padding-right: 12px; border: none;">(0.00)</td>
## </tr>
## <tr>
## <td style="border-top: 1px solid black;">AIC</td>
## <td style="border-top: 1px solid black;">993.69</td>
## <td style="border-top: 1px solid black;">947.45</td>
## <td style="border-top: 1px solid black;">941.48</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">BIC</td>
## <td style="padding-right: 12px; border: none;">1023.39</td>
## <td style="padding-right: 12px; border: none;">987.05</td>
## <td style="padding-right: 12px; border: none;">986.03</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">Log Likelihood</td>
## <td style="padding-right: 12px; border: none;">-490.84</td>
## <td style="padding-right: 12px; border: none;">-465.73</td>
## <td style="padding-right: 12px; border: none;">-461.74</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;">Deviance</td>
## <td style="padding-right: 12px; border: none;">981.69</td>
## <td style="padding-right: 12px; border: none;">931.45</td>
## <td style="padding-right: 12px; border: none;">923.48</td>
## </tr>
## <tr>
## <td style="border-bottom: 2px solid black;">Num. obs.</td>
## <td style="border-bottom: 2px solid black;">1043</td>
## <td style="border-bottom: 2px solid black;">1043</td>
## <td style="border-bottom: 2px solid black;">1043</td>
## </tr>
## <tr>
## <td style="padding-right: 12px; border: none;" colspan="5"><span style="font-size:0.8em"><sup style="vertical-align: 0px;">***</sup>p < 0.001, <sup style="vertical-align: 0px;">**</sup>p < 0.01, <sup style="vertical-align: 0px;">*</sup>p < 0.05</span></td>
## </tr>
## </table>
library(Zelig)
z.out <- zelig(survival ~ age + sex + pclass + fare, model = "logit", data = titanic, cite = F)
summary(z.out)
## Model:
##
## Call:
## z5$zelig(formula = survival ~ age + sex + pclass + fare, data = titanic)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6431 -0.6982 -0.4384 0.6706 2.3949
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.4784425 0.3759530 9.252 < 0.0000000000000002
## age -0.0343766 0.0063680 -5.398 0.0000000673
## sexmale -2.4878826 0.1671013 -14.888 < 0.0000000000000002
## pclass2nd -1.2510727 0.2552347 -4.902 0.0000009503
## pclass3rd -2.2555724 0.2655107 -8.495 < 0.0000000000000002
## fare 0.0003464 0.0017796 0.195 0.846
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1409.99 on 1042 degrees of freedom
## Residual deviance: 981.69 on 1037 degrees of freedom
## AIC: 993.69
##
## Number of Fisher Scoring iterations: 4
##
## Next step: Use 'setx' method
x.out <- setx(z.out)
s.out <- sim(z.out, x = x.out)
summary(s.out)
##
## sim x :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.09440502 0.01401926 0.09353095 0.06994333 0.1248277
## pv
## 0 1
## [1,] 0.913 0.087
plot(s.out)
z_tit <- zelig(survival ~ age + sex*pclass + fare, model = "logit", data = titanic, cite = F)
summary(z_tit)
## Model:
##
## Call:
## z5$zelig(formula = survival ~ age + sex * pclass + fare, data = titanic)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0634 -0.6641 -0.4943 0.4336 2.4941
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.8978215 0.6131092 7.988 0.00000000000000137
## age -0.0387245 0.0068044 -5.691 0.00000001262563680
## sexmale -3.8996177 0.5015659 -7.775 0.00000000000000755
## pclass2nd -1.5923247 0.6024844 -2.643 0.00822
## pclass3rd -4.1382715 0.5601819 -7.387 0.00000000000014976
## fare -0.0009058 0.0020509 -0.442 0.65874
## sexmale:pclass2nd -0.0600076 0.6372949 -0.094 0.92498
## sexmale:pclass3rd 2.5019110 0.5479901 4.566 0.00000498035051247
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1409.99 on 1042 degrees of freedom
## Residual deviance: 931.45 on 1035 degrees of freedom
## AIC: 947.45
##
## Number of Fisher Scoring iterations: 5
##
## Next step: Use 'setx' method
a.range = min(titanic$age):max(titanic$age)
x <- setx(z_tit, age = a.range)
s <- sim(z_tit, x = x)
ci.plot(s)
f.range = min(titanic$fare):max(titanic$fare)
x <- setx(z_tit, fare = f.range)
s <- sim(z_tit, x = x)
ci.plot(s)
x <- setx(z_tit, sex = "male")
x1 <- setx(z_tit, sex = "female")
s <- sim(z_tit, x = x, x1 = x1)
summary(s)
##
## sim x :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.1412239 0.01963338 0.1403358 0.1061013 0.1834262
## pv
## 0 1
## [1,] 0.852 0.148
##
## sim x1 :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## [1,] 0.3945576 0.04361669 0.3932977 0.30772 0.4831456
## pv
## 0 1
## [1,] 0.616 0.384
## fd
## mean sd 50% 2.5% 97.5%
## [1,] 0.2533337 0.04412305 0.2530665 0.1692193 0.3430113
plot(s)
fd <- s$get_qi(xvalue="x1", qi="fd")
summary(fd)
## V1
## Min. :0.1340
## 1st Qu.:0.2230
## Median :0.2531
## Mean :0.2533
## 3rd Qu.:0.2813
## Max. :0.3849
c1x <- setx(z_tit, sex = "male", pclass = "1st")
c1x1 <- setx(z_tit, sex = "female", pclass = "1st")
c1s <- sim(z_tit, x = c1x, x1 = c1x1)
c2x <- setx(z_tit, sex = "male", pclass = "2nd")
c2x1 <- setx(z_tit, sex = "female", pclass = "2nd")
c2s <- sim(z_tit, x = c2x, x1 = c2x1)
c3x <- setx(z_tit, sex = "male", pclass = "3rd")
c3x1 <- setx(z_tit, sex = "female", pclass = "3rd")
c3s <- sim(z_tit, x = c3x, x1 = c3x1)
plot(c1s)
plot(c2s)
plot(c3s)
d1 <- c1s$get_qi(xvalue="x1", qi="fd")
d2 <- c2s$get_qi(xvalue="x1", qi="fd")
d3 <- c3s$get_qi(xvalue="x1", qi="fd")
dfd <- as.data.frame(cbind(d1, d2, d3))
head(dfd)
## V1 V2 V3
## 1 0.5482568 0.7741060 0.2844886
## 2 0.4950737 0.7410563 0.2421762
## 3 0.5956448 0.7444520 0.1960172
## 4 0.4132343 0.6944552 0.2294445
## 5 0.4715408 0.7369695 0.2721561
## 6 0.5087496 0.6995612 0.2604529
library(tidyr)
tidd <- dfd %>%
gather(class, simv, 1:3)
head(tidd)
## class simv
## 1 V1 0.5482568
## 2 V1 0.4950737
## 3 V1 0.5956448
## 4 V1 0.4132343
## 5 V1 0.4715408
## 6 V1 0.5087496
library(dplyr)
tidd %>%
group_by(class)%>%
summarise(mean = mean(simv), sd = sd(simv))
## # A tibble: 3 x 3
## class mean sd
## <chr> <dbl> <dbl>
## 1 V1 0.522 0.0496
## 2 V2 0.750 0.0414
## 3 V3 0.255 0.0463
library(ggplot2)
ggplot(tidd, aes(simv)) + geom_histogram() + facet_grid(~class)
So I reimplemented the computations described on the lecture slides. I had a few issues with the library and getting the data to show but I realized I had to run some other packages first in order for the computations to work.