Dear Daniel,
There are several papers which use gamlss for centile estimation.
See www.gamlss.org/ and click MORE ON GAMLSS, Books & Articles.
In particular Rigby and Stasinopoulos (2004, 2005, 2006, 2013)
and Stasinopoulos and Rigby (2007).
m3 <- gamlss(BMI~age, family=BCT, data=adatok_fiu)
fits a linear model in age for the median of BMI.
A better model is probably given by smoothing all 4 parameters of the BCT
distribution in age,
m4 <- gamlss(BMI~pb(age), sigma.fo=~pb(age), nu.fo=~pb(age), tau.fo=~pb(age),
family=BCT, data=adatok_fiu)
[It may also be better to transform age to x=log(age) or x=sqrt(age).]
pb(age) chooses automatically the effective degrees of freedom for
smoothing using a method of local maximum likelihood estimation (Rigby and
Stasinopoulos, 2013).
Alternatively a local method of Generalized Akaike Information Criterion
(GAIC) can be used to choose the effective degrees of freedom for smoothing
by
replacing pb(age) by e.g. pb(age, method=GAIC, k=3).
Use k = 2 for AIC
or k = log(n) for SBC.
There is no general agreement on what value of k to use.
We have often found k=3 gives sensible results.
Global methods of choose the effective degrees of freedom for smoothing are
available, see Rigby and Stasinopoulos (2004, 2005, 2006)
and Stasinopoulos and Rigby (2007).
Robert Rigby
Il 08/12/2013 16.45, Daniel Kehl ha scritto:
> Dear Community,
>
> I am struggling with a growth curve estimation problem. It is a classic
BMI change with age calculation. I am not an expert of this field but have
some statistical experience in other fields.
> Of course I started reading classical papers related to the topic and
understood the concept of the LMS method. After that I thought it will be a
"piece of cake", R must have a package related to the topic, so I just need
the data and I am just a few lines of code from the results.
>
> I encountered some problems:
> - found at least three packages related to LMS: gamlss, VGAM and an old
one lmsqreg (I did not try this because it does not support my R version
(3.0.1.))
> - it was relatively easy to get plots of percentiles in both packages,
although they were far not the same (I also tried an other software,
LMSchartmaker it gave different results from the previous ones)
> - I tried to get tables of predicted values (with the predict function in
VGAM and with centiles.pre in gamlss) but without any success.
> - I tried to use the function gamlss() instead of lms() in gamlss but I
could not force them to give the same (or very similar results), but the
centiles.pred() function did work as expected for the model resulted from
galmss()
> - lms gives really different results if k is specified different ways,
which is "best"?
>
> Also I have a general question: some publications state they estimated
the centiles so that aroun 18 years of age the curves pass through certain
points. How is that possible?
>
> Thank you for any suggestions or insights about the methods or preferred
package!
>
> Here is my code (without data):
>
> #####gamlss
> library(gamlss)
> library(VGAM)
> library(foreign)
> adatok <- read.spss("MDSZ adatok.sav", to.data.frame=TRUE)
>
> adatok_fiu <- subset(adatok, adatok$gender == "Fi??k")[,2:3]
> row.names(adatok_fiu) <- NULL
> adatok_lany <- subset(adatok, adatok$gender == "L??nyok")[,2:3]
> row.names(adatok_lany) <- NULL
>
> m1 <- lms(BMI,age,data=adatok_fiu, cent=c(3,10,25,50,75,90,97),
families="BCCG")
> fittedPlot(m1, x=adatok_fiu$age)
> m1 <- lms(BMI,age,data=adatok_fiu, cent=c(3,10,25,50,75,90,97),
families="BCCG", method.pb="GAIC", k=log(1455))
> fittedPlot(m1, x=adatok_fiu$age)
> m1 <- lms(BMI,age,data=adatok_fiu, cent=c(3,10,25,50,75,90,97),
families="BCCG", method.pb="GAIC")
> fittedPlot(m1, x=adatok_fiu$age)
>
> m2 <- lms(BMI,age,data=adatok_lany, cent=c(3,10,25,50,75,90,97),
families="BCCG")
> m2 <- lms(BMI,age,data=adatok_lany, cent=c(3,10,25,50,75,90,97),
families="BCCG", method.pb="GAIC", k=log(1144))
> m2 <- lms(BMI,age,data=adatok_lany, cent=c(3,10,25,50,75,90,97),
families="BCCG", method.pb="GAIC")
>
> m3 <- gamlss(BMI~age, family=BCT, data=adatok_fiu)
> centiles(m3,xvar=adatok_fiu$age, cent=c(3,10,25,50,75,90,97))
>
> newx <- seq(12,20,.5)
>
> centiles.pred(m1, xname="age", xvalues=newx)
> centiles.pred(m3, xname="age", xvalues=newx)
> centiles(m1,adatok_fiu$age)
>
> #####VGAM
> library(foreign)
> library(VGAM)
> library(VGAMdata)
>
> adatok <- read.spss("MDSZ adatok.sav", to.data.frame=TRUE)
>
> adatok_fiu <- subset(adatok, adatok$gender == "Fi??k")
> adatok_lany <- subset(adatok, adatok$gender == "L??nyok")
>
> fit1 <- vgam(BMI ~ s(age), lms.bcn(percentiles = c(3, 10, 25, 50, 75, 90,
97)), adatok_fiu)
> fit2 <- vgam(BMI ~ s(age), lms.bcn(percentiles = c(3, 10, 25, 50, 75, 90,
97)), adatok_lany)
>
> qtplot(fit1, percentiles = c(3, 10, 25, 50, 75, 90, 97), xlim = c(10.5,
20.5), ylim=c(13,34),
> las = 1, ylab = "BMI", lcol = 4, pch=NA)
>
> qtplot(fit2, percentiles = c(3, 10, 25, 50, 75, 90, 97), xlim = c(10.5,
20.5), ylim=c(13,34),
> las = 1, ylab = "BMI", lcol = 3, add=TRUE, pch=NA, label=FALSE)
>
>
> Thank you:
> Daniel
Companies Act 2006 : http://www.londonmet.ac.uk/companyinfo
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
