@kc wrote:
I wonder, Why do you want a piecewise linear regression? How many breakpoints do you need? Do you know where the breakpoint(s) for the pieces is(are) (for example, at 6 months), or do you need to estimate the location(s) of the breakpoint(s)?
Well, from prior clinical knowledge, there is no significant effect of treatment on summary score beyond the 6 month timepoint. Therefore, there is need for only 2 breakpoints, one each at 1 month and 6 month. I included only one breakpoint at 6 month in my code as an example.
I'm thinking of a breakpoint as a value at which the slope changes, i.e., the boundary between the segments. You don't have any data prior to 1 month, so you can't have a breakpoint there.
Multiple methods (paired t-tests, ANCOVA, mixed effects models) are employed in analyzing these data. Also, an underlying assumption of longitudinal growth curve models is that the missing data is missing at random.
It is certainly convenient to assume that data are missing completely at random. But convenience does not necessarily make it true. If data are not MCAR, then any statistical method will be subject to bias.
So, any help with the syntax in running a piecewise regression to fill the table in my original post would be great. Let me know if more data (in CSV format this time) would be helpful in carrying out this task!
Here is some code to consider; I provide no guarantees so you'll want to understand it thoroughly. It includes some graphics that might help you understand what the model is doing and confirm visually that it might be doing what you want.
/* Create variable for breakpoint */
data plr;
set plr;
time_6 = max(time, 0.5); /* Breakpoint at 0.5 */
run;
proc tabulate data=plr;
class time time_6;
table time, time_6;
run;
/* Fit random coefficients model */
proc glimmix data=plr;
class subject_id group;
model summary_score = group|time group|time_6 / solution ;
random intercept time time_6 / subject=subject_id type=un g gcorr; /* random intercepts, random slopes */
output out=plr_out2 pred(noblup)=predpa pred=pred;
run;
proc sort data=plr_out2 out=plr_out2_srt;
by time;
/* Plot fitted regression for each subject and population-averaged regression */
proc sgpanel data=plr_out2_srt;
panelby group;
series x=time y=pred / group=subject_id markers;
series x=time y=predpa / lineattrs=(thickness=2 color=black);
run;
/* Plot population-averaged regression by group in one figure for comparison */
proc sgplot data=plr_out2_srt;
series x=time y=predpa / group=group;
run;
/* Plot observed data for each subject with its regression */
proc sgpanel data=plr_out2_srt;
panelby subject_id / columns=3;
series x=time y=pred / markers lineattrs=(thickness=2 color=black);
series x=time y=summary_score / markers;
run;
I'll give you a hint about acquiring the mean comparisons that you want: use the LSMEANS statement with the AT option.
I hope this helps.
... View more
