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

. * do-file for lecture 11 of VHM 802, Winter 2025
. version 18 /* works also with versions 14-17 */

. set more off

. set scheme stcolor_alt

. cd "r:\"
r:\

.
. * Example 16.5
. import delimited ch16ta1.csv, clear
(encoding automatically selected: ISO-8859-1)
(5 vars, 48 obs)

. anova errors anxiety##tension / subject|anxiety#tension anxiety##tension##type

                         Number of obs =         48    R-squared     =  0.9585
                         Root MSE      =    1.47432    Adj R-squared =  0.9188

                  Source | Partial SS         df         MS        F    Prob>F
 ------------------------+----------------------------------------------------
                   Model |  1205.8333         23   52.427536     24.12  0.0000
                         |
                 anxiety |  10.083333          1   10.083333      0.98  0.3517
                 tension |  8.3333333          1   8.3333333      0.81  0.3949
         anxiety#tension |  80.083333          1   80.083333      7.77  0.0237
 subject|anxiety#tension |       82.5          8     10.3125
 ------------------------+----------------------------------------------------
                    type |      991.5          3       330.5    152.05  0.0000
            anxiety#type |  8.4166667          3   2.8055556      1.29  0.3003
            tension#type |  12.166667          3   4.0555556      1.87  0.1624
    anxiety#tension#type |      12.75          3        4.25      1.96  0.1477
                         |
                Residual |  52.166667         24   2.1736111
 ------------------------+----------------------------------------------------
                   Total |       1258         47   26.765957

. regress

      Source |       SS           df       MS      Number of obs   =        48
-------------+----------------------------------   F(23, 24)       =     24.12
       Model |  1205.83333        23  52.4275362   Prob > F        =    0.0000
    Residual |  52.1666667        24  2.17361111   R-squared       =    0.9585
-------------+----------------------------------   Adj R-squared   =    0.9188
       Total |        1258        47  26.7659574   Root MSE        =    1.4743

------------------------+----------------------------------------------------------------
                  errors | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
------------------------+----------------------------------------------------------------
              2.anxiety |   -1.666667    1.474317     -1.13    0.269     -4.709508     1.376174
              2.tension |   -1.916667    1.474317     -1.30    0.206     -4.959508     1.126174
                        |
        anxiety#tension |
                   2 2  |    2.583333    2.084999      1.24    0.227     -1.719894      6.88656
                        |
subject|anxiety#tension |
                 2 1 1  |       -1.75      1.0425     -1.68    0.106     -3.901614     .4016136
                 2 1 2  |        -3.5      1.0425     -3.36    0.003     -5.651614    -1.348386
                 2 2 1  |       -1.75      1.0425     -1.68    0.106     -3.901614     .4016136
                 2 2 2  |          -1      1.0425     -0.96    0.347     -3.151614     1.151614
                 3 1 1  |        -4.5      1.0425     -4.32    0.000     -6.651614    -2.348386
                 3 1 2  |          -2      1.0425     -1.92    0.067     -4.151614     .1516136
                 3 2 1  |         -.5      1.0425     -0.48    0.636     -2.651614     1.651614
                 3 2 2  |       -2.25      1.0425     -2.16    0.041     -4.401614    -.0983864
                        |
                   type |
                     2  |          -5    1.203775     -4.15    0.000     -7.484469    -2.515531
                     3  |   -8.333333    1.203775     -6.92    0.000      -10.8178    -5.848864
                     4  |         -13    1.203775    -10.80    0.000     -15.48447    -10.51553
                        |
           anxiety#type |
                   2 2  |   -1.666667    1.702395     -0.98    0.337     -5.180237     1.846904
                   2 3  |   -2.333333    1.702395     -1.37    0.183     -5.846904     1.180237
                   2 4  |   -1.333333    1.702395     -0.78    0.441     -4.846904     2.180237
                        |
           tension#type |
                   2 2  |   -.3333333    1.702395     -0.20    0.846     -3.846904     3.180237
                   2 3  |   -1.18e-15    1.702395     -0.00    1.000      -3.51357      3.51357
                   2 4  |   -1.24e-14    1.702395     -0.00    1.000      -3.51357      3.51357
                        |
   anxiety#tension#type |
                 2 2 2  |           4     2.40755      1.66    0.110     -.9689387     8.968939
                 2 2 3  |           3     2.40755      1.25    0.225     -1.968939     7.968939
                 2 2 4  |    5.666667     2.40755      2.35    0.027       .697728     10.63561
                        |
                  _cons |    19.08333      1.0425     18.31    0.000      16.93172     21.23495
