. do "C:\Users\hstryhn\AppData\Local\Temp\STDdb8_000000.tmp"

. * do-file for lecture 3b of VHM 812/802, Winter 2022
. * two additions: invlogit function, and code with glm command for grouped data
. version 17 /* works also with versions 14-16 */

. set more off

. cd "r:\"
r:\

.
. use micegrp.dta, clear

. regress dead dose /* to illustrate extrapolation below zero in regression */

      Source |       SS           df       MS      Number of obs   =        12
-------------+----------------------------------   F(1, 10)        =     29.83
       Model |   49.616367         1   49.616367   Prob > F        =    0.0003
    Residual |   16.633633        10   1.6633633   R-squared       =    0.7489
-------------+----------------------------------   Adj R-squared   =    0.7238
       Total |       66.25        11  6.02272727   Root MSE        =    1.2897

-------------+----------------------------------------------------------------
         dead | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
        dose |     29.1441    5.336195      5.46    0.000      17.25431     41.03388
       _cons |   -2.107171    1.049163     -2.01    0.072     -4.444851     .2305092
-------------+----------------------------------------------------------------

. use mice.dta, clear

. regress dead dose /* same estimates of regression coefficients, up to scaling */

      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 | Coefficient   Std. err.       t    P>|t|      [95% con f. 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 | Coefficient   Std. err.       z    P>|z|      [95% con f. 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 | Coefficient   Std. err.       z    P>|z|      [95% con f. 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 */

Goodness-of-fit test after logistic model
Variable: dead

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
-------------+----------------------------------------------------------------
             |       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 H0:   .0855132     4.87   0.000
-------------+----------------------------------------------------------------
        diff = prop(x) - prop(y)                                  z =   4.8725
    H0: 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 i.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 | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
     1.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 i.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% con f. interval]
-------------+----------------------------------------------------------------
     1.dose2 |    8.653061    4.138022      4.51    0.000       3.38932     22.09159
       _cons |    .1320755    .0531142     -5.03    0.000      .0600497     .2904914
-------------+----------------------------------------------------------------


Note: _cons estimates baseline odds.


. display 1/(1+exp(2.024382)) /* baseline prob */
.11666664

. display invlogit(-2.024382) /* same with invlogit function */
.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 | Coefficient   Std. err.       z    P>|z|      [95% con f. 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% con f. interval]
-------------+----------------------------------------------------------------
        dose |     2273668     7577074      4.39    0.000      3311.843     1.56e+09
       _cons |    .0281632    .0198647     -5.06    0.000      .0070679     .1122204
-------------+----------------------------------------------------------------


Note: _cons estimates baseline odds.


. display 1/(1+exp(3.569738)) /* intercept prob */
.02739179

. display invlogit(-3.569738) /* same with invlogit function */
.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 | Coefficient   Std. err.       z    P>|z|      [95% con f. 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
-------------+----------------------------------------------------------------

. * see below for equivalent code with glm command (blogit is outdated)
. predict linpred, xb

. predict phat, p

. predict deadhat, n

. * illustration of calculation of X^2
. generate x2=(dead-deadhat)^2/(total*phat*(1-phat))

. total x2

Total estimation                            Number of obs    = 12

-------------+------------------------------------------------
             |      Total    Std. err.      [95% con f. interval]
-------------+------------------------------------------------
          x2 |    8.742145    3.243557       1.603125     15.88116
-------------+------------------------------------------------

. display chi2tail(12-1-1,8.742145)
.55673739

. * using instead glm command
. capture drop linpred phat deadhat

. glm dead dose, fam(binomial total)

Iteration 0:   log likelihood =  -18.60686
Iteration 1:   log likelihood = -18.516941
Iteration 2:   log likelihood = -18.516789
Iteration 3:   log likelihood = -18.516789

Generalized linear models                         Number of obs   =         12
Optimization     : ML                             Residual df     =         10
                                                  Scale parameter =          1
Deviance         =  10.25371404                   (1/df) Deviance =   1.025371
Pearson          =    8.7421447                   (1/df) Pearson  =   .8742145

Variance function: V(u) =  u*(1-u/total)          [Binomial]
Link function    : g(u) =  ln(u/(total-u))        [Logit]

                                                  AIC             =   3.419465
Log likelihood   = -18.51678904                   BIC             =  -14.59535

-------------+----------------------------------------------------------------
             |                 OIM
         dead | Coefficient   std. err.       z    P>|z|      [95% con f. 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 deadhat, mu

. generate phat=deadhat/total

. * Pearson goodness-of-fit test
. display chi2tail(12-1-1,8.7421447)
.55673742

. display chi2tail(12-1-1,e(deviance_p))
.55673742

.
end of do-file

. exit, clear