Repeated measures over space and time using hpmixed and glimmix

[TD21]

DATA Summary: Blocks = 3, Treatments (Treat) per block = 3, Locations (LOC) per treatment = 2, Depths (Depth) per LOC per Treat = 3. Date collected monthly on a given date over three years.

So, basically I have repeated measurements over space and time. I am using HPMIXED and GLIMMIX combination to minimize running time and memory issues. I am thinking of averaging over LOC and Depth and do the analysis with repeated measures overtime only. However, I would appreciate it if you could share ideas on how to appropriately do the analysis. Thank you.

DATA Work.Log_WQ2013;
SET Work.WQ2013;
LOG_WQ1 = log10(WQ1+0.01);
RUN;

PROC HPMIXED DATA = Work.Log_WQ2013;
CLASS Block Treat LOC Depth Year Season Date;
MODEL WQ1 = Treat|LOC|Depth|Year|Season/DIST = Lognormal DDFM = KR2;
RANDOM INT Treat/SUBJECT = Block;

REPEATED Date/SUBJECT = Block*Treat*LOC*Depth*Year*Seas TYPE = ar(1);
TEST Treat|LOC|Depth|Year|Seas;
ODS OUTPUT CovParms = HPMEstimate;
RUN;

%LET DEPVAR = WQ1;
PROC GLIMMIX DATA = Work.Log_WQ2013 method=rmpl noprofile;
CLASS Block Treat LOC Depth Year Season Date;
MODEL WQ1 = Treat|LOC|Depth|Year|Season/DIST = Lognormal DDFM = KR2;
RANDOM INT Treat/SUBJECT = Block;

RANDOM Date/SUBJECT = Block*Treat*Well*Depth*Year*Seas TYPE = ar(1);
PARMS/PDATA=HPMEstimate  noiter;
NLOPTIONS  technique=none;
LSMEANS Treat*Year Treat*Seas Treat*Year*Seas/ILINK ADJUST = TUKEY
SLICEDIFF = (Treat Year Seas)SLICEDIFFTYPE = ALL CL;
ODS OUTPUT LSMEANS = LSMNDS_&DEPVAR COVPARMS = COVPARMS_&DEPVAR;

RUN;

[lvm]

You will need to read more about HPMIXED. It is strictly for linear, normal, mixed models. Thus, there is no /DIST= option for a different distribution. Also, DDFM=KR or KR2 is not supported. Only residual df, or manually specified numbers. There are no easy sparse-matrix ways of doing more complex df calculations. Since you are assuming lognormal, you could model log_wq1 as the response variable in HPMIXED (and then pass the results to GLIMMIX, where you would be back to wq1 as the response). But note, you are adding 0.01 in the data step, presumably to avoid zeros. But if there are zeros, then your data cannot have a log-normal distribution (by definition).

With log-normal, there is a free scale parameter. Thus, in glimmix, your last random statment should have a residual option (after the /). Reason: with your code, you would be getting a residual nested under the repeated measure. Normally you want the residual term to account for the serial correlation. There may be other things you need to do also, but these are quick comments. By the way, based on the given information, I doubt if you would need HPMIXED, although it might help some for speed.

[SteveDenham]

What Larry said.

Steve Denham