------------------------+----------------------------------------------------------------

. * residual analysis
. predict pred, xb

. predict stdres, rstandard

. scatter stdres pred

. qnorm stdres

. swilk stdres

                   Shapiro–Wilk W test for normal data

    Variable |        Obs       W           V         z       Prob>z
-------------+------------------------------------------------------
      stdres |         48    0.97244      1.255     0.484    0.31434

.
. * random effect checks
. preserve /* run next 8 lines together with this one */

. collapse (mean) errors, by(anxiety tension subject)

. anova errors anxiety##tension

                         Number of obs =         12    R-squared     =  0.5442
                         Root MSE      =    1.60565    Adj R-squared =  0.3733

                  Source | Partial SS         df         MS        F    Prob>F
         ----------------+----------------------------------------------------
                   Model |     24.625          3   8.2083333      3.18  0.0845
                         |
                 anxiety |  2.5208333          1   2.5208333      0.98  0.3517
                 tension |  2.0833333          1   2.0833333      0.81  0.3949
         anxiety#tension |  20.020833          1   20.020833      7.77  0.0237
                         |
                Residual |     20.625          8    2.578125
         ----------------+----------------------------------------------------
                   Total |      45.25         11   4.1136364

. predict pred2, xb

. predict stdres2, rstandard

. scatter stdres2 pred2

. qnorm stdres2

. swilk stdres2

                   Shapiro–Wilk W test for normal data

    Variable |        Obs       W           V         z       Prob>z
-------------+------------------------------------------------------
     stdres2 |         12    0.97505      0.417    -1.705    0.95589

. restore

.
. * likelihood-based analysis
. egen id=group(anxiety tension subject)

. mixed errors anxiety##tension##type || id:, reml

Performing EM optimization ...

Performing gradient-based optimization:
Iteration 0:  Log restricted-likelihood = -72.845037
Iteration 1:  Log restricted-likelihood = -72.845037

Computing standard errors ...

Mixed-effects REML regression                        Number of obs    =     48
Group variable: id                                   Number of groups  =     12
                                                     Obs per group:
                                                                  min =      4
                                                                  avg =    4.0
                                                                  max =      4
                                                     Wald chi2(15)    = 481.04
Log restricted-likelihood = -72.845037               Prob > chi2      = 0.0000

---------------------+----------------------------------------------------------------
               errors | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
---------------------+----------------------------------------------------------------
           2.anxiety |   -.3333333    1.674979     -0.20    0.842     -3.616233     2.949566
           2.tension |   -1.666667    1.674979     -1.00    0.320     -4.949566     1.616233
                     |
     anxiety#tension |
                2 2  |           2    2.368779      0.84    0.398     -2.642721     6.642721
                     |
                type |
                  2  |          -5    1.203775     -4.15    0.000     -7.359355    -2.640645
                  3  |   -8.333333    1.203775     -6.92    0.000     -10.69269    -5.973978
                  4  |         -13    1.203775    -10.80    0.000     -15.35936    -10.64064
                     |
        anxiety#type |
                2 2  |   -1.666667    1.702395     -0.98    0.328     -5.003299     1.669966
                2 3  |   -2.333333    1.702395     -1.37    0.170     -5.669966     1.003299
                2 4  |   -1.333333    1.702395     -0.78    0.434     -4.669966     2.003299
                     |
        tension#type |
                2 2  |   -.3333333    1.702395     -0.20    0.845     -3.669966     3.003299
                2 3  |    2.52e-15    1.702395      0.00    1.000     -3.336632     3.336632
                2 4  |           0    1.702395      0.00    1.000     -3.336632     3.336632
                     |
