Date:         Mon, 14 Feb 2011 15:15:05 -0800
Reply-To:     Dale McLerran <stringplayer_2@YAHOO.COM>
Sender:       "SAS(r) Discussion" <SAS-L@LISTSERV.UGA.EDU>
From:         Dale McLerran <stringplayer_2@YAHOO.COM>
Subject:      Re: Linear Mixed Model Simulation: Random Intercepts and Slopes
In-Reply-To:  <AANLkTi=iuyUOM=327jsZ1BjT8gPq1Rov437-hzm0NtvC@mail.gmail.com>
Content-Type: text/plain; charset=iso-8859-1

Ryan,

In order to generate data in which the random effects are correlated,
you can take advantage of the Cholesky decomposition of a symmetric
positive definite matrix.  The Cholesky decomposition has the
property that:

      _                          _     _                        _
     |  a11    0    0   ...    0  |   |  a11  a21  a31  ...  ak1 |
     |  a21  a22    0   ...    0  |   |    0  a22  a32  ...  ak2 |
 V = |  a31  a32  a33   ...    0  | X |    0    0  a33  ...  ak3 |
     |  ...  ...  ...   ...  ...  |   |  ...  ...  ...  ...  ... |
     |  ak1  ak2  ak3   ...  akk  |   |    0    0    0  ...  akk |
      -                          -     -                        -

for V a symmetrix p.d. matrix.  For a 2x2 covariance matrix, we have

       _                             _
  V = |   a11*a11           0         |
      |   a21*a11   a21*a21 + a22*a22 |
       -                             -

So, if you want the random intercept variance to be V11,  the
variance of the slope to be V22, and the covariance to be
rho*sqrt(V11*V22), then you can construct

  V11 = <value>;
  V22 = <value>;
  rho = <value>;
  V21 = rho*sqrt(V11*V22);
  a11 = sqrt(V11);
  a21 = V21/a11;
  a22 = sqrt(V22 - a21*a21);

If you have two independent r.v.'s x0j and x1j and you construct

  u0j = a11*x0j;
  u1j = a21*x0j + a22*x1j;

then Cov(u0j,u1j)=V.

I don't have time right now to edit code.  But you should be
able to apply the above to get what you need.

Dale

---------------------------------------
Dale McLerran
Fred Hutchinson Cancer Research Center
mailto: dmclerra@nFoHsCpRaCm.org
Ph:  (206) 667-2926
Fax: (206) 667-5977
---------------------------------------


--- On Sun, 2/13/11, R B <ryan.andrew.black@GMAIL.COM> wrote:

> From: R B <ryan.andrew.black@GMAIL.COM>
> Subject: Re: Linear Mixed Model Simulation: Random Intercepts and Slopes
> To: SAS-L@LISTSERV.UGA.EDU
> Date: Sunday, February 13, 2011, 5:09 PM
> A minor correction to notation: The
> random slope, "Gamma01",
> technically should be "Gamma10". The code below should now
> be
> consistent with standard notation. If anyone happens to
> find other
> mistakes or perhaps a way to allow the random effects to
> be
> correlated, I'd welcome the feedback. -Ryan
>
> /* Generate Data */
> data response;
>
>  seed = 987624; /*seed for random # generator*/
>
>  do subject = 1 to 1000; /*N=1000 subjects*/
>
>  Gamma00 = 0.50; /*fixed intercept*/
>  Gamma10 = 0.30; /*fixed slope*/
>
>  u0j = sqrt(0.3)*rannor(seed); /*random intercept error
> term*/
>  u1j = sqrt(0.3)*rannor(seed); /*random slope error term*/
>
>  B0J = Gamma00 + u0j; /*random intercept*/
>  B1J = Gamma10 + u1j; /*random slope*/
>
>   do time=1 to 5; /*5 time points*/
>
>   eij = rannor(seed); /*error term*/
>
>   y = B0J + B1J*time + eij; /*full equation*/
>
>   output;
>  end;
>  end;
> run;
>
> /* Fit Model */
> proc mixed data=response;
> model y = time / solution;
> random intercept time / subject = subject;
> run;
>
> On Sun, Feb 13, 2011 at 12:34 PM, R B <ryan.andrew.black@gmail.com>
> wrote:
> > Hi,
> >
> > Below I generate data for a linear mixed model with
> random intercept
> > and slope terms which are considered a function of
> fixed intercept,
> > fixed slope, and error terms. As I understand it, the
> random slope
> > term suggests that the linear effect of time varies
> across subjects. I
> > believe the correlation between the random intercept
> and slope terms
> > are assumed to be zero in this simulation. That is, we
> are assuming
> > that as subject-specific intercepts change, the
> subject-specific
> > slopes do not change in a consistent way. Often, in
> real-world data,
> > one does observe such a correlation. I'm not sure how
> to adjust the
> > simulation code to allow for a non-zero correlation
> between these
> > random terms.
> >
> > As always, any thoughts on the simulation and MIXED
> code would be
> > appreciated. I'm particularly interested in finding
> out if there are
> > any mistakes in the simulation code.
> >
> > Ryan
> > --
> >
> > /* Generate Data */
> > data response;
> >
> >  seed = 987624; /*seed for random # generator*/
> >
> >  do subject = 1 to 1000; /*N=1000 subjects*/
> >
> >  Gamma00 = 0.50; /*fixed intercept*/
> >  Gamma01 = 0.30; /*fixed slope*/
> >
> >  u0j = sqrt(0.3)*rannor(seed); /*random intercept
> error term*/
> >  u1j = sqrt(0.3)*rannor(seed); /*random slope
> error term*/
> >
> >  B0J = Gamma00 + u0j; /*random intercept*/
> >  B1J = Gamma01 + u1j; /*random slope*/
> >
> >   do time=1 to 5; /*5 time points*/
> >
> >   eij = rannor(seed); /*error term*/
> >
> >   y = B0J + B1J*time + eij; /*full
> equation*/
> >
> >   output;
> >  end;
> >  end;
> > run;
> >
> > /* Fit Model */
> > proc mixed data=response;
> > model y = time / solution;
> > random intercept time / subject = subject;
> > run;
> >
>