BookmarkSubscribeRSS Feed
🔒 This topic is solved and locked. Need further help from the community? Please sign in and ask a new question.
lcgaras
Fluorite | Level 6

Hello,

 

I'm relatively unexperienced with SAS and statistics in general, but I've managed to figure out how to analyze most of my data until I got to the data set that include repeated samples over time from the same pigs.  

 

Study design includes 4 treatments, with 3 measurements taken over time from several pigs in each treatment (it's also unbalanced). I have two questions that I'm trying to answer with this analysis:

 

1) Is there a change over time within each of the treatments?

2) Is there a difference between the 4 treatments at t0? t1? t2?

 

I've tried to include "residual" in a Random statement in GLIMMIX (based on the equivalent in GLM), and I'm still struggling.  Here's what I've done so far:

 

Data DSS;
Title 'DSS DSS-hLZ RBC';
Input pig $ Trt $ Time Resp;
Cards;
28-5	TRT1	0	7.08
28-8	TRT1	0	7.34
29-12	TRT1	0	8.09
29-14	TRT1	0	7.98
30-4	TRT1	0	7.34
30-4	TRT1	0	8.07
30-7	TRT1	0	8.41
28-6	TRT2	0	7
29-10	TRT2	0	6.9
30-6	TRT2	0	8.8
31-1	TRT2	0	7.08
28-1	TRT3	0	7.61
28-7	TRT3	0	7.57
30-5	TRT3	0	7.23
31-6	TRT3	0	6.6
28-3	CTRL	0	7.16
28-9	CTRL	0	7.54
30-3	CTRL	0	6.66
31-4	CTRL	0	8.28
31-5	CTRL	0	6.68
28-5	TRT1	1	7.11
28-8	TRT1	1	8.16
29-12	TRT1	1	8.25
29-14	TRT1	1	7.47
30-4	TRT1	1	7.86
30-4	TRT1	1	8.06
30-7	TRT1	1	7.89
28-6	TRT2	1	7.17
29-10	TRT2	1	7.54
30-6	TRT2	1	9.75
31-1	TRT2	1	7.66
28-1	TRT3	1	7.6
28-7	TRT3	1	7.26
30-5	TRT3	1	7.45
31-6	TRT3	1	6.4
28-3	CTRL	1	7.37
28-9	CTRL	1	7.84
30-3	CTRL	1	7
31-4	CTRL	1	8.08
31-5	CTRL	1	7.18
28-5	TRT1	2	7.09
28-8	TRT1	2	8.19
29-12	TRT1	2	8.28
29-14	TRT1	2	7.25
30-4	TRT1	2	6.67
30-4	TRT1	2	6.7
30-7	TRT1	2	7.73
28-6	TRT2	2	7.49
29-10	TRT2	2	7.05
30-6	TRT2	2	7.89
31-1	TRT2	2	6.47
28-1	TRT3	2	7.27
28-7	TRT3	2	7
30-5	TRT3	2	7.61
31-6	TRT3	2	6.38
28-3	CTRL	2	6.81
28-9	CTRL	2	7.54
30-3	CTRL	2	6.32
31-4	CTRL	2	7.06
31-5	CTRL	2	6.84
;
Proc GLIMMIX Data=DSS;
Class pig trt time;
Model Resp = trt(time);
Random time / subject=pig type=ar(1) residual;
LSMEANS trt(time)/ diff lines;
Run; 
Quit;

The problem is when I include the "residual" in the Random statement, it gives me an error: "WARNING: Obtaining minimum variance quadratic unbiased estimates as starting values for the covariance parameters failed."

 

When I remove "residual" I get every pairwise comparisons, most of which are irrelevant (ex: Ctrl t0 vs Trt3 t2).

Proc GLIMMIX Data=CBC;
Class pig trt time;
Model Resp = trt(time);
Random time / subject=pig type=ar(1);
LSMEANS trt(time)/ diff lines;
Run; 
Quit;

If I change the Model (and LSMEANS) to just "trt" or "time" (or "trt time trt*time") I still don't get the correct comparisons.

 

This is just one of about 20 such variables, so once I figure it out, I can move forward!

 

Thanks in advance, any help is much appreciated!

1 ACCEPTED SOLUTION

Accepted Solutions
lvm
Rhodochrosite | Level 12 lvm
Rhodochrosite | Level 12

