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 [dwinsem...@comcast.net] Küldve: 2013. december 9. 21:40 To: Dániel Kehl Cc: r-help@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@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 ______________________________________________ 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.
