Programming the statistical procedures from SAS

Questions with mixed model

Reply
New Contributor
Posts: 4

Questions with mixed model

Hi,

I am trying to replicate the results from the meta analysis in "Meta-analysis of summary survival curve data" by Arends et al, in Statis. Med. 2008; 27:4381-4396. The data are a simulated set in that article. The model suggested is an iterated multivariate random effects model. The within study covariances are calculated as below

ln(ln(si jk))=b 0+b 1treati jk+b 2 ln(yeari jk)+b0+b1treati jk+b2 ln(yeari jk)+ eijk

bi∼N(0, D) and eijkN(0,Vi jk)

where sijk is the survival estimates from at yearijk the K-M table. i=1,..,m studies, j=1,2 treatment, k=1,..,T years (not necessarily the same or equally-spaced across the studies). Two studies are even missing one arm. Thus the model includes fixed and random effects for treatment and ln(year), apart from the intercept. V is a block diagonal matrix of known covariances with the blocks corresponding to each treatment study combination. For the iterated model, the initial V matrix (R matrix in SAS language) is calculated from the se of individual studies (in the K-M table) and the Sjik. In subsequent iterations, the estimated values for Sjk from the model are used. The method of calculating the covariance is given the paper. The D matrix is estimated by the model.

The complete data and the program are below. I deleted the observations with missing sijk and created a new variable NewTime which has gaps in it due to missing values. The year values are 0.5, 1, 1.5, 2, 2.5, 2.7 and 3. I created a partial V matrix.

data tmp1;

infile cards;

input StudyId Treat sijk1 sijk2 sijk3 sijk4 sijk5 sijk6 sijk7


/* sijk are the survival estimates and seijk are the corresponding se's */

datalines;

1    1 0.74 0.58 0.52 0.50 0.47  .    0.43 0.04 0.05 0.05 0.05 0.05  .    0.05

1    2 0.80 0.74 0.66 0.62 0.52  .    0.48 0.04 0.04 0.05 0.05 0.05  .    0.05

2    1  .    0.50  .    0.25  .      .    0.16  .    0.05  .    0.04    .      .    0.04

2    2  .    0.67  .    0.45  .      .    0.37  .    0.05  .    0.05    .      .    0.05

3    1 0.77  .    0.52  .    0.29  .      .    0.04  .    0.05  .      0.05  .      .

3    2 0.80  .    0.53  .    0.40  .      .    0.04  .    0.05  .      0.05  .      .

4    1  .    0.51  .      .      .      .    0.09  .    0.05  .      .        .      .    0.03

4    2  .    0.68  .      .      .      .    0.31  .    0.05  .      .        .      .    0.05

5    1  .    0.63  .    0.45  .      .      .      .    0.05  .    0.05    .      .      .

5    2  .      .      .      .      .      .      .      .      .      .      .        .      .      .

6    1  .    0.57  .      .      .      .    0.23  .    0.05  .      .        .      .    0.04

6    2  .    0.78  .    0.52  .      .    0.42  .    0.04  .    0.05    .      .    0.05

7    1  .      .      .      .      .      .      .      .      .      .      .        .      .      .

7    2 0.88 0.58 0.43 0.28 0.19  .      .    0.04 0.05 0.05 0.05  0.04  .      .

8    1  .    0.69  .      .      .      .    0.23  .    0.05  .      .      .      .    0.04

8    2  .    0.78  .      .      .      .    0.47  .    0.04  .      .      .      .    0.05

9    1  .    0.68  .      .      .      .      .      .    0.05  .      .      .      .      .

9    2  .      .      .      .      .      .    0.42  .      .      .      .      .      .    0.05

10  1  .      .      .      .      .    0.19  .      .      .      .      .      .    0.04  .

10  2  .      .      .      .      .    0.47  .      .      .      .      .      .    0.05  .

;

run;

data tmp2;

set tmp1;

array sijkt{*} sijk1-sijk7;

array sesijkt{*} sesijk1-sesijk7;

do i = 1 to 7;

year = 1/2;

if i = 6 then

year = 2.7;

if i = 7 then

year = 3;

Sijk = sijkt{i};

SeSijk = sesijkt{i};

if Sijk not = . then

do;

LnNegLnSijk = log(-log(Sijk));

LnYear = log(year);

end;

else

do;

LnNegLnSijk =.;

LnYear = .;

end;

output;

end;

drop i sijk1-sijk7 sesijk1-sesijk7;

run;

data tmp3;

set tmp2;

if SeSijk not = . then

    est = SeSijk*SeSijk;

if Sijk = . then

delete;

run;

data estvar;

    infile cards;

    input est;

cards;

0.05
;

run;

data vars;

    set tmp3;

    keep est;

run;

proc append base=estvar data=vars;

run;

proc sql noprint;

    select count(*)

          into :nobs

          from estvar;

quit;

Data sets tmp3 and estvar are attached;

I have two issues. How do I  code the parmsdata to provide the complete V matrix? In other words what should the structure of the V matrix look in an expanded form?

The second issue is I am not sure if I am coding the random and repeated statements correctly. I have not specified any covariance structure. If I run this code, I get the errorthat the Parmsdata should have 3 more or 10 less values.

proc mixed cl method=ml data=tmp3 ;

parms / parmsdata=estvar eqcons=2 to &nobs noiter;

/* the parms statement provides proc mixed with list of known within-study variances,

along with a starting value of 0.05 for the random between-studies variance component,

which is an unknown variance component to be estimated. */

     class StudyId Treat NewTime;

     model LnNegLnSijk = Treat LnYear / s CORRB COVB COVBI cl ddf = 1000, 1000;

     random int Treat lnyear / subject=StudyId ; /* Specifies that between-studies variance is a random

                                                                             effect to be estimated. */

     repeated NewTime / subject=Treat(StudyId) ; /*Specifies that each treatment arm study has its own unique variance. */

run;

Message was edited by: Srinivasan Rajagopalan I updated the post with the REPEATED and RANDOM statements that I think would be appropriate for the model described.

Attachment
Attachment
Ask a Question
Discussion stats
  • 0 replies
  • 260 views
  • 0 likes
  • 1 in conversation