*******************************************************************************
*** ***
*** Stata-program (do-file) associated with the Exercises on day one. ***
*** ***
*** ***
*******************************************************************************
* Specifying the path to where the data are, for example
cd "C:\ANOVA\DAY 1"
********** Exercise 1 (red cell folate levels) *************
* Reading in the data.
use folate.dta, clear
* 1. Calculate relevant summary statistics using the tabstat command. Make a
* scatter plot of the data. What would you expect to find in a statistical
* analysis of the data?
tabstat folate, by(group) stat(n mean sd)
* group | N mean sd
* ---------+------------------------------
* 1 | 9 315.7778 54.98358
* 2 | 8 249.875 33.63008
* 3 | 5 278 33.75648
* ---------+------------------------------
* Total | 22 283.2273 51.28439
dotplot folate, over(group) center
* Larger levels of red cell folate in group 1 followed by 3 and 2. Slightly
* larger variation in group 3. We would expect that the mean levels in all three
* groups are statistically significantly different from each other but that the
* variation is more or less the same in the three groups.
* 2. Test the hypothesis of equal variations in the three groups using the
* oneway command. Test for equal levels. How would you formulate the
* conclusions of the tests?
oneway folate group
* We conclude from the Bartlett-test that there is little evidence against the
* hypothesis of equal standard variations in the three groups (p=0.35). We also
* see in the output from oneway that we just reject the hypothesis of equal red
* cell folate levels in the three groups (p=0.02).
* 3. Calculate estimates and 95%-confidence intervals for the expected red cell
* folate levels using the margins command. Make a plot of the estimates with
* 95%-confidence intervals using the marginsplot command. Why are the
* intervals not of equal width?
anova folate group
margins, over(group)
* | Margin Std. Err. t P>|t| [95% Conf. Interval]
* ------------+----------------------------------------------------------------
* group |
* 1 | 315.7778 14.64201 21.57 0.000 285.1317 346.4239
* 2 | 249.875 15.5302 16.09 0.000 217.3699 282.3801
* 3 | 278 19.64432 14.15 0.000 236.884 319.116
marginsplot, recast(scatter)
* The increasing width of the 95%-confidence intervals with increasing group
* number reflects the decreasing number of observations in the groups. The group
* means are estimated with a greater precision the more observations you base
* the estimates on.
* 4. Make pairwise comparisons of the three groups using the pwcompare command.
* Do any two groups differ from each other?
pwcompare group, eff
* | Contrast Std. Err. t P>|t| [95% Conf. Interval]
* ------------+----------------------------------------------------------------
* group |
* 2 vs 1 | -65.90278 21.34422 -3.09 0.006 -110.5767 -21.22882
* 3 vs 1 | -37.77778 24.50077 -1.54 0.140 -89.05848 13.50292
* 3 vs 2 | 28.125 25.04169 1.12 0.275 -24.28786 80.53786
* Of course the groups differ from each other (they received different
* treatments). Also if you just consider the mean red cell folate level, it is
* our best guess that the three groups differ as the means are not identical.
* If the question is whether the mean red cell folate level are significantly
* different in two groups, then we would indeed conclude that they differ
* significantly in groups 1 and 2. The mean red cell folate level in group 3 is
* not signifficantly different from that in the other two groups.
* 5. Check the assumptions behind the model by plotting the residuals against
* percentiles from the normal distribution. Does the plot give you reason to
* doubt the appropriateness of the model to describe the data?
predict res, residual
qnorm res
* Looking at the QQ-plot there are no systematic deviations from a straight
* line. This means that there is not anything in this part of the model
* validation that leads us to reject the statistical model.
* 6. Write a short summary of statistical methods used in the analysis and the
* findings.
* The red cell folate levels were analyzed using a oneway ANOVA. It was accepted
* that the standard deviations were the same in the three groups (p=0.35), and
* the QQ-plot did not give rise to concern about the models ability to describe
* the data adequately. Overall there was evidence against the hypothesis of
* equal mean red cell folate levels in the three groups (p=0.02). The estimated
* mean levels were 315.8 [285.1 ; 346.4] ng/ml in group 1, 249.9 [217.4 ; 282.4]
* ng/ml in group 2, and 278.0 [236.9 ; 319.1] ng/ml in group 3. The mean level
* was significantly higher in group 1 compared to group 2 (65.9 ng/ml, 95%-CI:
* 21.2 - 110.6, p=0.006). The mean level in group 3 did not differ significantly
* from those in group 1 and 2 (1-3: 37.8 ng/ml, 95%-CI: -13.5 - 89.1, p=0.14,
* 3-2: 28.1 ng/ml, 95%-CI: -24.3 - 80.5, p=0.28).
********************************************************************************
********** Exercise 2 (results of a shot put competition) *************
* Reading in the data.
use shotput.dta, clear
* 1. Make a scatter plot of the data and calculate relevant summary statistics.
* Make interaction plots.
egen sexyear=group(sex year)
scatter distance sexyear, xlabel(1 "M, 1998" 2 "M, 1999" 3 "M, 2000" ///
4 "F, 1998" 5 "F, 1999" 6 "F, 2000")
table sex, by(year) c(mean distance sd distance)
anova distance sex##year
margins sex#year
marginsplot, xdimension(sex) title("") noci saving(g1, replace) xtitle("") ///
ytitle("Mean shot put distance (meters)") ///
legend(position(1) ring(0) row(3) region(lp(blank)))
marginsplot, xdimension(year) title("") noci saving(g2, replace) xtitle("") ///
ytitle("") legend(position(3) ring(0) row(2) region(lp(blank)))
graph combine g1.gph g2.gph
* On average the men consistently throw the shot put to a greater distance than
* the women. Is there a slight increase in the distance over the three year
* period? The variation (standard deviation) is more or less the same for all
* combinations of sex and year. Is there interaction between sex and year? Not
* much even though it looks as if the average distance decreases for the men
* between 1999 and 2000 whereas that is not the case for the women, but there
* is substantial variation in the data.
* 2. Compare the six groups using a one-way analysis of variance model.
* Include a test for equal variances in the six groups in the analysis.
oneway distance sexyear
* Analysis of Variance
* Source SS df MS F Prob > F
* ------------------------------------------------------------------------
* Between groups 66.4841615 5 13.2968323 7.20 0.0000
* Within groups 88.6447061 48 1.84676471
* ------------------------------------------------------------------------
* Total 155.128868 53 2.92695976
*
* Bartlett's test for equal variances: chi2(5) = 0.9557 Prob>chi2 = 0.966
* The Bartlett-test tells us that there is no evidence against the hypothesis
* of equal variation (standard deviation) in the six groups (combinations of
* sex and year), p=0.97. In the analysis of variance table we see that we
* clearly reject the hypythesis of equal mean shot put distance in the six
* groups, p<0.0001.
* 3. Compare the six groups using a two-way analysis of variance model.
* Make a QQ-plot of the residuals.
anova distance sex##year
predict res, residual
qnorm res
* Number of obs = 54 R-squared = 0.4286
* Root MSE = 1.35896 Adj R-squared = 0.3691
*
* Source | Partial SS df MS F Prob>F
* -----------+----------------------------------------------------
* Model | 66.484161 5 13.296832 7.20 0.0000
* |
* sex | 62.060221 1 62.060221 33.60 0.0000
* year | 3.0638374 2 1.5319187 0.83 0.4424
* sex#year | 1.3601033 2 .68005165 0.37 0.6939
* |
* Residual | 88.644706 48 1.8467647
* -----------+----------------------------------------------------
* Total | 155.12887 53 2.9269598
* Here we only look at the test corresponding to the interaction (sex#year).
* We conclude that there is not a significant interaction between sex and year,
* p=0.69. In order to assess whether there is a main effect of sex and year,
* we need to consider the additive model. From the QQ-plot we conclude that
* there is no reason to question the twoway ANOVA model for these data.
anova distance sex year
* Number of obs = 54 R-squared = 0.4198
* Root MSE = 1.34168 Adj R-squared = 0.3850
*
* Source | Partial SS df MS F Prob>F
* -----------+----------------------------------------------------
* Model | 65.124058 3 21.708019 12.06 0.0000
* |
* sex | 62.060221 1 62.060221 34.48 0.0000
* year | 3.0638374 2 1.5319187 0.85 0.4331
* |
* Residual | 90.004809 50 1.8000962
* -----------+----------------------------------------------------
* Total | 155.12887 53 2.9269598
* Here we see that there is not a significant main effect of year, p=0.43,
* but that there is a clearly significant main effect of sex, p<0.0001.
margins sex#year
* | Margin Std. Err. t P>|t| [95% Conf. Interval]
* -------------+--------------------------------------------------------------
* sex#year |
* Man#1998 | 10.53815 .3651581 28.86 0.000 9.804706 11.27159
* Man#1999 | 10.93481 .3651581 29.95 0.000 10.20137 11.66826
* Man#2000 | 11.10704 .3651581 30.42 0.000 10.3736 11.84048
* Woman#1998 | 8.394074 .3651581 22.99 0.000 7.660632 9.127516
* Woman#1999 | 8.790741 .3651581 24.07 0.000 8.057299 9.524182
* Woman#2000 | 8.962963 .3651581 24.55 0.000 8.229521 9.696405
pwcompare sex year, eff
* | Contrast Std. Err. t P>|t| [95% Conf. Interval]
* --------------+--------------------------------------------------------------
* sex |
* Woman vs Man | -2.144074 .3651581 -5.87 0.000 -2.877516 -1.410632
* |
* year |
* 1999 vs 1998 | .3966667 .4472255 0.89 0.379 -.5016123 1.294946
* 2000 vs 1998 | .5688889 .4472255 1.27 0.209 -.32939 1.467168
* 2000 vs 1999 | .1722222 .4472255 0.39 0.702 -.7260567 1.070501
* 4. Write a short summary of statistical methods used in the analysis and the
* findings.
* The data corresponding to distance of shot put for men and women over a three
* year period were analyzed using a twoway ANOVA model. Model validation was
* performed by testing for equal standard deviations in the six groups (a
* Bartlett test gave a p-value of 0.97) and by inspecting a QQ-plot of the
* residuals which did not raise any concerns about the validity of the analysis.
* The interaction between sex and year was not statistically significant,
* p=0.43, and a subsequent analysis including only the main effects of sex and
* year resulted in the conclusion that only sex had a significant effect,
* p<0.0001. The effect of year was not statistically significant, p=0.43. The
* estimated mean distances were:
*
* Men Women
* 1998 10.54 [ 9.80;11.27] 8.39 [7.66;9.13]
* 1999 10.93 [10.20;11.67] 8.79 [8.06;9.52]
* 2000 11.11 [10.37;11.84] 8.96 [8.23;9.70]
*
* On average men threw the shot put 2.14 meters further than women, and with
* 95% certainty between 1.41 and 2.88 meters further. There was a slight
* increase in the distance over the three year period (57 cm, 95%-CI:
* -33 - 147 cm from 1998 to 1999, p=0.38, and 17 cm, 95%-CI: -73 - 107 cm from
* 1999 to 2000, p=0.70), though that was not significant.
*******************************************************************************
********** Exercise 3 (weight and height in three groups) *************
* Reading in the data.
use wh.dta, clear
* 1. Compare group 1 and group 2 with respect to the weight and height
* measurements using two t-tests.
ttest weight if group==1 | group==2, by(group)
ttest height if group==1 | group==2, by(group)
* Group | Obs Mean Std. Err. Std. Dev. [95% Conf. Interval]
* ---------+-------------------------------------------------------------------
* 1 | 25 76.024 2.586956 12.93478 70.68478 81.36321
* 2 | 25 83.14 2.495249 12.47624 77.99006 88.28994
* ---------+-------------------------------------------------------------------
* combined | 50 79.582 1.849891 13.0807 75.86451 83.29949
* ---------+-------------------------------------------------------------------
* diff | -7.116001 3.594247 -14.34272 .1107165
* Pr(T < t) = 0.0267 Pr(|T| > |t|) = 0.0535 Pr(T > t) = 0.9733
* Group | Obs Mean Std. Err. Std. Dev. [95% Conf. Interval]
* ---------+-------------------------------------------------------------------
* 1 | 25 171.32 2.778321 13.8916 165.5858 177.0542
* 2 | 25 179.2 2.910899 14.5545 173.1922 185.2078
* ---------+-------------------------------------------------------------------
* combined | 50 175.26 2.06937 14.63266 171.1014 179.4186
* ---------+-------------------------------------------------------------------
* diff | -7.88 4.023978 -15.97075 .2107503
* Pr(T < t) = 0.0280 Pr(|T| > |t|) = 0.0560 Pr(T > t) = 0.9720
* We see that on average the people in group 2 are both heavier (7.1 kg,
* 95%-CI: -0.1 - 14.3 kg) and taller (7.9 cm, 95%-CI: -0.2 - 16.0 cm) than
* those in group 1. The two groups are not statistically different regarding
* mean weight (p=0.054) and height (p=0.056), though it is close.
* 2. Compare group 1 and group 2 with respect to the weight and height
* measurements using MANOVA.
mvtest mean weight height if group==1 | group==2, by(group)
* Test for equality of 2 group means, assuming homogeneity
* | Statistic F(df1, df2) = F Prob>F
* -----------------------+-----------------------------------------------
* Wilks' lambda | 0.9242 2.0 47.0 1.93 0.1569 e
* Pillai's trace | 0.0758 2.0 47.0 1.93 0.1569 e
* Lawley-Hotelling trace | 0.0820 2.0 47.0 1.93 0.1569 e
* Roy's largest root | 0.0820 2.0 47.0 1.93 0.1569 e
* -----------------------------------------------------------------------
* e = exact, a = approximate, u = upper bound on F
* The multivariate test gives that we accept the hypothesis of no difference
* between groups 1 and 2 based on the mean weight and height (p=0.16).
* 3. Make a scatter plot similar to that on slide 1-34. What are your comments?
egen weightmean = mean(weight), by(group)
egen heightmean = mean(height), by(group)
twoway ///
(scatter height weight if group==1, ms(O) mc(blue)) ///
(scatter heightmean weightmean if group==1, m(S) msiz(large) mc(blue)) ///
(scatter height weight if group==2, ms(O) mc(red)) ///
(scatter heightmean weightmean if group==2, m(S) msiz(large) mc(red)), ///
legend(label(1 "Group 1") label(2 "Mean") label(3 "Group 2") ///
label(4 "Mean") position(11) ring(0) row(2) region(lp(blank))) ///
xtitle("Weight (kg)") ytitle("Height (cm)") xscale(range(50 115)) ///
yscale(range(140 215))
* Comments: In both groups weight and height follow each other so that heavy
* people tend to be taller. The means are quite far apart (difference in mean
* weight of around 7 kg and in mean height of around 8 cm). We see, however,
* that the difference in means between the two groups is along the axis
* dictated by the correlation between the height and weight measurements.
* 4. Repeat questions 1 to 3 for group 1 and 3.
ttest weight if group==1 | group==3, by(group)
ttest height if group==1 | group==3, by(group)
* Two-sample t test with equal variances
* -----------------------------------------------------------------------------
* Group | Obs Mean Std. Err. Std. Dev. [95% Conf. Interval]
* ---------+-------------------------------------------------------------------
* 1 | 25 76.024 2.586956 12.93478 70.68478 81.36321
* 3 | 25 78.5 1.970525 9.852623 74.43304 82.56696
* ---------+-------------------------------------------------------------------
* combined | 50 77.262 1.618997 11.44803 74.0085 80.5155
* ---------+-------------------------------------------------------------------
* diff | -2.475999 3.25197 -9.014523 4.062525
* Pr(T < t) = 0.2251 Pr(|T| > |t|) = 0.4502 Pr(T > t) = 0.7749
* Two-sample t test with equal variances
* -----------------------------------------------------------------------------
* Group | Obs Mean Std. Err. Std. Dev. [95% Conf. Interval]
* ---------+-------------------------------------------------------------------
* 1 | 25 171.32 2.778321 13.8916 165.5858 177.0542
* 3 | 25 170.16 2.266186 11.33093 165.4828 174.8372
* ---------+-------------------------------------------------------------------
* combined | 50 170.74 1.776217 12.55975 167.1706 174.3094
* ---------+-------------------------------------------------------------------
* diff | 1.16 3.585341 -6.04881 8.36881
* Pr(T < t) = 0.6262 Pr(|T| > |t|) = 0.7477 Pr(T > t) = 0.3738
* We see that the estimated mean weight in group 3 is 2.5 kg (95%-CI: -4.1 -
* 9.0 kg) higher than in group 1. The estimated mean height, however, is 1.2 cm
* higher in group 1 (95%-CI: -6.0 - 8.4 cm) compared to group 3. Neither of
* these differences are significant on a 5% level (weight: p=0.45, height:
* p=0.75).
mvtest mean weight height if group==1 | group==3, by(group)
* Test for equality of 2 group means, assuming homogeneity
* | Statistic F(df1, df2) = F Prob>F
* -----------------------+-----------------------------------------------
* Wilks' lambda | 0.7828 2.0 47.0 6.52 0.0032 e
* Pillai's trace | 0.2172 2.0 47.0 6.52 0.0032 e
* Lawley-Hotelling trace | 0.2774 2.0 47.0 6.52 0.0032 e
* Roy's largest root | 0.2774 2.0 47.0 6.52 0.0032 e
* -----------------------------------------------------------------------
* e = exact, a = approximate, u = upper bound on F
* Comparing the two groups using MANOVA results in a clearly significant
* difference between the two groups when based on the mean weight and height
* (p=0.0032).
twoway ///
(scatter height weight if group==1, ms(O) mc(blue)) ///
(scatter heightmean weightmean if group==1, m(S) msiz(large) mc(blue)) ///
(scatter height weight if group==3, ms(O) mc(red)) ///
(scatter heightmean weightmean if group==3, m(S) msiz(large) mc(red)), ///
legend(label(1 "Group 1") label(2 "Mean") label(3 "Group 3") ///
label(4 "Mean") position(11) ring(0) row(2) region(lp(blank))) ///
xtitle("Weight (kg)") ytitle("Height (cm)") xscale(range(50 95)) ///
yscale(range(140 200))
* Also in group 3 we see a strong positive correlation between weight and
* height. The two group means are not that far apart, though, around 2 kg
* in weight and 1 cm in height. However, the means differ away from the axis
* dictated by the correlation between weight and height, and that is why when
* you take into account this correlation the difference between the two groups
* is clearly significant.
* 5. What are the assumptions behind the analyses? Do they seem reasonable in
* this situation?
* From slide 1-36 we note the following assumptions:
*
* 1. The observations are normally distributed
* 2. Independent observations for different subjects
* 3. Standard deviations and correlations are the same in the three groups
* Regarding 1. normality:
hist weight, normal by(group)
hist height, normal by(group)
* There are no strong deviations from normality in either group.
* Regarding 2. independence: As usual this has to be assessed from a
* description of how the data were retrieved (the design). In this example
* this is not an issue as different unrelated subjects enter the three groups.
* Regarding 3. Same standard deviations and correlations in the three groups:
bysort group: pwcorr weight height
mvtest covariances weight height, by(group)
* Test of equality of covariance matrices across 3 samples
* Modified LR chi2 = 6.541967
* Box F(6,129201.2) = 1.05 Prob > F = 0.3927
* Box chi2(6) = 6.28 Prob > chi2 = 0.3926
* Here we see that the correlations are very high in all three groups and
* almost the same (0.93-0.98). A test for the hypothesis of equal standard
* deviations and correlations in the three groups is accepted with a p-value
* of 0.39.
* 6. Write a short summary of statistical methods used in the analysis and the
* findings.
* Using simple t-tests we saw that it was borderline significant that group 1
* and 2 differed both with regard to weight (2-1: 7.1 kg, 95%-CI: -0.1 - 14.3
* kg, p=0.054) and height (2-1: 7.9 cm, 95%-CI: -0.2 - 16.0 cm, p=0.056). Based
* the MANOVA we concluded that taking into account the correlation between
* height and weight, group 1 and 2 were not significantly different, p=0.16.
* With regard to group 1 and 3, then they did not differ significantly with
* respect to weight (3-1: 2.5 kg, 95%-CI: -4.1 - 9.0 kg, p=0.45) and height
* (1-3: 1.2 cm, 95%-CI: -6.0 - 8.4 cm, p=0.75). However, combining weight and
* height in a MANOVA we concluded that there was a clearly significant
* difference between group 1 and 3, p=0.0032. A plot of the means with contour
* ellipses may be a good way to illustrate the findings and help to understand
* them.
ellip height weight, by(group) lco(blue red black) constant(t2) ///
plot((scatter heightmean weightmean if group==1, m(S) mc(blue)) ///
(scatter heightmean weightmean if group==2, m(S) mc(red)) ///
(scatter heightmean weightmean if group==3, m(S) mc(black))) ///
overlay legend(label(4 "Group 1") label(5 "Group 2") label(6 "Group 3") ///
order(- - - 4 5 6) colfirst position(11) ring(0) row(3) ///
region(lp(blank))) title("") subtitle("") note("")