*******************************************************************************
*** ***
*** Stata-program (do-file) associated with the Exercises on day two. ***
*** ***
*** ***
*******************************************************************************
* Specifying the path to where the data are, for example
cd "C:\ANOVA\DAY 2"
********** Exercise 4 (repeated leukocyte counts) *************
* Reading in the data.
use leukocyte.dta, clear
* 1. Plot the individual leukocyte curves, the grouped individual curves and
* the mean group curves.
reshape long plat leuk m, i(id) j(time)
scatter leuk time, connect(L) sort(time) by(id)
scatter leuk time, connect(L) sort(id time) by(grp)
* No very abnormal observations and more or less the same variation within and
* between individuals.
egen groupmean = mean(leuk), by(grp time)
sort grp time id
twoway (connected groupmean time if grp==1, msymbol(Oh)) ///
(connected groupmean time if grp==2, msymbol(O)), ///
legend(label(1 "CPB") label(2 "Sham")) ytitle("Mean leukocyte count")
* We see a decrease in the mean leukocyte curves in both groups followed by a
* steady increase. The increase happens at an earlier time-point in the sham
* group than in the CPB group.
* 2. Calculate the standard deviations and correlations corresponding to each
* group.
drop groupmean
reshape wide plat leuk m, i(id) j(time)
bysort grp: pwcorr leuk1 leuk2 leuk3 leuk4 leuk5 leuk6
tabstat leuk1 leuk2 leuk3 leuk4 leuk5 leuk6, nototal by(grp) stat(sd)
* -> grp = CPB
*
* | leuk1 leuk2 leuk3 leuk4 leuk5 leuk6
* -------------+------------------------------------------------------
* leuk1 | 1.0000
* leuk2 | 0.7093 1.0000
* leuk3 | 0.8598 0.6137 1.0000
* leuk4 | 0.2168 0.4490 0.4349 1.0000
* leuk5 | 0.3995 0.4965 0.5339 0.8705 1.0000
* leuk6 | 0.5274 0.7055 0.4901 0.6949 0.8671 1.0000
*
* -> grp = Sham
*
* | leuk1 leuk2 leuk3 leuk4 leuk5 leuk6
* -------------+------------------------------------------------------
* leuk1 | 1.0000
* leuk2 | 0.8106 1.0000
* leuk3 | 0.8511 0.8676 1.0000
* leuk4 | 0.8667 0.7891 0.9246 1.0000
* leuk5 | 0.6953 0.8050 0.9135 0.8225 1.0000
* leuk6 | 0.6183 0.6785 0.8193 0.7883 0.9270 1.0000
*
* Summary statistics: sd
* by categories of: grp (Treatment group)
*
* grp | leuk1 leuk2 leuk3 leuk4 leuk5 leuk6
* -----+------------------------------------------------------------
* CPB | 3.946463 3.568406 3.853421 5.336915 4.739228 4.678441
* Sham | 4.042359 3.710153 5.295524 5.354864 6.006996 6.27745
* We observe slightly higher standard deviations and correlations in the sham
* group compared to the CPB group.
* 3. Analyze the leukocyte count data using the multivariate repeated
* measurements ANOVA.
generate leuk21 = leuk2 - leuk1
generate leuk32 = leuk3 - leuk2
generate leuk43 = leuk4 - leuk3
generate leuk54 = leuk5 - leuk4
generate leuk65 = leuk6 - leuk5
mvtest means leuk21 leuk32 leuk43 leuk54 leuk65, by(grp)
* Test for equality of 2 group means, assuming homogeneity
*
* | Statistic F(df1, df2) = F Prob>F
* -----------------------+-----------------------------------------------
* Wilks' lambda | 0.3764 5.0 24.0 7.95 0.0002 e
* Pillai's trace | 0.6236 5.0 24.0 7.95 0.0002 e
* Lawley-Hotelling trace | 1.6569 5.0 24.0 7.95 0.0002 e
* Roy's largest root | 1.6569 5.0 24.0 7.95 0.0002 e
* -----------------------------------------------------------------------
* e = exact, a = approximate, u = upper bound on F
* We conclude that there is clear evidence against the hypothesis of equal
* evolvement over time in the two groups (Test 1: p=0.0002). Of course this
* comes as no surprise having seen the plot of the mean curves.
* Let us test the hypothesis of equal standard deviations and correlations (for
* the differences) in the two groups:
mvtest covariances leuk21 leuk32 leuk43 leuk54 leuk65, by(grp)
* Test of equality of covariance matrices across 2 samples
*
* Modified LR chi2 = 39.11329
* Box F(15,3156.6) = 2.10 Prob > F = 0.0078
* Box chi2(15) = 31.66 Prob > chi2 = 0.0072
* There is some evidence against the hypothesis of equal standard deviations
* and correlations for the differences in the two groups (p=0.0078). This is
* probably of little importance as the groups are of equal size (n=15 in each
* group), but let of repeat Test 1 allowing for heterogenous variation:
mvtest means leuk21 leuk32 leuk43 leuk54 leuk65, het by(grp)
* Test for equality of 2 group means, allowing for heterogeneity
*
* MNV F(5,22.4) = 7.87
* Prob > F = 0.0002
* As suspected this leads to almost the same p-value (p=0.0002) and exactly the
* same conclusion.
* Let us finish the analysis (before inspecting the residuals) by giving the
* group differences at the different time-points.
drop leuk21 leuk32 leuk43 leuk54 leuk65
reshape long plat leuk m, i(id) j(time)
xi: mixed leuk bn.time#bn.grp || id: i.time, cov(un)
contrast grp@time, eff nowald
* | Contrast Std. Err. z P>|z| [95% Conf. Interval]
* ------------------+----------------------------------------------------------
* leuk |
* grp@time |
* (Sham vs base) 1 | .4066668 1.409257 0.29 0.773 -2.355426 3.168759
* (Sham vs base) 2 |-.1399999 1.284092 -0.11 0.913 -2.656775 2.376775
* (Sham vs base) 3 | 5.32 1.633687 3.26 0.001 2.118033 8.521967
* (Sham vs base) 4 | 4.893333 1.885894 2.59 0.009 1.19705 8.589617
* (Sham vs base) 5 | 4.22 1.908627 2.21 0.027 .4791595 7.960841
* (Sham vs base) 6 | 4.82 1.952946 2.47 0.014 .9922971 8.647703
* We see more or less the same leukocyte level in the two groups at time-points
* 1 and 2, whereas the sham group has a significantly higher mean level at the
* later time-points.
* 4. Make a QQ-plot for the residuals.
predict res, rstandard
qnorm res
* The is no reason to distrust the conclusions of the analysis based on the
* multivariate repeated measurement ANOVA based on this plot, as there are no
* clear systematic deviations from a straight line in the QQ-plot.
* 5. Plot the individual m-curves, the grouped individual curves and the mean
* group curves. What result from the multivariate repeated measurements
* ANOVA based on m1,..., m6 would you expect?
scatter m time, connect(L) sort(time) by(id)
scatter m time, connect(L) sort(id time) by(grp)
* Very much the same conclusions as for the leukocyte counts: No very abnormal
* observations and more or less the same variation within and between
* individuals.
egen groupmean = mean(m), by(grp time)
sort grp time id
twoway (connected groupmean time if grp==1, msymbol(Oh)) ///
(connected groupmean time if grp==2, msymbol(O)), ///
legend(label(1 "CPB") label(2 "Sham")) ytitle("Mean m-values")
* We don't see the same pattern in the mean curves. It still looks as if the
* difference between the two groups depends on which time-point we consider,
* so we would expect to reject Test 1 also for the variable m.
* 6. Perform the analysis and compare with the results of the analysis in
* question 3 (more specifically compare the p-values).
drop groupmean _Itime_2 _Itime_3 _Itime_4 _Itime_5 _Itime_6 res
reshape wide plat leuk m, i(id) j(time)
generate m21 = m2 - m1
generate m32 = m3 - m2
generate m43 = m4 - m3
generate m54 = m5 - m4
generate m65 = m6 - m5
mvtest means m21 m32 m43 m54 m65, by(grp)
* From the output we see that we get exactly the same F-test statistic and
* thereby the same p-value as for the leukocyte data. As mentioned on slide
* 2-24 the multivariate repeated measurement ANOVA is not influenced by a
* reordering of the time-points (nor by a linear transformation of the
* observations). The tests and model validation are therefor the same for m as
* for the leukocyte counts. Of course a comparison of the two groups at the
* different time-points will change, but if you kept track of which time-points
* that were interchanged you would rediscover the group differences up to the
* scaling factor used to create m.
drop m21 m32 m43 m54 m65
reshape long plat leuk m, i(id) j(time)
xi: mixed m bn.time#bn.grp || id: i.time, cov(un)
contrast grp@time, eff nowald
* | Contrast Std. Err. z P>|z| [95% Conf. Interval]
* ------------------+----------------------------------------------------------
* m |
* grp@time |
* (Sham vs base) 1 | 26.6 8.16838 3.26 0.001 10.59027 42.60973
* (Sham vs base) 2 | 21.1 9.543242 2.21 0.027 2.395589 39.80441
* (Sham vs base) 3 | 2.033333 7.046059 0.29 0.773 -11.77669 15.84336
* (Sham vs base) 4 | 24.1 9.764783 2.47 0.014 4.961377 43.23862
* (Sham vs base) 5 | -.7 6.420568 -0.11 0.913 -13.28408 11.88408
* (Sham vs base) 6 | 24.46667 9.429483 2.59 0.009 5.98522 42.94811
* Notice that the p-values are exactly the same (the ordering is just
* different, 1=3, 2=5, 3=1, 4=6, 5=2 and 6=4). The differences and 95%-CI are
* also the same except for a factor of 5.
*
*******************************************************************************
***** Exercise 5 (crp in pigs undergoing CPB or a sham-operation) ******
* Reading in the data.
use crp.dta, clear
* 1. Plot the individual curves, the grouped individual curves and the mean
* group curves.
scatter crp time, connect(L) sort(time) by(id)
scatter crp time, connect(L) sort(id time) by(group)
* No very abnormal observations but somewhat larger between pig variation in
* the sham group perhaps confined to the time-points where the crp-levels are
* higher (3 and 4).
egen groupmean = mean(crp), by(group time)
sort group time id
twoway (connected groupmean time if group==1, msymbol(Oh)) ///
(connected groupmean time if group==2, msymbol(O)), ///
legend(label(1 "CPB") label(2 "Sham")) ytitle("Mean crp-values")
* Increasing mean crp-levels over time in both groups. Slightly larger increase
* in the sham-operated group for some reason. The average crp-level in the sham
* group was also higher at baseline!
* 2. Calculate the standard deviations and correlations corresponding to each
* group.
drop groupmean
reshape wide crp, i(id) j(time)
bysort group: pwcorr crp1 crp2 crp3 crp4
tabstat crp1 crp2 crp3 crp4, nototal by(group) stat(sd)
* -> group = 1
*
* | crp1 crp2 crp3 crp4
* -------------+------------------------------------
* crp1 | 1.0000
* crp2 | 0.7436 1.0000
* crp3 | 0.6539 0.3294 1.0000
* crp4 | 0.4691 0.3804 0.6606 1.0000
*
* -> group = 2
*
* | crp1 crp2 crp3 crp4
* -------------+------------------------------------
* crp1 | 1.0000
* crp2 | 0.7396 1.0000
* crp3 | 0.6019 0.7638 1.0000
* crp4 | 0.5056 0.5786 0.6916 1.0000
*
* Summary statistics: sd
* by categories of: group (Group)
*
* group | crp1 crp2 crp3 crp4
* ---------+----------------------------------------
* 1 | 1.304944 1.481093 1.39975 2.087349
* 2 | 2.020286 1.888297 3.980121 4.006938
* There is much larger standard deviations in the sham group compared to the
* CPB-group especially at time-points 3 and 4 in accordance with what we
* observed in the plot of the grouped individual curves. The correlations are
* more similar.
* 3. Analyze the crp data using the multivariate repeated measurements ANOVA.
generate crp21 = crp2 - crp1
generate crp32 = crp3 - crp2
generate crp43 = crp4 - crp3
mvtest means crp21 crp32 crp43, by(group)
* Test for equality of 2 group means, assuming homogeneity
*
* | Statistic F(df1, df2) = F Prob>F
* -----------------------+-----------------------------------------------
* Wilks' lambda | 0.5834 3.0 18.0 4.28 0.0190 e
* Pillai's trace | 0.4166 3.0 18.0 4.28 0.0190 e
* Lawley-Hotelling trace | 0.7141 3.0 18.0 4.28 0.0190 e
* Roy's largest root | 0.7141 3.0 18.0 4.28 0.0190 e
* -----------------------------------------------------------------------
* e = exact, a = approximate, u = upper bound on F
* We conclude that there is some evidence against the hypothesis of equal group
* differences at all four time-points (Test 1: p=0.0190). This is also what we
* would from looking at the mean curves as the increase in crp over time is
* higher in the sham-group compared to the CPB-group.
* Let us test the hypothesis of equal standard deviations and correlations (for
* the differences) in the two groups:
mvtest covariances crp21 crp32 crp43, by(group)
* Test of equality of covariance matrices across 2 samples
*
* Modified LR chi2 = 16.2534
* Box F(6,2613.3) = 2.26 Prob > F = 0.0356
* Box chi2(6) = 13.58 Prob > chi2 = 0.0347
* There is some evidence against the hypothesis of equal standard deviations
* and correlations in the two groups, p=0.0356. As the two groups are almost
* of the same size (12 in CPB-group and 10 in the sham-group), this is probably
* not that critical.
mvtest means crp21 crp32 crp43, het by(group)
* Test for equality of 2 group means, allowing for heterogeneity
*
* MNV F(3,11.8) = 3.67
* Prob > F = 0.0444
* Making Test 1 without assuming equal standard deviations and correlations in
* the two groups only changes the p-value marginally and not the conclusion:
* There is some evidence against the hypothesis of equal development over time
* in the two groups.
* Model validation:
drop crp21 crp32 crp43
reshape long crp, i(id) j(time)
xi: mixed crp bn.time#bn.group || id: i.time, cov(un) technique(bfgs)
predict res, rstandard
qnorm res
* No clear deviations from a straight line in the QQ-plot so no reason here to
* reject the findings based on the multivariate repeated measurement ANOVA.
contrast group@time, eff nowald
* | Contrast Std. Err. z P>|z| [95% Conf. Interval]
* ---------------+-------------------------------------------------------------
* crp |
* group@time |
* (2 vs base) 1 | 1.476667 .6797322 2.17 0.030 .1444161 2.808917
* (2 vs base) 2 | 2.32 .6844832 3.39 0.001 .9784375 3.661562
* (2 vs base) 3 | 3.450833 1.169364 2.95 0.003 1.158922 5.742745
* (2 vs base) 4 | 5.754167 1.266337 4.54 0.000 3.272191 8.236142
* Increasing differences between the two groups with increasing time. Even at
* baseline the sham-group has a significantly higher mean crp-level than the
* CPB-group (1.48, 95%-CI: 0.14 - 2.81, p=0.03).
* 4. Repeat questions 1-3 on the log-transformed crp data and compare with the
* results of the analysis of the data on the original scale.
drop _Itime_2 _Itime_3 _Itime_4 res
generate logcrp = log(crp)
scatter logcrp time, connect(L) sort(time) by(id)
scatter logcrp time, connect(L) sort(id time) by(group)
* No very abnormal observations but and more similar within and between
* pig variation in the two groups compared to the analysis of the original
* crp-data.
egen groupmean = mean(logcrp), by(group time)
sort group time id
twoway (connected groupmean time if group==1, msymbol(Oh)) ///
(connected groupmean time if group==2, msymbol(O)), ///
legend(label(1 "CPB") label(2 "Sham")) ytitle("Mean log(crp-values)")
* In the plot of the mean crp-curves we see very parallel curves on the
* log-scale. Of course it is still a little strange that the sham-group starts
* out on a higher level at baseline, and any difference between the two groups
* cannot be attributed to the treatment.
drop groupmean
reshape wide crp logcrp, i(id) j(time)
bysort group: pwcorr logcrp1 logcrp2 logcrp3 logcrp4
tabstat logcrp1 logcrp2 logcrp3 logcrp4, nototal by(group) stat(sd)
* The conclusions here have not changed much: We still see a larger variation
* and correlation in the sham-group compared to the CPB-group.
generate logcrp21 = logcrp2 - logcrp1
generate logcrp32 = logcrp3 - logcrp2
generate logcrp43 = logcrp4 - logcrp3
mvtest means logcrp21 logcrp32 logcrp43, by(group)
* Test for equality of 2 group means, assuming homogeneity
*
* | Statistic F(df1, df2) = F Prob>F
* -----------------------+-----------------------------------------------
* Wilks' lambda | 0.9305 3.0 18.0 0.45 0.7218 e
* Pillai's trace | 0.0695 3.0 18.0 0.45 0.7218 e
* Lawley-Hotelling trace | 0.0746 3.0 18.0 0.45 0.7218 e
* Roy's largest root | 0.0746 3.0 18.0 0.45 0.7218 e
* -----------------------------------------------------------------------
* e = exact, a = approximate, u = upper bound on F
* For the log-transformed crp-data we clearly accept the hypothesis of equal
* development over time in the two groups (p=0.72). Having accepted that the
* mean curves are parallel, we can proceed and test if the mean log crp-curves
* are the same in the two groups (Test 2):
generate ave = (logcrp1+logcrp2+logcrp3+logcrp4)/4
sdtest ave, by(group)
ttest ave, by(group)
* This t-test gives a p-value of 0.0029, so we conclude that there is clear
* evidence against the hypothesis of coinciding mean log-crp curves in the two
* groups.
* What about the variations and correlations in the two groups?
mvtest covariances logcrp21 logcrp32 logcrp43, by(group)
* Test of equality of covariance matrices across 2 samples
*
* Modified LR chi2 = 8.567527
* Box F(6,2613.3) = 1.19 Prob > F = 0.3088
* Box chi2(6) = 7.16 Prob > chi2 = 0.3066
* As we suspected from individual and grouped individual curves there is no
* real evidence against the hypothesis of equal standard deviations and
* correlations (p=0.31) even though this is probably not so important in this
* situation.
drop logcrp21 logcrp32 logcrp43
reshape long crp logcrp, i(id) j(time)
xi: mixed logcrp bn.time bn.group || id: i.time, cov(un) technique(bfgs)
predict res, rstandard
qnorm res
* Again we see no clear deviations from a straight line in the QQ-plot so no
* reason here to reject the findings based on the multivariate repeated
* measurement ANOVA for the log-transformed crp-values.
contrast group, eff nowald eform
* | exp(b) Std. Err. z P>|z| [95% Conf. Interval]
* -------------+--------------------------------------------------------------
* logcrp |
* group |
* (2 vs base) | 1.859195 .2661316 4.33 0.000 1.404368 2.461325
* Conclusion (regarding the difference between the two groups): Based on an
* analysis of the log-transformed crp-data we conclude that the median level of
* crp is 86% (95%-CI: 40-146%, p<0.0001) higher in the sham-group compared to
* the CPB-group. This holds at all time-points. There are only marginal
* differences between the analyses based on the original data and the
* log-transformed data (based on the model-validation performed), but it is
* much simpler to present the results regarding the treatment effect based on
* the analysis on the log-scale, and that is a valid argument for choosing that
* analysis.
* 5. What are the assumptions behind the analysis in question 3? Check the
* validity of the assumptions.
* The assumptions are given on slide 2-23:
*
* 1. Independent observations for different pigs
* 2. The observations are multivariate normally distributed
* 3. Standard deviations and correlations are the same in the two groups
*
* We have checket the validity of the assumptions by inspecting QQ-plots of the
* residuals and by testing that the standard deviations and correlations are
* the same in the two groups. Both for the original data and the
* log-transformed data we saw QQ-plots without any clear deviations from a
* straight line. Regarding the test for equal standard deviations and
* correlations, we got a p-value of 0.0356 for the original data and 0.31 for
* the log-transformed data. At a first glance you would think that it is much
* better to proceed with the analysis of the log-scale, but with almost equal
* group sizes, we would not expect any problems with trusting the results based
* on an analysis of the original data.
*
*******************************************************************************
***** Exercise 6 (crp in pigs undergoing CPB or a sham-operation) ******
* Reading in the data.
use crp.dta, clear
* 1. Analyze the log-transformed data using a univariate repeated measurements
* ANOVA.
generate logcrp=log(crp)
anova logcrp group/id time time#group, bse(id) repeated(time)
* Number of obs = 88 R-squared = 0.8056
* Root MSE = .340988 Adj R-squared = 0.7181
*
* Source | Partial SS df MS F Prob>F
* -----------+----------------------------------------------------
* Model | 28.901513 27 1.0704264 9.21 0.0000
* |
* group | 1.8394974 1 1.8394974 2.93 0.1024
* id | 12.552354 20 .62761772
* -----------+----------------------------------------------------
* time | 9.0776701 3 3.02589 26.02 0.0000
* time#group | .19688601 3 .06562867 0.56 0.6406
* |
* Residual | 6.976377 60 .11627295
* -----------+----------------------------------------------------
* Total | 35.87789 87 .41238954
*
* Between-subjects error term: id
* Levels: 22 (20 df)
* Lowest b.s.e. variable: id
* Repeated variable: time
* Huynh-Feldt epsilon = 1.0657
* *Huynh-Feldt epsilon reset to 1.0000
* Greenhouse-Geisser epsilon = 0.8723
* Box's conservative epsilon = 0.3333
*
* ------------ Prob > F ------------
* Source | df F Regular H-F G-G Box
* -----------+----------------------------------------------------
* time | 3 26.02 0.0000 0.0000 0.0000 0.0001
* time#group | 3 0.56 0.6406 0.6406 0.6177 0.4612
* Residual | 60
* ----------------------------------------------------------------
* Also when we analyze the log-crp measurements using the univariate repeated
* measurements ANOVA we conclude that there is no evidence against the
* hypothesis of equal development over time (p=0.64). The different corrections
* to this p-value does not change this conclusion.
mixed logcrp b1.time b1.group || id: , reml
contrast group, eff nowald eform
* | exp(b) Std. Err. z P>|z| [95% Conf. Interval]
* -------------+---------------------------------------------------------------
* logcrp |
* group |
* (2 vs base) | 1.779089 .3017422 3.40 0.001 1.275939 2.480651
* Based on the univariate repeated measurements ANOVA we get more or less the
* same conclusion as in the situation with the multivariate repeated
* measurements ANOVA: The median crp-level is 78% (95%-CI: 28-148%, p=0.001)
* higher in the sham-group compared to the CPB-group. Again this holds for all
* time-points as we have accepted the hypothesis of equal differences between
* the mean log-crp levels in the two groups at all time-points (Test 1).
* 2. Estimate the within animal variation, the between animal variation, and
* the correlation between measurements on the same pig.
mixed logcrp b1.time b1.group || id: , nofetable stddev reml
estat icc
* Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
* -----------------------------+-----------------------------------------------
* id: Identity |
* sd(_cons) | .358384 .0695847 .2449499 .5243484
* -----------------------------+-----------------------------------------------
* sd(Residual) | .3374334 .030061 .2833721 .4018085
* Level | ICC Std. Err. [95% Conf. Interval]
* ----------------------------+------------------------------------------------
* id | .5300821 .1104516 .3211272 .7289978
* From this we get:
*
* Betwee pig variation: sigmaB = 0.3584
* Within pig variation: sigmaW = 0.3374
*
* The corelation between any two measurements on the same pig: rho=0.53
*
* The total variation is obtained from the quantities above:
*
* Total variation: sigmaT = 0.3584/sqrt(0.53) = 0.4922
* 3. What are the assumptions behind the analysis in question 1? Check the
* validity of the assumptions.
* The assumptions behind the univariate repeated measurements ANOVA are as for
* the multivariate repeated measurements ANOVA except that we also assume that
* the correlation between any two measurements on the same pig are identical
* and that the standard deviations are the same at all time-points.
predict fit if e(sample), xb
predict res if e(sample), rstandard
scatter res fit, yline(0) name(q1,replace)
qnorm res, title ("residual probability-plot") name(q2,replace)
graph combine q1 q2
* There are no clear patterns in the plot of the residuals against the fitted
* values, and no obvious deviations from a sraight line in the QQ-plot.
* 4. Compare with the findings in Exercise 5 question 4.
* We get very similar findings: Higher median crp-levels in the sham-group
* compared to the CPB-group:
*
* Multivariate repeated measurements ANOVA: 86%, 95%-CI: 40-146%, p<0.0001
* Univariate repeated measurements ANOVA: 78%, 95%-CI: 28-148%, p=0.001
*
* This holds for all time-points.
*
******************************************************************************
********** Exercise 7 (Cardiac output for CHF and healthy) *************
* Reading in the data.
use cardiac.dta, clear
egen subj = group(group subject)
* The variable subject is the person number within group. In order to get a
* unique identifier, we have created the variable subj.
* 1. Plot the individual curves, the grouped individual curves and the mean
* group curves.
scatter co workload, connect(L) sort(workload) by(subj)
scatter co workload, connect(L) sort(subj workload) by(group)
* No very abnormal observations but a larger between subject variation in
* the healthy group compared to the CHF group. The variation between subjects
* at workload 0 is particularly small in the CHF group. The variation within
* subjects is more or less the same i the two groups.
egen groupmean = mean(co), by(group workload)
sort group workload subj
twoway (connected groupmean workload if group==1, msymbol(Oh)) ///
(connected groupmean workload if group==2, msymbol(O)), ///
legend(label(1 "Healthy") label(2 "CHF")) ///
ytitle("Mean cardiac output (l/min)")
* We see increasing cardiac output with increasing workload. The healthy
* volunteer group has a higher mean cardiac output than the cronic heart
* failure patients at all workloads. The mean curves look very parallel, though
* the cardiac output at rest (workload 0) seems to be more similar than during
* exercise.
* 2. Calculate the standard deviations and correlations corresponding to each
* group.
drop groupmean subject
reshape wide co, i(subj) j(workload)
bysort group: pwcorr co0 co30 co60 co90 co120
tabstat co0 co30 co60 co90 co120, nototal by(group) stat(sd)
* -> group = Healthy
*
* | co0 co30 co60 co90 co120
* -------------+---------------------------------------------
* co0 | 1.0000
* co30 | 0.9447 1.0000
* co60 | 0.9387 0.9877 1.0000
* co90 | 0.9235 0.9591 0.9780 1.0000
* co120 | 0.9203 0.9447 0.9700 0.9716 1.0000
*
* -> group = CHF
*
* | co0 co30 co60 co90 co120
* -------------+---------------------------------------------
* co0 | 1.0000
* co30 | 0.8395 1.0000
* co60 | 0.7188 0.9315 1.0000
* co90 | 0.7791 0.9338 0.9586 1.0000
* co120 | 0.7206 0.9127 0.9525 0.9708 1.0000
*
* Summary statistics: sd
* by categories of: group (Healthy/CHF-patient)
*
* group | co0 co30 co60 co90 co120
* --------+--------------------------------------------------
* Healthy | 2.649228 3.209361 2.604309 2.522445 2.075598
* CHF | .7856189 1.534571 1.867667 1.950991 1.710101
* We see that the standard deviations are somewhat higher in the healthy group
* compared to the CHF group. This is the case at all workloads and especially
* at workload 0. The correlations are very high and similar in the two groups
* except the correlations with workload 0 in the CHF group, where the
* correlations are smaller. This is completely in accordance with what was
* observed in the plot of the grouped individual curves where the cardiac
* output measurements at workload 0 were VERY similar for CHF-patients.
* 3. Analyze the data using the univariate repeated measurements ANOVA.
reshape long co, i(subj) j(workload)
anova co group/subj workload workload#group, bse(subj) repeated(workload)
* Number of obs = 114 R-squared = 0.9729
* Root MSE = .685875 Adj R-squared = 0.9626
*
* Source | Partial SS df MS F Prob>F
* ---------------+----------------------------------------------------
* Model | 1383.7271 31 44.636357 94.89 0.0000
* |
* group | 93.025 1 93.025 4.44 0.0467
* subj | 460.66091 22 20.939132
* ---------------+----------------------------------------------------
* workload | 526.06546 4 131.51637 279.57 0.0000
* workload#group | 14.482621 4 3.6206552 7.70 0.0000
* |
* Residual | 38.574782 82 .47042417
* ---------------+----------------------------------------------------
* Total | 1422.3019 113 12.586742
*
* Between-subjects error term: subj
* Levels: 24 (22 df)
* Lowest b.s.e. variable: subj
*
* Repeated variable: workload
* Huynh-Feldt epsilon = 0.5586
* Greenhouse-Geisser epsilon = 0.4875
* Box's conservative epsilon = 0.2500
*
* ------------ Prob > F ------------
* Source | df F Regular H-F G-G Box
* ---------------+----------------------------------------------------
* workload | 4 279.57 0.0000 0.0000 0.0000 0.0000
* workload#group | 4 7.70 0.0000 0.0009 0.0016 0.0115
* Residual | 82
* --------------------------------------------------------------------
* We see clear evidence against the hypothesis of parallel curves in the two
* groups (Test 1: p<0.0001). This is also the conclusion if you consider the
* different adjustments to the F-test statstic that Stata provides.
mixed co bn.workload#bn.group, noconst || subj: , reml
predict fit if e(sample), fit
predict res if e(sample), rstandard
scatter res fit, yline(0) name(q1,replace)
qnorm res, title ("residual probability-plot") name(q2,replace)
graph combine q1 q2
* Looking at the plot of the residuals against the fitted values maybe you
* would argue that there is some pattern indicating that this analysis perhaps
* is with a too simple model. The QQ-plot looks fine.
contrast group@workload, nowald eff
* | Contrast Std. Err. z P>|z| [95% Conf. Interval]
* -------------------+---------------------------------------------------------
* co |
* group@workload |
* (CHF vs base) 0 | -2.366667 .8936099 -2.65 0.008 -4.11811 -.6152235
* (CHF vs base) 30 | -3.663861 .8985599 -4.08 0.000 -5.425006 -1.902716
* (CHF vs base) 60 | -3.675937 .8960896 -4.10 0.000 -5.43224 -1.919633
* (CHF vs base) 90 | -4.010868 .896087 -4.48 0.000 -5.767166 -2.254569
* (CHF vs base) 120 | -4.490822 .8990593 -5.00 0.000 -6.252946 -2.728698
* We see significant differences in cardiac output between CHF patients and
* healthy volunteers that seems to increase with increasing workload. At all
* workloads the healthy volunters have a higher mean cardiac output than the
* CHF-patients:
*
* Healthy - CHF: 0 2.4, 95%-CI: 0.6 - 4.1, p=0.008
* 30 3.7, 95%-CI: 1.9 - 5.4, p<0.0001
* 60 3.7, 95%-CI: 1.9 - 5.4, p<0.0001
* 90 4.0, 95%-CI: 2.3 - 5.8, p<0.0001
* 120 4.5, 95%-CI: 2.7 - 6.3, p<0.0001
* 4. Analyze the data using the multivariate repeated measurements ANOVA. Is it
* possible to simplify to the univariate repeated measurement model?
xi: mixed co bn.workload#bn.group, noconst || subj: i.workload, cov(un) ///
technique(bfgs)
estimates store mod1
xi: mixed co bn.workload bn.group || subj: i.workload, cov(un) technique(bfgs)
estimates store mod2
lrtest mod1 mod2
* Likelihood-ratio test LR chi2(4) = 19.45
* (Assumption: mod2 nested in mod1) Prob > chi2 = 0.0006
* Based on the multivariate repeated measurements ANOVA we make the same
* conclusion: We clearly reject the hypothesis of the same difference in mean
* cardiac output between healthy and CHF-patients at all time-points
* (p=0.0006).
xi: mixed co bn.workload#bn.group, noconst || subj: i.workload, cov(un) reml
predict fit1 if e(sample), fit
predict res1 if e(sample), rstandard
scatter res1 fit1, yline(0) name(q1,replace)
qnorm res1, title ("residual probability-plot") name(q2,replace)
graph combine q1 q2
* Looking at the plot of the residuals against the fitted values we would
* conclude that the patterns that we observed in the univariate analysis are
* not present here in the multivariate analysis. The QQ-plot is still fine.
contrast group@workload, nowald eff
* | Contrast Std. Err. z P>|z| [95% Conf. Interval]
* -------------------+---------------------------------------------------------
* co |
* group@workload |
* (CHF vs base) 0 | -2.366667 .79766 -2.97 0.003 -3.930052 -.8032819
* (CHF vs base) 30 | -3.639698 1.00533 -3.62 0.000 -5.610107 -1.669288
* (CHF vs base) 60 | -3.685215 .9164765 -4.02 0.000 -5.481476 -1.888954
* (CHF vs base) 90 | -4.001846 .9466934 -4.23 0.000 -5.85733 -2.146361
* (CHF vs base) 120 | -4.589373 .7986746 -5.75 0.000 -6.154746 -3.023999
* We get exactly the same conclusions as in the univariate analysis:
*
* Healthy - CHF: 0 2.4, 95%-CI: 0.8 - 3.9, p=0.003
* 30 3.6, 95%-CI: 1.7 - 5.6, p<0.0001
* 60 3.7, 95%-CI: 1.9 - 5.5, p<0.0001
* 90 4.0, 95%-CI: 2.1 - 5.9, p<0.0001
* 120 4.6, 95%-CI: 3.0 - 6.2, p<0.0001
* We conclude by testing whether we can reduce the multivariate repeated
* measurements ANOVA to the univariate:
xi: mixed co bn.workload#bn.group, noconst || subj: i.workload, cov(un) reml
estimates store multi
mixed co bn.workload#bn.group, noconst || subj: , reml
estimates store uni
lrtest multi uni
* Likelihood-ratio test LR chi2(14) = 47.20
* (Assumption: uni nested in multi) Prob > chi2 = 0.0000
* We get a clear rejection of the hypothesis corresponding to a reduction of
* the multivariate repeated measurements ANOVA to the univariate, p<0.0001.