anxiety#tension#type |
              2 2 2  |           4     2.40755      1.66    0.097     -.7187106     8.718711
              2 2 3  |           3     2.40755      1.25    0.213     -1.718711     7.718711
              2 2 4  |    5.666667     2.40755      2.35    0.019       .947956     10.38538
                     |
               _cons |          17    1.184389     14.35    0.000      14.67864     19.32136
---------------------+----------------------------------------------------------------

-----------------------------+------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
id: Identity                 |
                  var(_cons) |   2.034724   1.298573      .5824564    7.107999
-----------------------------+------------------------------------------------
               var(Residual) |   2.173611   .6274673      1.234415    3.827388
-----------------------------+------------------------------------------------
LR test vs. linear model: chibar2(01) = 8.69          Prob >= chibar2 =  0.0016

. * residual analysis
. predict pred2, fitted

. predict stdres2, rstandard

. scatter stdres2 pred2

. qnorm stdres2

. predict pred_ref, reffects

. bysort id: gen within_n=_n

. qnorm pred_ref if within_n==1

. swilk pred_ref if within_n==1

                   Shapiro–Wilk W test for normal data

    Variable |        Obs       W           V         z       Prob>z
-------------+------------------------------------------------------
    pred_ref |         12    0.97505      0.417    -1.705    0.95589

.
. * Example 16.7
. import delimited ch16ta2.csv, clear
(encoding automatically selected: ISO-8859-9)
(7 vars, 48 obs)

. anova biomass n table / n#table n##weed / n#table#weed n##weed##clip

                         Number of obs =         48    R-squared     =  0.9989
                         Root MSE      =     1.0336    Adj R-squared =  0.9957

                  Source | Partial SS         df         MS        F    Prob>F
            -------------+----------------------------------------------------
                   Model |  11602.747         35   331.50704    310.30  0.0000
                         |
                       n |   3197.055          3    1065.685     11.46  0.0377
                   table |  14.300847          1   14.300847      0.15  0.7211
                 n#table |  278.95078          3   92.983593
            -------------+----------------------------------------------------
                    weed |  7001.2553          2   3500.6277    555.45  0.0000
                  n#weed |  929.51621          6   154.91937     24.58  0.0001
            n#table#weed |  50.418337          8   6.3022922
            -------------+----------------------------------------------------
                    clip |  125.45333          1   125.45333    117.43  0.0000
                  n#clip |    .734999          3   .24499967      0.23  0.8742
               weed#clip |  .24541515          2   .12270758      0.11  0.8925
             n#weed#clip |  4.8162422          6   .80270703      0.75  0.6203
                         |
                Residual |  12.819997         12   1.0683331
            -------------+----------------------------------------------------
                   Total |  11615.567         47   247.13971

. predict pred, xb

. predict stdres, rstandard

. scatter stdres pred

. qnorm stdres

. * computation of means for model checking:
. * collapse (mean) biomass, by(table n weed)
. * collapse (mean) biomass, by(table n)
.
. * likelihood-based analysis
. mixed biomass i.table n##weed##clip || tray: || wetland:, reml

Performing EM optimization ...

Performing gradient-based optimization:
Iteration 0:  Log restricted-likelihood = -56.754758
Iteration 1:  Log restricted-likelihood = -56.754691
Iteration 2:  Log restricted-likelihood = -56.754691

Computing standard errors ...

Mixed-effects REML regression                          Number of obs  =      48

        Grouping information
        ----------------+--------------------------------------------
                        |     No. of       Observations per group
         Group variable |     groups    Minimum    Average    Maximum
        ----------------+--------------------------------------------
                   tray |          8          6        6.0          6
                wetland |         24          2        2.0          2
        ----------------+--------------------------------------------

                                                       Wald chi2(24)  = 1415.79
