re: st: Re: Contrasts and interpretation

[David Airey <david.airey@Vanderbilt.Edu>]

Hi Jaye,

Here is a worked example I had on file on a 3 way factorial analysis  
with posthoc tests for the 3 way and 2 way significant interactions  
using the cell means anova. I have some gripes with Stata making it a  
little difficult to do this kind of thing. For example, the command  
test does not make it easy to test contrasts using an error term  
other than the model error, which in any moderately complicated mixed  
model will need doing. Compare what you have to do below versus what  
you have to do in JMP or SPSS, for example. Hope this helps.

// Table 8.12 from ISBN0-8058-3718-3 but labeled differently
// This text only covers solutions by SPSS and SAS, so here is Stata
clear
input y a b c
170 1 2 1
175 1 2 1
165 1 2 1
180 1 2 1
160 1 2 1
158 1 2 1
161 1 2 2
173 1 2 2
157 1 2 2
152 1 2 2
181 1 2 2
190 1 2 2
186 2 2 1
194 2 2 1
201 2 2 1
215 2 2 1
219 2 2 1
209 2 2 1
164 2 2 2
166 2 2 2
159 2 2 2
182 2 2 2
187 2 2 2
174 2 2 2
180 3 2 1
187 3 2 1
199 3 2 1
170 3 2 1
204 3 2 1
194 3 2 1
162 3 2 2
184 3 2 2
183 3 2 2
156 3 2 2
180 3 2 2
173 3 2 2
173 1 1 1
194 1 1 1
197 1 1 1
190 1 1 1
176 1 1 1
198 1 1 1
164 1 1 2
190 1 1 2
169 1 1 2
164 1 1 2
176 1 1 2
175 1 1 2
189 2 1 1
194 2 1 1
217 2 1 1
206 2 1 1
199 2 1 1
195 2 1 1
171 2 1 2
173 2 1 2
196 2 1 2
199 2 1 2
180 2 1 2
203 2 1 2
202 3 1 1
228 3 1 1
190 3 1 1
206 3 1 1
224 3 1 1
204 3 1 1
205 3 1 2
199 3 1 2
170 3 1 2
160 3 1 2
179 3 1 2
179 3 1 2
end

// 3 x 2 x 2 factorial design
table a b c, contents(count y mean y sd y)

// overparameterized anova model
anova y a b c a*b a*c b*c a*b*c
anovaplot

// cell means anova model
egen z = group(a b c)
// this table is key to setting up all matrices below so refer back
table a b c, contents(mean z) // table rowvar [colvar [supercolvar]]
anova y z, noconstant
regress, noheader

// 6 omnibus effects from the cell means anova are exactly
// the same as from the overparameterized model above

mat a = (1, 1, 1, 1, 0, 0, 0, 0, -1, -1, -1, -1 \  ///
         0, 0, 0, 0, 1, 1, 1, 1, -1, -1, -1, -1 \  ///
         1, 1, 1, 1, -1, -1, -1, -1, 0, 0, 0, 0)

mat b = (1,  1, -1, -1, 1,  1, -1, -1, 1, 1, -1, -1)

mat c = (1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1)

// element by element multiplication of matrix a and b
forvalues i = 1/3 {
	forvalues j = 1/1 {
		mat axb = nullmat(axb) \ vecdiag(a[`i',1...]'*b[`j',1...])
	}
}
mat list axb

// element by element multiplication of matrix a and c
forvalues i = 1/3 {
	forvalues j = 1/1 {
		mat axc = nullmat(axc) \ vecdiag(a[`i',1...]'*c[`j',1...])
	}
}
mat list axc

// element by element multiplication of matrix b and c
forvalues i = 1/1 {
	forvalues j = 1/1 {
		mat bxc = nullmat(bxc) \ vecdiag(b[`i',1...]'*c[`j',1...])
	}
}
mat list bxc

// element by element multiplication of matrix a and bxc
forvalues i = 1/3 {
	forvalues j = 1/1 {
		mat axbxc = nullmat(axbxc) \ vecdiag(a[`i',1...]'*bxc[`j',1...])
	}
}
mat list axbxc

test, test(a)
test, test(b)
test, test(c)
test, test(axb)
test, test(axc)
test, test(bxc)
test, test(axbxc)

// Since the 3 way interaction is significant,
// we can look at AB interaction at each level of C.

mat c1 = (1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0)
mat c2 = (0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1)

forvalues i = 1/3 {
	forvalues j = 1/1 {
		mat axb_c1 = nullmat(axb_c1) \ vecdiag(axb[`i',1...]'*c1[`j',1...])
		mat axb_c2 = nullmat(axb_c2) \ vecdiag(axb[`i',1...]'*c2[`j',1...])
	}
}

mat list axb_c1
mat list axb_c2

test, test(axb_c1)
test, test(axb_c2)

// Since a*b is significant when c = 1
// but not when c = 2, subsequent post-hoc tests in each
// will be different. When c = 1, we look at "simple,
// simple main effects", or b with a1c1, b within a2c1, b within a3c1,
// as well as a within b1c1 and a within b2c1.
// These comparisons are called "simple, simple" because
// two factors are fixed at one level.

mat b1c1 = (1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0)
mat b2c1 = (0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0)

mat a1c1 = (1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0)
mat a2c1 = (0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0)
mat a3c1 = (0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0)

forvalues i = 1/3 {
	forvalues j = 1/1 {
		mat a_b1c1 = nullmat(a_b1c1) \ vecdiag(a[`i',1...]'*b1c1[`j',1...])
		mat a_b2c1 = nullmat(a_b2c1) \ vecdiag(a[`i',1...]'*b2c1[`j',1...])
	}
}

forvalues i = 1/1 {
	forvalues j = 1/1 {
		mat b_a1c1 = nullmat(b_a1c1) \ vecdiag(b[`i',1...]'*a1c1[`j',1...])
		mat b_a2c1 = nullmat(b_a2c1) \ vecdiag(b[`i',1...]'*a2c1[`j',1...])
		mat b_a3c1 = nullmat(b_a3c1) \ vecdiag(b[`i',1...]'*a3c1[`j',1...])
	}
}

test, test(a_b1c1)
test, test(a_b2c1)
test, test(b_a1c1)
test, test(b_a2c1)
test, test(b_a3c1)

// Now, since b has two levels we need not test further.
// We should test further for a_b1c1 and a_b2c1, since a has 3 levels

mat a12_b1c1 = (1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0)
mat a13_b1c1 = (1, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0)
mat a23_b1c1 = (0, 0, 0, 0, 1, 0, 0, 0, -1, 0, 0, 0)

mat a12_b2c1 = (0, 0, 1, 0, 0, 0, -1, 0, 0, 0, 0, 0)
mat a13_b2c1 = (0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, 0)
mat a23_b2c1 = (0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0)

test, test(a12_b1c1)
test, test(a13_b1c1)
test, test(a23_b1c1)
test, test(a12_b2c1)
test, test(a13_b2c1)
test, test(a23_b2c1)

// Now, we go back to the a*b layout when c = 2. a*b was
// not significant, so we test a_c2 and b_c2:

mat list a
mat list b
mat list c2

forvalues i = 1/3 {
	forvalues j = 1/1 {
		mat a_c2 = nullmat(a_c2) \ vecdiag(a[`i',1...]'*c2[`j',1...])
	}
}

forvalues i = 1/1 {
	forvalues j = 1/1 {
		mat b_c2 = nullmat(b_c2) \ vecdiag(b[`i',1...]'*c2[`j',1...])
	}
}

test, test(a_c2)
test, test(b_c2) // further testing is not necessary since b has 2  
levels



*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/