2
$\begingroup$

Newbie with linear mixed modelling here. Need some advice/help on my problem detailed below.

Dataset below with assessments taken at different visits: baseline (visit=1), and then at 30 days, 60 days and 90 days post-baseline), and the following variables:

DISEASE - the outcome, is an ordinal variable with 4 levels (1=normal, 2=mild, 3=moderate, 4=severe); BSL_DISEASE - the baseline value of DISEASE; VISIT - the visit; BNP - brain natriuretic peptide lab value (continuous); SEVERITY - binary variable derived based on DISEASE, i.e., if DISEASE in (0,1) then SEVERITY=0 (not severe), else SEVERITY=1 (severe);

Want to investigate the correlation between BNP (as predictor) and the DISEASE (as outcome). I basically want to look at change in serum concentrations of the BNP as marker for change in DISEASE. Here is what I have done so far to explore the correlation between BNP and DISEASE from baseline (visit=1) to Day 90 (visit=4) by using repeated measures logistic regression, implemented via PROC GLIMMIX. I fit the following model:

data have;
input ID$ DISEASE$ AGEGRP$ VISIT$ BNP SEVERITY$ BSL_DISEASE$;
datalines;
a001 1 1 1 1997.02 0 1
a001 1 1 2 1275.52 0 1
a001 4 1 3 180.23 1 1
a001 2 1 4 735.91 0 1
a002 1 2 1 454.16 0 1
a002 1 2 3 1776.52 0 1
a002 3 2 4 73.15 1 1
a003 1 2 1 1700.26 0 1
a003 3 2 2 1621.32 1 1
a003 2 2 4 850.65 0 1
a004 2 3 1 1963.25 0 2
a004 2 3 2 544.87 0 2
a004 4 3 3 768.54 1 2
a004 2 3 4 780.16 0 2
a005 1 2 1 655.24 0 1
a005 2 2 4 722.14 0 1
a006 1 1 1 1472.06 0 1
a006 1 1 4 749.78 0 1
a007 2 1 1 848.88 0 2
a007 2 1 2 1482.78 0 2
a007 3 1 4 735.26 1 2
a008 1 1 1 1752.35 0 1
a008 1 1 2 1698.82 0 1
a008 3 1 3 1871.25 1 1
a008 4 1 4 587.35 1 1
a009 1 3 1 1549.89 0 1
a009 3 3 3 785.52 1 1
a009 1 3 4 384.72 0 1
a010 3 3 1 1211.95 1 3
a010 3 3 4 1596.38 1 3
a011 4 1 1 1785.45 1 4
a011 4 1 4 644.12 1 4
a012 3 3 1 798.28 1 3
a012 3 3 2 742.69 1 3
a012 3 3 3 1423.59 1 3
a012 3 3 4 1089.47 1 3
;
run;
proc glimmix data=have noclprint; 
class ID VISIT (ref='1'); 
model DISEASE (event='1')= BNP VISIT/ dist=mult link=clogit solution; 
random VISIT/subject=ID residual type=CS; 
random INT/subject=ID type=CS;
output out=FITDAT pred(ilink noblup)=predprob; 
NLOPTIONS tech=NRRIDG Maxiter=1000; 
run;

But I get an error message that "R side random effects are not supported for the multinomial" so I deleted the random VISIT statement and it converges now but my questions are:

  1. Is this model the correct one to fit to the data in order to address my objective?
  2. Don't I need a random VISIT statement? My understanding is that I need to impose some sort of covariance structure on visit, otherwise we're just assuming that the values at the various visits are not correlated which I'm not sure is accurate?
$\endgroup$

1 Answer 1

6
$\begingroup$

I'm not sure I can give an adequate answer to your first question, which is quite general. I'm not sure you can ever tell if a model is the correct one unless you also know the data generating mechanism exactly. Both of your models address the fact that observations within subjects are not independent. They do not have the same covariance structure however and the one including VISIT is closer to what I would expect to be called 'repeated measures', allowing for covariance across different visits.

To address the problem with the random statement(s): you can fit the exact same model through a G-side random effect by removing the / residual option. Presumably the multinomial model won't let you increase the dimension of R because the it already has multiple intercepts, one per response category. The problem then becomes that the default G-side denominator degrees of freedom estimation method does not leave any for the BNP effect, so I'd update that to a more sensible one by including / ddfm=kr (or a method of your choice) in the model statement.

There's also no need to have separate random statements, you can get the same parameters estimated through a single one:

random intercept visit / subject=id type=vc;

This might not matter much for small datasets, but it should be more efficient to do it like this. I'm quite sure SAS does not handle them exactly the same internally, because there are tiny numerical differences in the final V matrix.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.