Log restricted-likelihood = -56.754691                 Prob > chi2   =  0.0000

-------------+----------------------------------------------------------------
      biomass | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
     2.table |   -1.091667    2.783638     -0.39    0.695     -6.547498     4.364163
             |
           n |
          2  |   -2.549995    4.258194     -0.60    0.549      -10.8959     5.795911
          3  |   -3.099995    4.258194     -0.73    0.467      -11.4459     5.245911
          4  |   -6.099998    4.258194     -1.43    0.152      -14.4459     2.245908
             |
        weed |
          2  |      -14.95    1.919717     -7.79    0.000     -18.71257    -11.18742
          3  |   -12.09999    1.919717     -6.30    0.000     -15.86257    -8.337419
             |
      n#weed |
        2 2  |   -6.800005     2.71489     -2.50    0.012     -12.12109    -1.478919
        2 3  |   -9.050007     2.71489     -3.33    0.001     -14.37109    -3.728921
        3 2  |   -15.50001     2.71489     -5.71    0.000     -20.82109    -10.17892
        3 3  |   -17.50001     2.71489     -6.45    0.000     -22.82109    -12.17892
        4 2  |       -24.3     2.71489     -8.95    0.000     -29.62109    -18.97892
        4 3  |       -24.4     2.71489     -8.99    0.000     -29.72109    -19.07892
             |
      2.clip |    1.950005    1.033602      1.89    0.059     -.0758178     3.975827
             |
      n#clip |
        2 2  |    2.399994    1.461734      1.64    0.101     -.4649515     5.264939
        3 2  |    .3999939    1.461734      0.27    0.784     -2.464952     3.264939
        4 2  |    1.549995    1.461734      1.06    0.289      -1.31495     4.414941
             |
   weed#clip |
        2 2  |    2.299995    1.461734      1.57    0.116       -.56495     5.164941
        3 2  |    1.049992    1.461734      0.72    0.473     -1.814954     3.914937
             |
 n#weed#clip |
      2 2 2  |   -3.599993    2.067204     -1.74    0.082     -7.651637      .451652
      2 3 2  |   -2.349991    2.067204     -1.14    0.256     -6.401636     1.701654
      3 2 2  |   -2.049994    2.067204     -0.99    0.321     -6.101638     2.001651
      3 3 2  |    .4000111    2.067204      0.19    0.847     -3.651634     4.451656
      4 2 2  |   -2.549995    2.067204     -1.23    0.217      -6.60164     1.501649
      4 3 2  |    -.899992    2.067204     -0.44    0.663     -4.951637     3.151653
             |
       _cons |    83.24583    3.317117     25.10    0.000       76.7444     89.74726
-------------+----------------------------------------------------------------

-----------------------------+------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
tray: Identity               |
                  var(_cons) |    14.4469   12.66438       2.59178    80.52881
-----------------------------+------------------------------------------------
wetland: Identity            |
                  var(_cons) |    2.61698   1.590593      .7951379    8.613078
-----------------------------+------------------------------------------------
               var(Residual) |   1.068333   .4361449        .47996    2.377979
-----------------------------+------------------------------------------------
LR test vs. linear model: chi2(2) = 32.98                 Prob > chi2 = 0.0000


Note: LR test is conservative and provided only for reference.


. testparm n#weed#clip


 ( 1)  [biomass]2.n#2.weed#2.clip = 0


 ( 2)  [biomass]2.n#3.weed#2.clip = 0


 ( 3)  [biomass]3.n#2.weed#2.clip = 0


 ( 4)  [biomass]3.n#3.weed#2.clip = 0


 ( 5)  [biomass]4.n#2.weed#2.clip = 0


 ( 6)  [biomass]4.n#3.weed#2.clip = 0


           chi2(  6) =    4.51
         Prob > chi2 =    0.6082

