Bookmark and Share

Notice: On March 31, it was announced that Statalist is moving from an email list to a forum. The old list will shut down on April 23, and its replacement, statalist.org is already up and running.


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: st: Convert SAS code to STATA


From   "Joseph Coveney" <jcoveney@bigplanet.com>
To   <statalist@hsphsun2.harvard.edu>
Subject   Re: st: Convert SAS code to STATA
Date   Sun, 23 Sep 2012 20:31:07 +0900

Marquis Hawkins wrote:

My goal was to get change per year.  with SAS using the contrast "-0.333 0.333 0
0 0 0 ", this compares the average change in GFR over three years between the
two groups.  I'm not sure how to do this in STATA.

In SAS I using the following:

PROC mixed data= habc.habclong covtest;
class habcid ckd1 time;
model gfr = ckd1*time/solution ;
repeated time/type=un subject=habcid; 
estimate 'avg change/year between baseline and year 3 for grp 1' ckd1*time
-0.333 0.333 0 0 0 0 ; 
estimate 'avg change/year between year 3 and year 10 for grp 1' ckd1*time 0
-0.143 0.143 0 0 0 ; 
estimate 'avg change/year between baseline and year 3 for grp 2' ckd1*time 0 0 0
-0.333 0.333 0 ; 
estimate 'avg change/year between year 3 and year 10 for grp 2' ckd1*time 0 0 0
0 -0.143 0.143 ; 
where inclusion=1; 
run;

When I use the following STATA code:

xi: xtmixed gfr totpa3cat##time if inclusion==1 ///
 || habcid: , noconstant residuals(uns, t(time)) nolog reml
 
margins totpa3cat time totpa3cat##time if inclusion==1, contrast

It does not give me the same information.  It only gives me the mean GFR at each
time point by CKD group.  I can't put in the contrast here the same as I do it
in SAS.  Any advice?

--------------------------------------------------------------------------------

Try the following.  The first part just sets up a fake dataset for illustration.
The model and contrasts start at "Begin here".  

You can use the fully parameterized model (main effects and interaction term)
like you show for Stata, but it might be a little simpler to compute the desired
annual rate of GFR change by excluding the main effects and include just the
interaction term alone as you do for your SAS code.

It helps (even for setting up the vector for SAS's ESTIMATE and CONTRAST
statements) to see the beta vector (vector of regression coefficients and
covariance parameters) for orientation.  I like the way that Stata labels the
coefficients--it makes life easier when you go to set up the contrast (or, in
this case, SAS ESTIMATE statement / Stata -lincom- postestimation command).

Joseph Coveney

. version 11.2

. 
. set more off

. local linesize `c(linesize)'

. set linesize 80

. 
. drawnorm gfr0 gfr3 gfr10, double mean(105 105 105) ///
>         sd(20 20 20) corr(1 0.7 0.4 \ 0.7 1 0.7 \ 0.4 0.7 1) ///
>         seed(`=date("2012-09-23", "YMD")') n(150) clear
(obs 150)

. generate byte ckd1 = 1

. tempfile tmpfil0

. quietly save `tmpfil0'

. 
. drawnorm gfr0 gfr3 gfr10, double mean(105 90 45) ///
>         sd(20 20 20) corr(1 0.7 0.4 \ 0.7 1 0.7 \ 0.4 0.7 1) ///
>         n(150) clear
(obs 150)

. generate byte ckd1 = 2

. append using `tmpfil0'

. 
. generate int habcid = _n

. 
. generate byte inclusion = 1

. quietly replace inclusion = 0 if runiform() < 0.05

. 
. quietly reshape long gfr, i(habcid) j(time)

. 
. *
. * Begin here
. *
. xtmixed gfr i.ckd1#i.time if inclusion == 1 || habcid: , noconstant reml ///
>         residuals(unstructured, t(time)) nolrtest nolog

Mixed-effects REML regression                   Number of obs      =       861
Group variable: habcid                          Number of groups   =       287

                                                Obs per group: min =         3
                                                               avg =       3.0
                                                               max =         3


                                                Wald chi2(5)       =   1480.43
Log restricted-likelihood = -3584.0104          Prob > chi2        =    0.0000

