GAMLSS包——LMS生长曲线

参考来自http://psy-data.github.io/2015/theory_page6.html

具体实现是我同学的数据，数据就不贴了，主要放代码

library(dplyr) #好像没有使用到dplyr的函数，可以不载入
library(gamlss)
# 计算儿童的年龄
age=difftime(as.Date("2019-01-01"),gamdata$BIRTH,units = "days")
age=round(age/365) #取整了

# 看一下儿童年龄的范围
max(newdata$age)
min(newdata$age)
gamdata$age=age

# 将Hb1A与年龄合并成新数据集，原始集里有其他变量
dia=gamdata$tanghua
bing=cbind.data.frame(age,dia)

# 缺失值可视化，
library(VIM)
aggr(bing,prop=FALSE,numbers=TRUE)
View(bing)
newdata=na.omit(bing) #删除缺失行
newdata=newdata[newdata$dia!=999,] #删除缺失值999的行

# 模型测试
gammodel=gamlss(dia~cs(age,df=3),family = BCT,data = newdata) #对Hb1A（dia）进行了box-cox变换，这里没有对sigma，tau，nu进行设定，对mu采用了三次样条拟合cs(age,df=3)
plot(gammodel) #绘制模型最后训练的参数的残差图，QQ图等
summary(gammodel) #观测每个参数最终的形式以及他们的估计

#绘制分位数曲线
centiles(gammodel,newdata$age,cent = c(1,99)) #绘制分位数曲线，cent设置图中需要的分位数

# 输出各个年龄下的Hb1A的1,99分位数的具体数值
range(newdata$age)  #观测年龄的取值范围

newx=seq(3,17,2) #将按年龄划分为从3-17岁的间隔2岁作为新的自变量

rest=centiles.pred(gammodel,xname = "age",xvalues = newx,type = "centiles",cent = c(1,99))
rest #1，99分位数下，3到17岁的儿童Hb1A的值

分男童女童分别绘图

male=gamdata[gamdata$GENDER==1,] #提取了男童的数据
dim(male)
aggr(male)
male=na.omit(male)
aggr(male)
range(male$tanghua)
male=male[male$tanghua!=13.1,] #删除了异常值

# 建立模型
malem=gamlss(tanghua~cs(age,df=3),family=BCT,data=male)
plot(malem)
summary(malem)

# 分位数曲线，没有给出图例，因为没找到更改图例位置的参数，legend=TRUE有图例，ylim限定了y轴范围
centiles(malem,male$age,cent = c(1,3,5,10,25,50,75,90,95,97,99),legend = FALSE,ylab = "HbA1c（%）",xlab = "Age",main = "Boys",ylim=c(3,10))

# 分位数数值表格，年龄划分为3到7岁，间隔1岁，直接存到工作目录路径下
result_male=centiles.pred(malem,xname = "age",xvalues = seq(3,17,1),cent =c(1,3,5,10,25,50,75,90,95,97,99) )
write.csv(result_male,"result_male.csv")

# 后面是女孩的
female=gamdata[gamdata$GENDER==2,]
aggr(female)
female=na.omit(female)
aggr(female)
range(female$tanghua)
female=female[female$tanghua<10&female$tanghua>2,]


femalem=gamlss(tanghua~cs(age,df=3),family = BCT,data=female)
plot(femalem)

centiles(femalem,female$age,cent =  c(1,3,5,10,25,50,75,90,95,97,99),legend = FALSE,ylab = "HbA1c（%）",xlab = "Age",main = "Girls")
result_girls=centiles.pred(femalem,xname = "age",xvalues = seq(3,17,1),cent = c(1,3,5,10,25,50,75,90,95,97,99) )
write.csv(result_girls,"result_girls.csv")

昂，不想解释结果，大家随便看看就好，如有需要调参，sigma.fo设定标准差分布，nu.fo设定偏度分布，tau.fo设定峰度分布。

可以用的设定

gamlss(dia~cs(age,df=3),sigma.fo=log(age),nu.fo=~1,tau.fo=~1,family=BCT,data=newdata）

family还有其他很多很多，可以具体看上面链接
下图为参考文献的绘制的儿童生长曲线，我的代码没有对图的样式进行过多的修改，因为太懒了。。。。