You definitely want the "residual" in the model statement: this gives you the correlated R-side pattern that would be appropriate for your dataset. That warning message can be caused by many things, and can be data-dependent. Your model is not apparently over-parameterized, so some of the common causes would not apply. For a similar problem (reported elsewhere on the web), the issue was the detection of singularity in solving the mixed-model equations. One solution from SAS was to add the following option to the GLIMMIX statement:

singular=1e-9

As in proc glimmix data=.... singular=1e-9;

 

Another fix is to switch to PROC MIXED. There you use REPEATED statement instead of the random statement.

 

It seems that you have a different question related to the multiple comparisons of means. To get only the contrasts you want, you could use ESTIMATE, CONTRAST, or newer LSMESTIMATE statements. There is a lot on the web about learning these. In glimmix, you could also use SLICE command:

slice time*trt / sliceby=time sliceby=trt lines;

This gives comparisons of treatments at each level of time, and comparisons of times at each level of treatment.

 

 

View solution in original post

5 REPLIES 5
lvm
Rhodochrosite | Level 12 lvm
Rhodochrosite | Level 12

You definitely want the "residual" in the model statement: this gives you the correlated R-side pattern that would be appropriate for your dataset. That warning message can be caused by many things, and can be data-dependent. Your model is not apparently over-parameterized, so some of the common causes would not apply. For a similar problem (reported elsewhere on the web), the issue was the detection of singularity in solving the mixed-model equations. One solution from SAS was to add the following option to the GLIMMIX statement:

singular=1e-9

As in proc glimmix data=.... singular=1e-9;

 

Another fix is to switch to PROC MIXED. There you use REPEATED statement instead of the random statement.

 

It seems that you have a different question related to the multiple comparisons of means. To get only the contrasts you want, you could use ESTIMATE, CONTRAST, or newer LSMESTIMATE statements. There is a lot on the web about learning these. In glimmix, you could also use SLICE command:

slice time*trt / sliceby=time sliceby=trt lines;

This gives comparisons of treatments at each level of time, and comparisons of times at each level of treatment.

 

 

lcgaras
Fluorite | Level 6

Your suggestions solved *part* of my problem! Unfortunately, it didn't solve the error when I add in:

Proc GLIMMIX Data=DSS singular=1e-9;

I still can't run the code with "residual" in the Random statement without getting the same error mentioned in my original post; however, with some tweaking of the Model statement, and including the Slice command, I was able to get the analyses I was looking for (over time w/in a treatment, and at each time point b/w treatments):

Proc GLIMMIX Data=DSS singular=1e-9;
Class pig trt time;
Model Resp = trt time trt*time;
Random time / subject=pig type=ar(1);
slice time*trt / sliceby=time sliceby=trt lines;
LSMEANS trt time/ diff lines;
Run; 

So, that still leaves me with maybe not the most kosher analysis of data. Also, I can't remember why I switched from Proc Mixed (where I tried to use Repeated) to Glimmix, but I think it had to deal with the unbalanced design + multiple treatments + subsamples (for a different assay), and I can see the statistical groupings using the "lines" command.

 

Any other insight/suggestions for using the "residual" term w/o generating an error?

 

Thanks again!

lvm
Rhodochrosite | Level 12 lvm
Rhodochrosite | Level 12

You have to remember, you are now fitting a different model from what you intend. Without the residual option, you are getting the auto-regression correlation structure AND a residual, which would be analogous to a nugget effect that is available in MIXED. By the way, you said you got a WARNING, not error in the log. With a warning, I assume you still got output. Or not??? Warnings do not always mean stop.

 

You should try MIXED, just to see if it works.

Proc MIXED Data=DSS ;
Class pig trt time;
Model Resp = trt time trt*time;
repeated time / subject=pig type=ar(1);
LSMEANS trt time/ diff lines;
store results;
Run;
proc plm restore=results;
slice time*trt / sliceby=time sliceby=trt lines;
run;

As demonstrated here, you can use PROC PLM to get the slices. It is possible that SLICE has been added to MIXED in the most recent version which I don't have (I haven't looked). By the way, the slices I demonstrate are used when there is an interaction.

 

Also, be careful with your coding. You code would be appropriate if pig index crosses all treatments. For instance, pigs 1-5 for trt 1, pigs 6-10 for trt 2, piges 11-120 for trt 3, etc. However, If you use the same pig indices for each treatment (first pig of first trt, first pig of second trt), then your glimmix code is wrong, and that would generate an error (usually an infinite likelihood error). Incorrect coding would be pigs 1-5 for trt 1, pigs 1-5 for trt 2, etc. If you did the latter, change statement to:

random time / subject=pig(trt) type=ar(1) residual;

 

lcgaras
Fluorite | Level 6

Well, I kind of feel like an idiot right now, and just deleted everything I had previously typed in response...

 

While I was typing this up: "Just to clarify the data, each "pig" is an individual animal (the number is a unique identifier), and each pig is only in one treatment, it never received more than one treatment."

 

And that's when I realized that one of the notes that I got in the log while trying your MIXED / Repeated code said this:

 

"NOTE: An infinite likelihood is assumed in iteration 0 because of a nonpositive definite estimated R matrix for pig 30-4."

 

Which is when I realized that TWO pigs were labeled as "30-4", so it registered two t0, t1, t2 for that animal.  I fixed that, deleted the line "singular=1e-9", and then the problem with the "residual" went away.

 

I now have output separated by treatment and by time.  

 

And my appologies for wasting your time/energy due to a simple typo, but I do appreciate the enormous help figuring out the Slice code, I don't even want to know how long it would have taken me to figure that out myself. As a biologist, I wish I'd taken more statistics courses and gotten more experience with the software, rather than mucking around with pigs all day!

 

 

Data DSS;
Title 'DSS DSS-hLZ RBC';
Input pig $ Trt $ Time Resp;
Cards;
28-5	TRT1	0	7.08
28-8	TRT1	0	7.34
29-12	TRT1	0	8.09
29-14	TRT1	0	7.98
30-4a	TRT1	0	7.34
30-4b	TRT1	0	8.07
30-7	TRT1	0	8.41
28-6	TRT2	0	7
29-10	TRT2	0	6.9
30-6	TRT2	0	8.8
31-1	TRT2	0	7.08
28-1	TRT3	0	7.61
28-7	TRT3	0	7.57
30-5	TRT3	0	7.23
31-6	TRT3	0	6.6
28-3	CTRL	0	7.16
28-9	CTRL	0	7.54
30-3	CTRL	0	6.66
31-4	CTRL	0	8.28
31-5	CTRL	0	6.68
28-5	TRT1	1	7.11
28-8	TRT1	1	8.16
29-12	TRT1	1	8.25
29-14	TRT1	1	7.47
30-4a	TRT1	1	7.86
30-4b	TRT1	1	8.06
30-7	TRT1	1	7.89
28-6	TRT2	1	7.17
29-10	TRT2	1	7.54
30-6	TRT2	1	9.75
31-1	TRT2	1	7.66
28-1	TRT3	1	7.6
28-7	TRT3	1	7.26
30-5	TRT3	1	7.45
31-6	TRT3	1	6.4
28-3	CTRL	1	7.37
28-9	CTRL	1	7.84
30-3	CTRL	1	7
31-4	CTRL	1	8.08
31-5	CTRL	1	7.18
28-5	TRT1	2	7.09
28-8	TRT1	2	8.19
29-12	TRT1	2	8.28
29-14	TRT1	2	7.25
30-4a	TRT1	2	6.67
30-4b	TRT1	2	6.7
30-7	TRT1	2	7.73
28-6	TRT2	2	7.49
29-10	TRT2	2	7.05
30-6	TRT2	2	7.89
31-1	TRT2	2	6.47
28-1	TRT3	2	7.27
28-7	TRT3	2	7
30-5	TRT3	2	7.61
31-6	TRT3	2	6.38
28-3	CTRL	2	6.81
28-9	CTRL	2	7.54
30-3	CTRL	2	6.32
31-4	CTRL	2	7.06
31-5	CTRL	2	6.84
;
Proc GLIMMIX Data=DSS;
Class pig trt time;
Model Resp = trt time trt*time;
Random time / subject=pig type=ar(1) residual;
slice time*trt / sliceby=time sliceby=trt lines;
LSMEANS trt time/ diff lines;
Run; 

 

sas-innovate-2024.png

Join us for SAS Innovate April 16-19 at the Aria in Las Vegas. Bring the team and save big with our group pricing for a limited time only.

Pre-conference courses and tutorials are filling up fast and are always a sellout. Register today to reserve your seat.

 

Register now!

What is ANOVA?

ANOVA, or Analysis Of Variance, is used to compare the averages or means of two or more populations to better understand how they differ. Watch this tutorial for more.

Find more tutorials on the SAS Users YouTube channel.

Discussion stats
  • 5 replies
  • 2813 views
  • 5 likes
  • 2 in conversation