For Lab Meeting 20180223

rmANOVA.do

/* 
    This script is based on the instruction provided by IDRE UCLA
    https://stats.idre.ucla.edu/stata/seminars/repeated-measures-analysis-with-stata/ 
*/


/**** Data ****/
use https://stats.idre.ucla.edu/stat/data/repeated_measures, clear
/* 
    There are a total of 
        8 subjects measured at 
        4 time points each. 
    These data are in wide format where 
        y1 is the response at time 1, 
        y2 is the response at time 2, and so on. 
    The subjects are divided into 
        2 groups of 4 subjects 
        using the the variable trt. 
*/

d, s  // describe, short
li, sep(4)  // list, separator(4)

sum y1-y4  // summarize y1-y4
bysort trt: sum y1-y4  // summary stats of each group (here 'trt')
tabstat y1-y4, by(trt) stat(n mean sd var)  // this is more flexible

/* A graphical summary of this dataset, using a user package */
net install profileplot.pkg
profileplot y1-y4, by(trt)  // plot the means, at each time point
profileplot y1-y4, by(trt) median  // medians

/* Let's check the correlation and covariance matrix */
corr y1-y4
corr y1-y4, cov  // correlate y1-y4, cov
/* 
    For checking whether if covariance matrix have compound symmetry (exchangable).
        _             _
        |  1.           |
        |  ρ  1.        |
        |  ρ  ρ  1.     |
        |_ ρ  ρ  ρ  1. _|
    If not, we shall check the corrected F-tests. 
*/

/* Reshape the data into long format. */
reshape long y, i(id) j(time)
/* 
	:long: {'long' or 'wide'} convert to long/wide form.
	:y: prefixs or variables list (depends on data shape) of the columns that stores measurements
	:i(): (variable) for identifying subjects
	:j(): (variable) for identifying time point
*/

/* Check the data again */
li, sep(4)


/**** Repeated Measures ANOVA ****/
/* 
    In general, rule of thumb is that there is 
        one single error term for all of the between-subject effects and 
        a separate error term for each of the within-subject factors and 
            for the interaction of within-subject factors.
*/

anova y trt / id|trt time trt#time, repeated(time)
/* 
	:y: variable of outcome measurement
	:trt: main effect from treatment
	: / id|trt: slash indicates the id|trt is the error term for 'trt';
                pipe indcates 'id' is nested in 'trt'
	:time: main effect from time
	:trt#time: interaction term
	:repeated(): (variable) indicate repeated measure, at different 'time'
*/

/* Check the pooled within-subject covariance matrix */
matrix list e(Srep)  // pooled within-subject covariance matrix


/**** Test of Simple Effects ****/

/* The effect of time at each treatment */
contrast time@trt, effect
margins time, at(trt=1) pwcompare(effects) noestimcheck  // pairwise ones
margins time, at(trt=2) pwcompare(effects) noestimcheck  
/*
	:time: variable of interest
	:at(): (condition) to be estimate (trt=1 here)
	:pwcompare(effects): option that indicates pairwise comparison, 
                         'effects' option shows confidence and p-values
	:noestimcheck: option that suppress estimability checks
*/

/* The effects of treatment at each time point */ 
/* 
    This is estimated with pooled error, recall the command before:
    anova y trt / id|trt time trt#time, repeated(time)
*/
anova y trt time trt#time  // shorthand: anova y trt##time
contrast trt@time, effect
margins trt#time
marginsplot, x(time)


/**** An Alternative: Mixed Model ****/
xtmixed y trt##time || id:, var reml
/*
    :y: variable of outcome measurements
    :trt##time: trt, time, and interaction, using ## expansion is preferred here
    :id: variable for identifying individual
    :var: option to request covariance be reported
    :reml: option to request restricted-ML, to produce resuls similiar to ANOVA
*/


/**** End ****/