. * do-file for lecture 10 of VHM 802, Winter 2014 . version 13 . set more off . cd "f:\vhm\vhm802\data_csv" f:\vhm\vhm802\data_csv . . * Example 16.5 . insheet using 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.83333 23 52.4275362 24.12 0.0000 | anxiety | 10.0833333 1 10.0833333 0.98 0.3517 tension | 8.33333333 1 8.33333333 0.81 0.3949 anxiety#tension | 80.0833333 1 80.0833333 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.41666667 3 2.80555556 1.29 0.3003 tension#type | 12.1666667 3 4.05555556 1.87 0.1624 anxiety#tension#type | 12.75 3 4.25 1.96 0.1477 | Residual | 52.1666667 24 2.17361111 ------------------------+---------------------------------------------------- Total | 1258 47 26.7659574 . 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.20833333 3.18 0.0845 | anxiety | 2.52083333 1 2.52083333 0.98 0.3517 tension | 2.08333333 1 2.08333333 0.81 0.3949 anxiety#tension | 20.0208333 1 20.0208333 7.77 0.0237 | Residual | 20.625 8 2.578125 ----------------+---------------------------------------------------- Total | 45.25 11 4.11363636 . 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" . * 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 | 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 regression: 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 . insheet using 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.7465 35 331.507043 310.30 0.0000 | n | 3197.05502 3 1065.68501 11.46 0.0377 table | 14.3008467 1 14.3008467 0.15 0.7211 n#table | 278.950779 3 92.9835932 -------------+---------------------------------------------------- weed | 7001.25534 2 3500.62767 555.45 0.0000 n#weed | 929.516211 6 154.919368 24.58 0.0001 n#table#weed | 50.4183374 8 6.30229218 -------------+---------------------------------------------------- clip | 125.453328 1 125.453328 117.43 0.0000 n#clip | .734998995 3 .244999665 0.23 0.8742 weed#clip | .245415154 2 .122707577 0.11 0.8925 n#weed#clip | 4.81624219 6 .802707032 0.75 0.6203 | Residual | 12.8199971 12 1.06833309 -------------+---------------------------------------------------- Total | 11615.5665 47 247.139713 . * 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 regression: chi2(2) = 32.98 Prob > chi2 = 0.0000 Note: LR test is conservative and provided only for reference. . * residuals 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.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