PROC MIXED for ANCOVA: Comparing Slopes for Each Group

[jolich1]

Hi, I'm very new to SAS, and I'm trying to get a means separation test done for the slopes of an ANCOVA, but I'm unclear how to do so. 

Sample data is attached as a csv.

My data has the following variables: Height_cm (Continuous quantitative), Diff (Continuous quantitative) Type (a factor variable with 6 levels), Treatment (a factor variable with 2 levels), State (a random effect), blockid (blocking factor), groupid (blocking factor nested in blockid).

For my ANCOVA, I am running a regression of  Diff = a + b(Height_cm)

My code for extracting estimates at a value of 25 currently looks like this:

proc mixed empirical data=sampledata;
    class Type Treatment blockid groupid State;
    model Diff = Height_cm|Type|Treatment / solution;
    random State blockid groupid(blockid) / subject = groupid;
    lsmeans Type*Treatment / at Height_cm = 25 adjust=tukey;
    store BlockAnalysis / label='PLM: Getting Started';
run;
ods trace on;
proc plm restore=blockanalysis;
   lsmeans Type*Treatment / at Height_cm = 25 adjust=tukey lines;
   ods output LSMLines=LSMLines;
run;
ods trace off

Ideally I would like code that produces the same tukey lines output, but in comparison of regression slopes for each grouping of Type*Treatment

For those familiar with R, this would replicate the use of lmer, emmeans, and cld commands

model = lmer(Diff ~ Height_cm * Type * Treatment + (1|State) + (1|blockid) +(1|blockid:groupid), data = sampledata)
marginal = emtrends(model,pairwise ~ Type * Treatment, var = "Height_cm")
CLD = cld(marginal, alpha = 0.05, Letters = letters, adjust = "tukey")
CLD

[Ksharp]

Your MODEL statement did not looks like you are fitting a ANCOVA model.

And your RANDOM statement is messy and unreadable.Should like:

random int State blockid/ subject = groupid;

And check this to compare the slope between groups:

https://support.sas.com/kb/38/384.html

https://support.sas.com/kb/24/177.html

https://support.sas.com/kb/70/756.html

[jolich1]

Huge thanks to Ksharp for the direction. Posting as a reply to help anyone else that stumbles on this post.

A helpful link I found for understanding linear combinations of model parameters that is used in the estimate statement: https://stats.oarc.ucla.edu/sas/faq/how-do-i-write-an-estimate-statement-in-proc-glm/

Instead of proc mixed, I've used proc glimmix, which allows me to provide an adjustment for multiple tests.

proc glimmix empirical data=iheight;
    class Type Treatment blockid groupid State;
    model Diff = Height_cm|Type|Treatment / solution;
    random State blockid groupid(blockid) / subject = groupid;
    lsmeans Type*Treatment / at Height_cm = 25 adjust=tukey;
    estimate 'a open' Height_cm 1 Height_cm*Treatment 1 0 Height_cm*Type 1 0 0 0 0 0 Height_cm*Type*Treatment 1 0 0 0 0 0 0 0 0 0 0 0,
    'b open' Height_cm 1 Height_cm*Treatment 1 0 Height_cm*Type 0 1 0 0 0 0 Height_cm*Type*Treatment 0 0 1 0 0 0 0 0 0 0 0 0,
    'c open'Height_cm 1 Height_cm*Treatment 1 0 Height_cm*Type 0 0 1 0 0 0 Height_cm*Type*Treatment 0 0 0 0 1 0 0 0 0 0 0 0,
    'd open'Height_cm 1 Height_cm*Treatment 1 0 Height_cm*Type 0 0 0 1 0 0 Height_cm*Type*Treatment 0 0 0 0 0 0 1 0 0 0 0 0,
    'e open'Height_cm 1 Height_cm*Treatment 1 0 Height_cm*Type 0 0 0 0 1 0 Height_cm*Type*Treatment 0 0 0 0 0 0 0 0 1 0 0 0,
    'f open'Height_cm 1 Height_cm*Treatment 1 0 Height_cm*Type 0 0 0 0 0 1 Height_cm*Type*Treatment 0 0 0 0 0 0 0 0 0 0 1 0,
    'a closed' Height_cm 1 Height_cm*Treatment 0 1 Height_cm*Type 1 0 0 0 0 0 Height_cm*Type*Treatment 0 1 0 0 0 0 0 0 0 0 0 0,
    'b closed' Height_cm 1 Height_cm*Treatment 0 1 Height_cm*Type 0 1 0 0 0 0 Height_cm*Type*Treatment 0 0 0 1 0 0 0 0 0 0 0 0,
    'c closed'Height_cm 1 Height_cm*Treatment 0 1 Height_cm*Type 0 0 1 0 0 0 Height_cm*Type*Treatment 0 0 0 0 0 1 0 0 0 0 0 0,
    'd closed'Height_cm 1 Height_cm*Treatment 0 1 Height_cm*Type 0 0 0 1 0 0 Height_cm*Type*Treatment 0 0 0 0 0 0 0 1 0 0 0 0,
    'e closed'Height_cm 1 Height_cm*Treatment 0 1 Height_cm*Type 0 0 0 0 1 0 Height_cm*Type*Treatment 0 0 0 0 0 0 0 0 0 1 0 0,
    'f closed'Height_cm 1 Height_cm*Treatment 0 1 Height_cm*Type 0 0 0 0 0 1 Height_cm*Type*Treatment 0 0 0 0 0 0 0 0 0 0 0 1 / adjust = BON;
    store BlockAnalysis / label='PLM: Getting Started';
run;

I can then plot the estimates and standard errors to compare means separations (similar to a lines command)

As for if my model is correct, still unsure. My model statement (I hope) is testing a three-way interaction (Factorial ANCOVA, creating a regression of Diff by Height_cm, then comparing intercepts and slopes of subsequent combinations of Type and Treatment.)

[SteveDenham]

On your ESTIMATE statements, I assume you wish to find marginal values over Type, Treatment and the Type by Treatment interaction, as there are no terms included in the statement for those factors.

Have you considered using LSMESTIMATE statements, with the AT= option for various interesting values for the covariate?

SteveDenham