Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: st: Convert SAS code to STATA
From
"Joseph Coveney" <[email protected]>
To
<[email protected]>
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/