. * do-file for lecture 10 of VHM 802, Winter 2016 . version 14 /* works also with version 13 */ . set more off . cd "h:\vhm\vhm802" h:\vhm\vhm802 . . * Example 16.5 . import delimited data_csv\ch16ta1.csv, clear (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 | Coef. Std. Err. t P>|t| [95% Conf. 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 7 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 . end of do-file . do "C:\Users\DEFAUL~1.SID\AppData\Local\Temp\STD00000000.tmp" . 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 | Coef. Std. Err. z P>|z| [95% Conf. 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 | 8.88e-15 1.702395 0.00 1.000 -3.336632 3.336632 2 4 | 1.07e-14 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.108 -----------------------------+------------------------------------------------ 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 data_csv\ch16ta2.csv, clear (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 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 ------------------------------------------------------------- | 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 | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- 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 | 84.3375 5.147868 16.38 0.000 74.24786 94.42713 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ tray: Identity | var(_cons) | 14.4469 12.66438 2.59178 80.5288 -----------------------------+------------------------------------------------ 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% Conf. 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 ------------------------------------------------------------------------------ . * 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