------------------------------------------------------------------------------
         gfr |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
   ckd1#time |
       1  3  |    .928403   1.258088     0.74   0.461    -1.537404     3.39421
       1 10  |   1.760579   1.844002     0.95   0.340    -1.853598    5.374756
       2  0  |  -.5753795   2.451356    -0.23   0.814     -5.37995    4.229191
       2  3  |  -16.28983   2.381158    -6.84   0.000    -20.95681   -11.62284
       2 10  |  -58.79057   2.389539   -24.60   0.000    -63.47398   -54.10716
             |
       _cons |   104.5397   1.730348    60.42   0.000     101.1483    107.9312
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
habcid:              (empty) |
-----------------------------+------------------------------------------------
Residual: Unstructured       |
                      sd(e0) |   20.76418   .8697152      19.12766    22.54071
                      sd(e3) |    19.5612   .8193256      18.01949     21.2348
                     sd(e10) |   19.70681   .8254263      18.15363    21.39287
                 corr(e0,e3) |   .7212095   .0284242      .6607045    .7724036
                corr(e0,e10) |    .403059   .0496116      .3015267    .4955525
                corr(e3,e10) |   .6976497   .0304043      .6331215     .752548
------------------------------------------------------------------------------

. 
. matrix list e(b)

e(b)[1,13]
             gfr:          gfr:          gfr:          gfr:          gfr:
         1b.ckd1#      1b.ckd1#      1b.ckd1#       2.ckd1#       2.ckd1#
         0b.time        3.time       10.time       0b.time        3.time
y1             0     .92840295     1.7605791    -.57537955    -16.289825

             gfr:          gfr:      lnsig_e:   r_atr1_1_2:   r_atr1_1_3:
          2.ckd1#                                                        
         10.time         _cons         _cons         _cons         _cons
y1    -58.790573     104.53974     3.0332294       .910161     .42729588

     r_lns1_2ose:   r_atr1_2_3:  r_lns1_3ose:
                                             
           _cons         _cons         _cons
y1    -.05968157     .86270682    -.05226528

. 
. display in smcl as text "avg change/year between baseline and year 3 for grp 1
> "
avg change/year between baseline and year 3 for grp 1

. lincom _b[1.ckd1#3.time] / 3

 ( 1)  .3333333*[gfr]1b.ckd1#3.time = 0

------------------------------------------------------------------------------
         gfr |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |   .3094677   .4193626     0.74   0.461     -.512468    1.131403
------------------------------------------------------------------------------

. 
. display in smcl as text "avg change/year between year 3 and year 10 for grp 1"
avg change/year between year 3 and year 10 for grp 1

. lincom (_b[1.ckd1#10.time] - _b[1.ckd1#3.time]) / 7

 ( 1)  - .1428571*[gfr]1b.ckd1#3.time + .1428571*[gfr]1b.ckd1#10.time = 0

------------------------------------------------------------------------------
         gfr |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |   .1188823   .1817676     0.65   0.513    -.2373757    .4751403
------------------------------------------------------------------------------

. 
. display in smcl as text "avg change/year between baseline and year 3 for grp 2
> "
avg change/year between baseline and year 3 for grp 2

. lincom _b[2.ckd1#3.time] / 3

 ( 1)  .3333333*[gfr]2.ckd1#3.time = 0

------------------------------------------------------------------------------
         gfr |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |  -5.429942   .7937192    -6.84   0.000    -6.985603   -3.874281
------------------------------------------------------------------------------

. 
. display in smcl as text "avg change/year between year 3 and year 10 for grp 2"
avg change/year between year 3 and year 10 for grp 2

. lincom (_b[2.ckd1#10.time] - _b[2.ckd1#3.time]) / 7

 ( 1)  - .1428571*[gfr]2.ckd1#3.time + .1428571*[gfr]2.ckd1#10.time = 0

------------------------------------------------------------------------------
         gfr |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |  -6.071535   .1824021   -33.29   0.000    -6.429037   -5.714034
------------------------------------------------------------------------------

. 
. set linesize `linesize'

. 
. exit

end of do-file


*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/


© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   Site index