Re: st: repost mixed effects fractional polynomial

[Anders Alexandersson <andersalex@gmail.com>]

I meant Patrick Royston and the link to the presentation is
http://www.stata.com/meeting/uk10/UKSUG10.Royston.ppt

Anders

On Wed, Feb 8, 2012 at 4:59 PM, Anders Alexandersson
<andersalex@gmail.com> wrote:
> 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/