. margins n#weed clip

 Predictive margins                                         Number of obs    = 48


Expression:  Linear prediction, fixed portion, predict()


-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
      n#weed |
        1 1  |      83.675    2.966315     28.21    0.000      77.86113     89.48887
        1 2  |      69.875    2.966315     23.56    0.000      64.06113     75.68887
        1 3  |        72.1    2.966315     24.31    0.000      66.28613     77.91387
        2 1  |      82.325    2.966315     27.75    0.000      76.51113     88.13887
        2 2  |      59.925    2.966315     20.20    0.000      54.11113     65.73887
        2 3  |      60.525    2.966315     20.40    0.000      54.71113     66.33887
        3 1  |      80.775    2.966315     27.23    0.000      74.96113     86.58887
        3 2  |       50.45    2.966315     17.01    0.000      44.63613     56.26387
        3 3  |        51.9    2.966315     17.50    0.000      46.08613     57.71387
        4 1  |       78.35    2.966315     26.41    0.000      72.53613     84.16387
        4 2  |      38.975    2.966315     13.14    0.000      33.16113     44.78887
        4 3  |      41.925    2.966315     14.13    0.000      36.11113     47.73887
             |
        clip |
          1  |    62.61667    1.399792     44.73    0.000      59.87313     65.36021
          2  |       65.85    1.399792     47.04    0.000      63.10646     68.59354
-------------+----------------------------------------------------------------

. margins n#weed

 Predictive margins                                         Number of obs    = 48


Expression:  Linear prediction, fixed portion, predict()


-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
      n#weed |
        1 1  |      83.675    2.966315     28.21    0.000      77.86113     89.48887
        1 2  |      69.875    2.966315     23.56    0.000      64.06113     75.68887
        1 3  |        72.1    2.966315     24.31    0.000      66.28613     77.91387
        2 1  |      82.325    2.966315     27.75    0.000      76.51113     88.13887
        2 2  |      59.925    2.966315     20.20    0.000      54.11113     65.73887
        2 3  |      60.525    2.966315     20.40    0.000      54.71113     66.33887
        3 1  |      80.775    2.966315     27.23    0.000      74.96113     86.58887
        3 2  |       50.45    2.966315     17.01    0.000      44.63613     56.26387
        3 3  |        51.9    2.966315     17.50    0.000      46.08613     57.71387
        4 1  |       78.35    2.966315     26.41    0.000      72.53613     84.16387
        4 2  |      38.975    2.966315     13.14    0.000      33.16113     44.78887
        4 3  |      41.925    2.966315     14.13    0.000      36.11113     47.73887
-------------+----------------------------------------------------------------

. marginsplot, noci


Variables that uniquely identify margins: n weed


. * residuals analysis
. capture drop pred stdres

. predict pred, xb

. predict stdres, rstandard

. scatter stdres pred

. qnorm stdres

. swilk stdres

                   Shapiro–Wilk W test for normal data

    Variable |        Obs       W           V         z       Prob>z
-------------+------------------------------------------------------
      stdres |         48    0.98467      0.698    -0.764    0.77756

. predict pred_ref*, reffects

. bysort tray: gen within1_n=_n

. bysort wetland: gen within2_n=_n

. qnorm pred_ref1 if within1_n==1

. swilk pred_ref1 if within1_n==1

                   Shapiro–Wilk W test for normal data

    Variable |        Obs       W           V         z       Prob>z
-------------+------------------------------------------------------
   pred_ref1 |          8    0.99541      0.064    -3.373    0.99963

. qnorm pred_ref2 if within2_n==1

. swilk pred_ref2 if within2_n==1

                   Shapiro–Wilk W test for normal data

    Variable |        Obs       W           V         z       Prob>z
-------------+------------------------------------------------------
   pred_ref2 |         24    0.96355      0.983    -0.035    0.51382

.
end of do-file

. exit, clear