wald test in heterogenous variance model using proc glimmix

[Tomchoi]

when I run heterogeneous variance model using proc glimmix to test homogeneity of variance assumption I get this

DATA a;
input id 1 batch 1 ref 5-7 new 9-11;
diff=new-ref;
datalines;
1 116 119
2 115 121
3 103 121
4 135 130
5 108 116
6 113 122
7 115 117
8 127 124
;
keep batch ref new diff;

*wide-format data to long-format data;
DATA b;
SET a;
ARRAY trials[2] ref new;
DO TIME = 1 TO 2;
y = trials[TIME];
trial = TIME-1;
OUTPUT; *wihtout this time=3 values overwrite all;
END;
KEEP batch y trial;
proc print data=c (OBS=12) NOOBS; run;

*run glm;
PROC GLIMMIX data=b plots=residualpanel;
CLASS trial (ref=first);
MODEL y=trial /solution ddfm=KR;
random _residual_ /group=trial type=VC;
covtest /wald;
LSMEANS trial /alpha=0.05 cl diff plot=diffplot;
run;

Covariance Parameter Estimates

Cov Parm	Group	Estimate	Standard Error	Z Value	Pr > Z
Residual (VC)	trial 1	19.3571	10.3468	1.87	0.0307
Residual (VC)	trial 0	103.43	55.2849	1.87	0.0307

my question is: why z and p values for variance components are fixed to be equal to each other? is this some sort of model constraints? If the sample do not pass the homogeneity variance test so I end up keeping heterogenous variance in the model, can i use and interpret these p-values? Thanks!

Best,

Tom

[SteveDenham]

Hi @Tomchoi ,

I can't reproduce what you have, given your set-up code (The GLIMMIX part looks OK for now).  When I run everything down to the PROC PRINT, the dataset WORK.b is missing the y variable for all instances of Trial 0.  Could you double check what is going on?  Once we have values for Trial 0, then we can see what GLIMMIX is up to.

SteveDenham

[Tomchoi]

sorry there was mistake with data procedure. pls use this code:

DATA a;
input batch 1 ref 3-5 new 6-9;
diff=new-ref;
datalines;
1 116 119
2 115 121
3 103 121
4 135 130
5 108 116
6 113 122
7 115 117
8 127 124
;
proc print data=a; run;

[Tomchoi]

this is full code:

/*Entering data manually*/
DATA a;
input batch 1 ref 3-5 new 6-9;
diff=new-ref;
datalines;
1 116 119
2 115 121
3 103 121
4 135 130
5 108 116
6 113 122
7 115 117
8 127 124
;
proc print data=a; run;

*wide-format data to long-format data;
DATA b;
SET a;
ARRAY trials[2] ref new;
DO TIME = 1 TO 2;
y = trials[TIME];
trial = TIME-1;
OUTPUT; *wihtout this time=3 values overwrite all;
END;
KEEP batch y trial;
proc print data=b (OBS=12) NOOBS; run;

*run glm;
PROC GLIMMIX data=b plots=residualpanel;
CLASS trial (ref=first);
MODEL y=trial /solution ddfm=KR;
random _residual_ /group=trial type=VC;
covtest /wald;
LSMEANS trial /alpha=0.05 cl diff plot=diffplot;
run;

[SteveDenham]

Thanks, @Tomchoi .

I have poked and prodded at this, and it seems that no matter how I reparameterize, the Wald tests come out the same.  All I can guess is that the denominators (asymptotic standard errors) are wholly dependent on the estimate, so that algebraically you get the same ratio, and thus the same Z value.

SteveDenham

[Tomchoi]

Hi. Thanks for your honest answer 🙂 do you still think I can report & interpret the t values and p values?

[SteveDenham · Accepted solution]

The fixed effect t values and p values certainly should be fine. I have some hesitation about the Wald chi-square and p values for the random effects with only 8 subjects, as these are asymptotically correct as the sample size approaches infinity, and while 8 is a large number for tests about fixed effects, it certainly is not for tests about variance components.

SteveDenham