* do-file for lecture 3b of VHM 802/812, Winter 2016 . version 14 /* works also with version 13 */ . set more off . cd "h:\vhm\vhm802\data_stata" h:\vhm\vhm802\data_stata . . use mice.dta, clear . regress dead dose /* to illustrate extrapolation below zero in regression */ Source | SS df MS Number of obs = 120 -------------+---------------------------------- F(1, 118) = 27.41 Model | 4.9616367 1 4.9616367 Prob > F = 0.0000 Residual | 21.3633633 118 .181045452 R-squared = 0.1885 -------------+---------------------------------- Adj R-squared = 0.1816 Total | 26.325 119 .221218487 Root MSE = .42549 ------------------------------------------------------------------------------ dead | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- dose | 2.91441 .5567137 5.24 0.000 1.811965 4.016855 _cons | -.2107171 .1094569 -1.93 0.057 -.4274715 .0060373 ------------------------------------------------------------------------------ . . * logistic regression using logit command . logit dead dose Iteration 0: log likelihood = -75.669723 Iteration 1: log likelihood = -64.201632 Iteration 2: log likelihood = -63.944944 Iteration 3: log likelihood = -63.944713 Iteration 4: log likelihood = -63.944713 Logistic regression Number of obs = 120 LR chi2(1) = 23.45 Prob > chi2 = 0.0000 Log likelihood = -63.944713 Pseudo R2 = 0.1549 ------------------------------------------------------------------------------ dead | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- dose | 14.6369 3.332533 4.39 0.000 8.10526 21.16855 _cons | -3.569738 .7053432 -5.06 0.000 -4.952185 -2.187291 ------------------------------------------------------------------------------ . predict linpred, xb /* predictions on logit scale */ . twoway (scatter linpred dose) (line linpred dose) . lincom _cons+dose*0.1 ( 1) .1*[dead]dose + [dead]_cons = 0 ------------------------------------------------------------------------------ dead | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- (1) | -2.106048 .4015256 -5.25 0.000 -2.893023 -1.319072 ------------------------------------------------------------------------------ . predict phat, p /* predicted probabilities */ . twoway (scatter phat dose) (line phat dose) . estat gof /* Pearson goodness-of-fit test */ Logistic model for dead, goodness-of-fit test number of observations = 120 number of covariate patterns = 12 Pearson chi2(10) = 8.74 Prob > chi2 = 0.5567 . . * analysis with dichotomous dose . prtesti 60 32 60 7, count /* 2 independent proportions */ Two-sample test of proportions x: Number of obs = 60 y: Number of obs = 60 ------------------------------------------------------------------------------ Variable | Mean Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- x | .5333333 .0644061 .4070997 .659567 y | .1166667 .0414438 .0354382 .1978951 -------------+---------------------------------------------------------------- diff | .4166667 .0765881 .2665567 .5667766 | under Ho: .0855132 4.87 0.000 ------------------------------------------------------------------------------ diff = prop(x) - prop(y) z = 4.8725 Ho: diff = 0 Ha: diff < 0 Ha: diff != 0 Ha: diff > 0 Pr(Z < z) = 1.0000 Pr(|Z| < |z|) = 0.0000 Pr(Z > z) = 0.0000 . csi 32 7 28 53, or /* relative risk and odds-ratio */ | Exposed Unexposed | Total -----------------+------------------------+------------ Cases | 32 7 | 39 Noncases | 28 53 | 81 -----------------+------------------------+------------ Total | 60 60 | 120 | | Risk | .5333333 .1166667 | .325 | | | Point estimate | [95% Conf. Interval] |------------------------+------------------------ Risk difference | .4166667 | .2665567 .5667766 Risk ratio | 4.571429 | 2.191203 9.53721 Attr. frac. ex. | .78125 | .5436296 .8951475 Attr. frac. pop | .6410256 | Odds ratio | 8.653061 | 3.441938 21.63764 (Cornfield) +------------------------------------------------- chi2(1) = 23.74 Pr>chi2 = 0.0000 . generate dose2=(dose>0.16) if dose~=. . tabulate dead dose2, chi2 column exact lrchi2 /* 2x2 table analysis */ +-------------------+ | Key | |-------------------| | frequency | | column percentage | +-------------------+ | dose2 dead | 0 1 | Total -----------+----------------------+---------- 0 | 53 28 | 81 | 88.33 46.67 | 67.50 -----------+----------------------+---------- 1 | 7 32 | 39 | 11.67 53.33 | 32.50 -----------+----------------------+---------- Total | 60 60 | 120 | 100.00 100.00 | 100.00 Pearson chi2(1) = 23.7417 Pr = 0.000 likelihood-ratio chi2(1) = 25.2010 Pr = 0.000 Fisher's exact = 0.000 1-sided Fisher's exact = 0.000 . logit dead dose2 /* logistic regression */ Iteration 0: log likelihood = -75.669723 Iteration 1: log likelihood = -63.522292 Iteration 2: log likelihood = -63.070143 Iteration 3: log likelihood = -63.06923 Iteration 4: log likelihood = -63.06923 Logistic regression Number of obs = 120 LR chi2(1) = 25.20 Prob > chi2 = 0.0000 Log likelihood = -63.06923 Pseudo R2 = 0.1665 ------------------------------------------------------------------------------ dead | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- dose2 | 2.157913 .4782148 4.51 0.000 1.220629 3.095197 _cons | -2.024382 .4021506 -5.03 0.000 -2.812582 -1.236181 ------------------------------------------------------------------------------ . logistic dead dose2 /* same as logit with option or */ Logistic regression Number of obs = 120 LR chi2(1) = 25.20 Prob > chi2 = 0.0000 Log likelihood = -63.06923 Pseudo R2 = 0.1665 ------------------------------------------------------------------------------ dead | Odds Ratio Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- dose2 | 8.653061 4.138022 4.51 0.000 3.38932 22.09159 _cons | .1320755 .0531142 -5.03 0.000 .0600497 .2904914 ------------------------------------------------------------------------------ . display 1/(1+exp(2.024382)) /* baseline prob */ .11666664 . logit dead dose Iteration 0: log likelihood = -75.669723 Iteration 1: log likelihood = -64.201632 Iteration 2: log likelihood = -63.944944 Iteration 3: log likelihood = -63.944713 Iteration 4: log likelihood = -63.944713 Logistic regression Number of obs = 120 LR chi2(1) = 23.45 Prob > chi2 = 0.0000 Log likelihood = -63.944713 Pseudo R2 = 0.1549 ------------------------------------------------------------------------------ dead | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- dose | 14.6369 3.332533 4.39 0.000 8.10526 21.16855 _cons | -3.569738 .7053432 -5.06 0.000 -4.952185 -2.187291 ------------------------------------------------------------------------------ . logistic dead dose /* same as logit with option or */ Logistic regression Number of obs = 120 LR chi2(1) = 23.45 Prob > chi2 = 0.0000 Log likelihood = -63.944713 Pseudo R2 = 0.1549 ------------------------------------------------------------------------------ dead | Odds Ratio Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- dose | 2273668 7577074 4.39 0.000 3311.843 1.56e+09 _cons | .0281632 .0198647 -5.06 0.000 .0070679 .1122204 ------------------------------------------------------------------------------ . display 1/(1+exp(3.569738)) /* intercept prob */ .02739179 . . * demonstration of analysis using grouped data . use micegrp.dta, clear . blogit dead total dose Logistic regression for grouped data Number of obs = 120 LR chi2(1) = 23.45 Prob > chi2 = 0.0000 Log likelihood = -63.944713 Pseudo R2 = 0.1549 ------------------------------------------------------------------------------ _outcome | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- dose | 14.6369 3.332533 4.39 0.000 8.10526 21.16855 _cons | -3.569738 .7053432 -5.06 0.000 -4.952185 -2.187291 ------------------------------------------------------------------------------ . predict linpred, xb . predict phat, p . predict deadhat, n