[R] growth curve estimation

[Dániel Kehl]

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

________________________________________
Feladó: David Winsemius [dwinsemius at comcast.net]
Küldve: 2013. december 9. 21:40
To: Dániel Kehl
Cc: r-help at r-project.org
Tárgy: Re: [R] growth curve estimation

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