topic Proc glimmix contrasts in Statistical Procedures

Proc glimmix contrasts

GavinB — Thu, 12 Oct 2017 10:08:07 GMT

Dear Reader,

In my model, where treatment is fixed and block is random, there seems to be a different treatment effect for a different set of blocks. Attached you see code where I estimate treatment effects per set of random blocks, and code where I estimate the difference for a set of blocks. The estimates do NOT coincide, meaning that the difference per treatment does not flow from the estimates of a treatment. Any help on explaining this difference would be much appreciated.

proc glimmix data=test_Per plots=all;
class TT block;
model Per=TT / s cl ddfm=kenwardroger;
random block / s cl g gcorr;
estimate 	'trt=1 @ blocks 1-6' 		intercept 6 TT  6   0	| block 1 1 1 1 1 1 0 0 0 0 0 0 / divisor=6;
estimate 	'trt=2 @ blocks 1-6' 		intercept 6 TT  0 	6 	| block 1 1 1 1 1 1 0 0 0 0 0 0 / divisor=6;
estimate 	'trt=1 vs trt=2 @ blocks 1-6' 		TT -6 	6 	0	0	0	0	0	0	0	0  	| block 1 1 1 1 1 1 0 0 0 0 0 0, 
			'trt=1 vs trt=3 @ blocks 1-6' 		TT -6 	0 	6	0	0	0	0	0	0  	0	| block 1 1 1 1 1 1 0 0 0 0 0 0,
			'trt=1 vs trt=4 @ blocks 1-6' 		TT -6 	0 	0	6	0	0	0	0	0  	0	| block 1 1 1 1 1 1 0 0 0 0 0 0,
			'trt=1 vs trt=5 @ blocks 1-6' 		TT -6 	0  	0	0	6	0	0	0	0	0	| block 1 1 1 1 1 1 0 0 0 0 0 0,	
			'trt=1 vs trt=6 @ blocks 1-6' 		TT -6  	0  	0	0	0	6	0	0	0	0	| block 1 1 1 1 1 1 0 0 0 0 0 0,
			'trt=1 vs trt=7 @ blocks 1-6' 		TT -6  	0 	0	0	0	0	6	0	0 	0	| block 1 1 1 1 1 1 0 0 0 0 0 0,
			'trt=1 vs trt=8 @ blocks 1-6' 		TT -6  	0  	0	0	0	0	0	6	0	0	| block 1 1 1 1 1 1 0 0 0 0 0 0, 	
			'trt=1 vs trt=9 @ blocks 1-6' 		TT -6  	0  	0	0	0	0	0	0	6	0	| block 1 1 1 1 1 1 0 0 0 0 0 0,	
			'trt=1 vs trt=10 @ blocks 1-6' 		TT -6  	0  	0	0	0	0	0	0	0	6	| block 1 1 1 1 1 1 0 0 0 0 0 0 /	divisor=6 e adjdfe=row adjust=simulate(nsamp=100000);
estimate 	'trt=1 vs trt=2 @ blocks 7-12' 		TT -6 	6 	0	0	0	0	0	0	0	0	| block 0 0 0 0 0 0 1 1 1 1 1 1,  	
			'trt=1 vs trt=3 @ blocks 7-12' 		TT -6 	0 	6	0	0	0	0	0	0 	0	| block 0 0 0 0 0 0 1 1 1 1 1 1, 	
			'trt=1 vs trt=4 @ blocks 7-12' 		TT -6 	0   0	6	0	0	0	0	0 	0	| block 0 0 0 0 0 0 1 1 1 1 1 1,	
			'trt=1 vs trt=5 @ blocks 7-12' 		TT -6  	0	0	0	6	0	0	0	0	0	| block 0 0 0 0 0 0 1 1 1 1 1 1, 	
			'trt=1 vs trt=6 @ blocks 7-12' 		TT -6  	0	0	0	0	6	0	0	0	0	| block 0 0 0 0 0 0 1 1 1 1 1 1, 	
			'trt=1 vs trt=7 @ blocks 7-12' 		TT -6  	0	0	0	0	0	6	0	0	0	| block 0 0 0 0 0 0 1 1 1 1 1 1, 	
			'trt=1 vs trt=8 @ blocks 7-12' 		TT -6  	0	0	0	0	0	0	6	0	0	| block 0 0 0 0 0 0 1 1 1 1 1 1, 		
			'trt=1 vs trt=9 @ blocks 7-12' 		TT -6  	0	0	0	0	0	0	0	6	0	| block 0 0 0 0 0 0 1 1 1 1 1 1, 
			'trt=1 vs trt=10 @ blocks 7-12' 	TT -6  	0	0	0	0	0	0	0	0	6	| block 0 0 0 0 0 0 1 1 1 1 1 1 /	divisor=6	e adjdfe=row adjust=simulate(nsamp=100000);
store glimmix_Per;
run;

Kind regards,

Gavin

Re: Proc glimmix contrasts

sld — Tue, 07 Nov 2017 04:23:19 GMT

I looked at this not too long after you posted, but didn't see any obvious problem. I've been hoping someone (else) would respond.

If you are inclined, try posting again, including your data set or some variant that produces the same problem. Maybe you'll get a nibble.

If you have not already, you could consult

https://www.crcpress.com/Generalized-Linear-Mixed-Models-Modern-Concepts-Methods-and-Applications/Stroup/p/book/9781439815120

Edit: Or maybe you solved it already! Then you can post your own solution 🙂

Re: Proc glimmix contrasts

GavinB — Tue, 07 Nov 2017 07:26:06 GMT

Well, we managed to figure out what the problem was after quite some time (and a lot of headaches).

The way the code was setup was simply not correct to show the results that I wanted. The design was a randomized complete block design and because of the way a mixed model works (using the entire dataset to estimate variance) – the difference between treatments will always be the same in each of the blocks if you code it like I did above. In order to actually see if treatment differences change due to the positioning in (for this particular dataset) barn, you need to include an interaction effect. This will give the replicated estimates of the treatments that I wanted.

That basically solves the issue