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
[[alternative HTML version deleted]]
dear Daniel, yet another package performing growth modelling is quantregGrowth. It uses quantile regressions with B-splines and quadratic penalties to ensure flexible estimation with additional noncrossing and monotonicity (optional) constraints. The paper underlying the package is here: http://link.springer.com/article/10.1007/s10651-012-0232-1 best, vito Il 08/12/2013 16.45, D?niel 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 > > [[alternative HTML version deleted]] > > > > ______________________________________________ > R-help at 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. >-- =============================================Vito M.R. Muggeo Dip.to Sc Statist e Matem `Vianelli' Universit? di Palermo viale delle Scienze, edificio 13 90128 Palermo - ITALY tel: 091 23895240 fax: 091 485726 http://dssm.unipa.it/vmuggeo 28th IWSM International Workshop on Statistical Modelling July 8-12, 2013, Palermo http://iwsm2013.unipa.it
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 classicBMI 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 andunderstood 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 oldone 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 inVGAM and with centiles.pre in gamlss) but without any success.> - I tried to use the function gamlss() instead of lms() in gamlss but Icould 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 estimatedthe 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 preferredpackage!>> 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:> DanielCompanies Act 2006 : http://www.londonmet.ac.uk/companyinfo [[alternative HTML version deleted]]
On Dec 8, 2013, at 7:45 AM, D?niel Kehl wrote:> 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. >You might want to look at a more recent discussion: http://www.who.int/entity/childgrowth/standards/velocity/tr3chap_2.pdf (The WHO centre has published their R code.)> 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.Don't see any code or error messages.> - 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"?Won't that depend on the amount and distribution of the data?> > 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 > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at 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.David Winsemius Alameda, CA, USA
Dear Daniel,
The centiles.pred is currently set up for a gamlss object fitted using the
function gamlss, not function lms
(but we will aim to change this soon).
The lms function uses calibration=TRUE as default
(this automatically adjusts the centile curves so that the sample %'s below
each centile curve equals the value of cent,
see Matsumoto et al. (2012) Energy Policy, 48, 820-828).
The gamlss function does not automatically do this.
The warning about 'safe' predictions can be caused by the prediction
range
of x being wider than the data range of x
(see section 3.4 of the gamlss manual, Stasinopoulos et al. (2008)
'Instructions on how to use the gamlss package in R, 2nd edition',
or see Chambers and Hastie (1992) 'Statistical Models in S').
Robert Rigby
Date: Tue, 10 Dec 2013 18:27:45 +0000
From: D?niel Kehl <kehld@ktk.pte.hu>
To: David Winsemius <dwinsemius@comcast.net>
Cc: "r-help@r-project.org" <r-help@r-project.org>
Subject: Re: [R] growth curve estimation
Message-ID:
<33D76D77E9AC4B438DA38B348ED6890D0B90B0A3@EMAIL.ktkdom.pte.hu>
Content-Type: text/plain; charset="iso-8859-2"
Dear Vito, Robert and David,
thank you for your replies, although I made some step forward on my own,
your answers helped a lot and gave me more insight.
The only question I have left is why this code gives me an error (and that
is what David asked for):
m1 <- lms(BMI,age,data=adatok_fiu, cent=c(3,10,25,50,75,90,97),
families="BCCG")
centiles.pred(m1, xname="age", xvalues=seq(10,20,.5),
cent=c(3,10,25,50,75,90,97))
Error in X[onlydata, , drop = FALSE] :
(subscript) logical subscript too long
but
m3 <- gamlss(BMI~pb(age), sigma.formula=~pb(age), nu.formula = ~pb(age),
tau.formula = ~pb(age), family=BCT, data=adatok_fiu)
centiles.pred(m3, xname="age", xvalues=seq(10,20,.5),
cent=c(3,10,25,50,75,90,97))
works as expected (although I get a warning message
Warning message:
In predict.gamlss(obj, what = "tau", newdata = newx, type =
"response", :
There is a discrepancy between the original and the re-fit
used to achieve 'safe' predictions
I had the feeling centiles.pred should work for the result of an lms()
function just like in case of gamlss().
Thank you all again!
daniel
Companies Act 2006 : http://www.londonmet.ac.uk/companyinfo
[[alternative HTML version deleted]]