if (FALSE) {
#
library(gamlss)
boys <- boys7482
# calculate initial M curve
data <- na.omit(boys[,1:2])
f0154 <- gamlss(hgt~cs(age,df=15,c.spar=c(-1.5,2.5)),
sigma.formula=~cs(age,df=4,c.spar=c(-1.5,2.5)),
data=data,family=NO,
control=gamlss.control(n.cyc=3))
# calculate transformed age
t.age <- fitted(lm(data$age~fitted(f0154)))
t.age <- t.age - min(t.age)
data.t <- data.frame(data,t.age=t.age)
# calculate final solution
f0106r <- gamlss(hgt~cs(t.age,df=10,c.spar=c(-1.5,2.5)),
sigma.formula=~cs(t.age,df=6,c.spar=c(-1.5,2.5)),
data=data.t,family=NO,
control=gamlss.control(n.cyc=3))
# extract the LMS reference table in the 'classic' age grid
nl4.hgt.boys <- extractLMS(fit = f0106r, data=data.t, grid="compact",
dec = c(0,2,5))
nl4.hgt.boys
# flatten the reference beyond age 20Y (not very useful in this data)
nl4.hgt.boys.flat <- extractLMS(fit = f0106r, data=data.t, flatAge=20)
nl4.hgt.boys.flat
# use log age transformation
data.t <- data.frame(data, t.age = log(data$age))
f0106rlog <- gamlss(hgt~cs(t.age,df=10,c.spar=c(-1.5,2.5)),
sigma.formula=~cs(t.age,df=6,c.spar=c(-1.5,2.5)),
data=data.t,family=NO,
control=gamlss.control(n.cyc=1))
nl4.hgt.boys.log <- extractLMS(fit = f0106rlog, data=data.t)
nl4.hgt.boys.log
}
Run the code above in your browser using DataLab