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: repost mixed effects fractional polynomial


From   Anders Alexandersson <andersalex@gmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: repost mixed effects fractional polynomial
Date   Wed, 8 Feb 2012 16:59:18 -0500

I came across the 3-week old posting below. William's code seems
modified from Royson's presentation "DIY fractional polynomials" at
the 2010 UK Stata Users Group meeting. DIY fractional polynomials is
an interesting topic. But the problem was badly worded as "this [no
code, no output] doesn't work - any suggestions?".
Will, do you still have the problem? I suggest you first use a set of
code for power of degree 1 (xtmixed with random effects) and then a
second set of code for power of degree 2 (xtmixed without random
effects), not the forvalues loop.

Anders Alexandersson
andersalex@gmail.com


On Jan 24, William Johnson <wojohnso@umn.edu> wrote:
> I am using the code at the end of this post to fit a two degree mixed
> effects fractional polynomial model.
>
> I can't get the code to allow ONLY the first degree to have mixed effects.
> My thought was that (names) should be replaced with (p1) on the random
> effects side of the xtmixed code, because the power combinations that
> fracpolypowers creates are called r(p1) and r(p2), but this doesn't work.
>
> Does anyone have any suggestions?
>
> Will Johnson
>
> // Store FP2 powers in local macros
> fracpoly_powers, degree(2)
> local np = r(np)
> forvalues j = 1 / `np' {
>                local p`j' `r(p`j')'
> }
> // Compute deviance for each FP model with covariate age gen x = age gen y =
> bmi local devmin 1e30 forvalues j = 1 / `np' {
>                qui fracgen x `p`j'', replace adjust(mean)
>                qui xtmixed y `r(names)' || id: `r(names)',
> covariance(unstructured) mle
>                local dev = -2 * e(ll)
>                if `dev' < `devmin' {
>                                local p `p`j''
>                                local devmin `dev'
>                }
>       di "powers = `p`j''" _col(20) " deviance = " %9.3f `dev'
> }
> di _n "Best model has powers `p', deviance = " `devmin'

*
*   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