• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    公众号

HLM。RCS。限制性立方条。R语言实现RCS生成限制性立方条,画图,源码 ...

原作者: [db:作者] 来自: [db:来源] 收藏 邀请


#HLM

library(readxl)
library(xlsx)
mydata <- as.data.frame(read_excel("mydata_2019_12.xlsx"))
View(mydata)

#create new dataframe
mydata_new_constructed <- data.frame()
k=1
for (i in 1:nrow(mydata)){
  for (j in 2:35) {
    if (!is.na(mydata[i,j])) {
      mydata_new_constructed[k,1]<-mydata[i,1]
      mydata_new_constructed[k,2]<-as.numeric(colnames(mydata[j]))
      mydata_new_constructed[k,3]<-mydata[i,j]
      mydata_new_constructed[k,4]<-mydata[i,36]
      k<-k+1
    }
  }
}

colnames(mydata_new_constructed)<-c('id','age_week','weight_gain','BMI_group')

###############################################################generate the RCS variables based on age_week
RCS=as.data.frame(rcspline.eval(mydata_new_constructed$age_week, nk=6))
mydata_new_RCS = cbind(mydata_new_constructed,RCS)
write.table(mydata_new_RCS,'mydata_new_RCS.csv')

#############################################################begin to deal with the model
library(readxl)
mydata <- as.data.frame(read_excel("mydata_new_RCS.xlsx"))
names(mydata)<-c('id','age_week','weight_gain','BMI_group','age1','age2','age3','age4')

#delete the outliers original num = 9432, 
mydata<-mydata[mydata$weight_gain<32,]#丢去后9405
mydata<-mydata[mydata$weight_gain>-30,]#丢去后9423

mydata_new_constructed<-mydata_new_constructed[mydata_new_constructed$weight_gain<32,]
mydata_new_constructed<-mydata_new_constructed[mydata_new_constructed$weight_gain>-30,]


#if need to create the dependant variable
for (i in 1:nrow(mydata)){
if (mydata[i,4]==1) {
mydata[i,'LOGWTGAINKG']<-log(mydata[i,3]+4.5)
} else if (mydata[i,4]==2) {
mydata[i,'LOGWTGAINKG']<-log(mydata[i,3]+5.9)
} else if (mydata[i,4]==3) {
mydata[i,'LOGWTGAINKG']<-log(mydata[i,3]+7)
}
}

#build up the model
library(lme4)
Model3.1 = lmer(weight_gain ~ age1+age2+age3+age4+ (1|id) +(1+age1|BMI_group)+(1+age2|BMI_group)+(1+age3|BMI_group)+(1+age4|BMI_group), data=mydata,REML=FALSE)


summary(Model3.1)
coef(Model3.1)
Model3.1
anova(Model3.1)
VarCorr(Model3.1)
diag(VarCorr(Model3.1)$index)
attr(VarCorr(Model3.1),"sc")

#var example: 0.207 + 0.0002(GA Spline term1^2)+ 2(GA Spline term1)(-0.005) + 0.013

RCS_data <- unique(mydata[c(2,5,6,7,8)])
for (i in 1:nrow(RCS_data)){
  RCS_data[i,'mean1'] <- 2.867073+ 3.597614*RCS_data[i,2] -12.54554*RCS_data[i,3] +19.79477*RCS_data[i,4] -22.73170*RCS_data[i,5]
  RCS_data[i,'mean2'] <- 2.867073+ 3.568180*RCS_data[i,2]-12.49818*RCS_data[i,3] +19.79483*RCS_data[i,4]-22.73137*RCS_data[i,5]
  RCS_data[i,'mean3'] <- 2.867073+ 3.357253*RCS_data[i,2]-12.09662*RCS_data[i,3] +19.79483*RCS_data[i,4]-22.73169*RCS_data[i,5]
  RCS_data[i,'sd1'] <- sqrt((14.71+1.578) + 0.04079*(RCS_data[i,2]^2)+4.255)
  RCS_data[i,'sd2'] <- sqrt((14.71+1.578+0.0001188) + (0.04079+0.0003057)*(RCS_data[i,2]^2)+4.255)
  RCS_data[i,'sd3'] <- sqrt((14.71+1.578+0.0001188+0.06352) + (0.04079+0.0003057+0.004636)*(RCS_data[i,2]^2)+4.255)
}

write.table(RCS_data,'mean_sd_2019_12.csv')

#BMI 1
plot(mydata$age_week[mydata$BMI_group==1], mydata$weight_gain[mydata$BMI_group==1], col = "red")
points(RCS_data$age_week, RCS_data$mean1, col = "blue", cex = 1.5)

#group1 
group1<-mydata[mydata$BMI_group==1,]

#quantiles
table1<-data.frame(matrix(NA,ncol = 5))
for (gtage in sort(unique(group1$age_week))){
  temp<-quantile(group1$weight_gain[group1$age_week==gtage],probs = c(3,10,50,90,97)/100)
  table1<-rbind(table1,temp)
}
#sd
result_SD<-matrix(ncol = 2)
result_sd <- c()
for (gtage in sort(unique(group1$age_week))){
    temp_sum<-0
    for (observation in group1$weight_gain[group1$age_week==gtage]){
    temp_sum<-temp_sum+(observation - RCS_data$mean1[RCS_data$age_week==gtage])^2
    }
    temp_sd<-sqrt(temp_sum/length(group1$weight_gain[group1$age_week==gtage]))
    result_sd<-c(result_sd,temp_sd)
}

#result1
result1<-as.data.frame(cbind(sort(unique(group1$age_week)),result_sd,NA))
i = 1
for (gtage in sort(unique(group1$age_week))){
  result1[i,3]<-RCS_data$mean1[RCS_data$age_week==gtage]
  i<-i+1
}
names(result1)<-c('age','sd','mean')

#plot1
library(ggplot2)
ggplot(data =RCS_data ,aes(x=age_week,y=mean1))+ geom_smooth(se=F,
                                                      method="loess")+
  geom_point(mapping=aes(x=age_week, y=weight_gain), data=mydata_new_constructed)+
  geom_smooth(mapping=aes(x=age_week, y=mean1+sd1), data=RCS_data,se=F,method="loess",linetype="longdash")+
  geom_smooth(mapping=aes(x=age_week, y=mean1-sd1), data=RCS_data,se=F,method="loess",linetype="longdash")+
  geom_smooth(mapping=aes(x=age_week, y=mean1+2*sd1), data=RCS_data,se=F,method="loess",linetype="dotted")+
  geom_smooth(mapping=aes(x=age_week, y=mean1-2*sd1), data=RCS_data,se=F,method="loess",linetype="dotted")

#group2
group2<-mydata[mydata$BMI_group==2,]

#n
as.data.frame(table(group2$age_week))

#mean
RCS_data[order(RCS_data$age_week),c(1,7)]

#sd
RCS_data[order(RCS_data$age_week),c(1,10)]

#quantiles
quantile2<-matrix(NA,nrow = 1,ncol = 7)
result2<-as.data.frame(RCS_data[order(RCS_data$age_week),c(1,7,10)])
for (i in 1:nrow(result2)){
  temp <- qnorm(c(0.03,0.05,0.1,0.5,0.9,0.95,0.97),mean = result2$mean2[i],sd = result2$sd2[i])
  quantile2<-rbind(quantile2,temp)
}

#result2
result2$posi_sd<-result2$mean2+result2$sd2
result2$nega_sd<-result2$mean2-result2$sd2
result2$posi2_sd<-result2$mean2+2*result2$sd2
result2$nega2_sd<-result2$mean2-2*result2$sd2
result2[4:7]

#plot2
library(ggplot2)
ggplot(data =result2 ,aes(x=age_week,y=mean2))+ geom_smooth(se=F,
                                                      method="loess",col='black')+
  geom_point(mapping=aes(x=age_week, y=weight_gain), data=group2)+
  geom_smooth(mapping=aes(x=age_week, y=mean2+sd2), data=result2,se=F,method="loess",linetype="longdash",col='black')+
  geom_smooth(mapping=aes(x=age_week, y=mean2-sd2), data=result2,se=F,method="loess",linetype="longdash",col='black')+
  geom_smooth(mapping=aes(x=age_week, y=mean2+2*sd2), data=result2,se=F,method="loess",linetype="dotted",col='black')+
  geom_smooth(mapping=aes(x=age_week, y=mean2-2*sd2), data=result2,se=F,method="loess",linetype="dotted",col='black')+ 
  xlab("Gestational age (wk)") + 
  ylab('Gestational weight gain (kg)') +theme_bw() 

#covered points
tot_point<-c()
ok_point<-c()
for (age in result2$age_week){
  tot_point_temp<-length(group2$weight_gain[group2$age_week==age])
  tot_point<-c(tot_point,tot_point_temp)
  
  ok_point_temp<-sum(group2$weight_gain[group2$age_week==age]>result2$nega_sd[result2$age_week==age] & group2$weight_gain[group2$age_week==age]<result2$posi_sd[result2$age_week==age])
  ok_point<-c(ok_point,ok_point_temp)
}

coverage<-sum(ok_point)/sum(tot_point)
coverage

#group3
group3<-mydata[mydata$BMI_group==3,]

#n
as.data.frame(table(group3$age_week))

#quantiles
table3<-data.frame(matrix(NA,ncol = 5))
for (gtage in sort(unique(group3$age_week))){
  temp<-quantile(group3$weight_gain[group3$age_week==gtage],probs = c(3,10,50,90,97)/100)
  table3<-rbind(table3,temp)
}

#mean
RCS_data[order(RCS_data$age_week),c(1,8)]

#sd
result_sd3 <- c()
for (gtage in sort(unique(group3$age_week))){
  temp_sum<-0
  for (observation in group3$weight_gain[group3$age_week==gtage]){
    temp_sum<-temp_sum+(observation - RCS_data$mean3[RCS_data$age_week==gtage])^2
  }
  temp_sd<-sqrt(temp_sum/length(group3$weight_gain[group3$age_week==gtage]))
  result_sd3<-c(result_sd3,temp_sd)
}
as.data.frame(result_sd3)

#result3
result3<-as.data.frame(cbind(sort(unique(group3$age_week)),result_sd3,NA))
i = 1
for (gtage in sort(unique(group3$age_week))){
  result3[i,3]<-RCS_data$mean3[RCS_data$age_week==gtage]
  i<-i+1
}
names(result3)<-c('age','sd','mean')

#plot3
library(ggplot2)
ggplot(data =result3 ,aes(x=age,y=mean))+ geom_smooth(se=F,
                                                      method="loess")+
  geom_point(mapping=aes(x=age_week, y=weight_gain), data=group3)+
  geom_smooth(mapping=aes(x=age, y=mean+sd), data=result3,se=F,method="loess",linetype="longdash")+
  geom_smooth(mapping=aes(x=age, y=mean-sd), data=result3,se=F,method="loess",linetype="longdash")+
  geom_smooth(mapping=aes(x=age, y=mean+2*sd), data=result3,se=F,method="loess",linetype="dotted")+
  geom_smooth(mapping=aes(x=age, y=mean-2*sd), data=result3,se=F,method="loess",linetype="dotted")


#get the quantiles for group1
a=qnorm(p = 0.03,mean = result1$mean,sd = result1$sd)
b=qnorm(p = 0.1,mean = result1$mean,sd = result1$sd)
c=qnorm(p = 0.5,mean = result1$mean,sd = result1$sd)
d=qnorm(p = 0.9,mean = result1$mean,sd = result1$sd)
e=qnorm(p = 0.97,mean = result1$mean,sd = result1$sd)
cbind(a,b,c,d,e)

a=qnorm(p = 0.03,mean = result2$mean,sd = result2$sd)
b=qnorm(p = 0.1,mean = result2$mean,sd = result2$sd)
c=qnorm(p = 0.5,mean = result2$mean,sd = result2$sd)
d=qnorm(p = 0.9,mean = result2$mean,sd = result2$sd)
e=qnorm(p = 0.97,mean = result2$mean,sd = result2$sd)
cbind(a,b,c,d,e)

a=qnorm(p = 0.03,mean = result3$mean,sd = result3$sd)
b=qnorm(p = 0.1,mean = result3$mean,sd = result3$sd)
c=qnorm(p = 0.5,mean = result3$mean,sd = result3$sd)
d=qnorm(p = 0.9,mean = result3$mean,sd = result3$sd)
e=qnorm(p = 0.97,mean = result3$mean,sd = result3$sd)
cbind(a,b,c,d,e)


#############################################################3
##############################################################
######################for BMI = 2#############################
Model3.2 = lmer(weight_gain ~ age1+age2+age3+age4+age1:id+(1+age1|id)+(1+age2|id)+(1+age3|id)+(1+age4|id), data=group2,REML=FALSE)
summary(Model3.2)
coef(Model3.2)

qqnorm(mydata$weight_gain)
qqline(mydata$weight_gain)

qqnorm(log(group2$weight_gain))
qqline(log(group2$weight_gain))

library(fBasics)
ksnormTest(mydata$weight_gain)

ks.test(mydata$weight_gain,'pnorm')


########################################################################################

#labels
test <-mydata_new_constructed[1:200,]
test$Centiles<-factor(test$BMI_group,levels = c(1,2,3),labels = c("Mean", "1SD", "2SD"))


 s

plot+theme(legend.title = element_text(color="black", size=16, face="bold"))+theme(legend.text = element_text(color=c('black'), size = 14)) + scale_fill_manual(name = 'Centiles')
  scale_fill_discrete(name="Centiles",
                    breaks=c(1, 2, 3),
                    labels=c("50th", "10th and 90th", "3rd and 97th"))


#plot1
library(ggplot2)
ggplot(data =RCS_data ,aes(x=age_week,y=mean1))+ geom_smooth(se=F,
                                                             method="loess",col = 'red')+
  geom_point(mapping=aes(x=age_week, y=weight_gain), data=mydata_new_constructed[mydata_new_constructed$BMI_group=='1',])+
  geom_smooth(mapping=aes(x=age_week, y=mean1+sd1), data=RCS_data,se=F,method="loess",linetype="dotted",col='green')+
  geom_smooth(mapping=aes(x=age_week, y=mean1-sd1), data=RCS_data,se=F,method="loess",linetype="dotted",col='green')+
  geom_smooth(mapping=aes(x=age_week, y=mean1+2*sd1), data=RCS_data,se=F,method="loess",linetype="dashed",col= 'blue')+
  geom_smooth(mapping=aes(x=age_week, y=mean1-2*sd1), data=RCS_data,se=F,method="loess",linetype="dashed",col='blue') + 
  xlab('Gestational age (weeks)') + ylab('Gestational weight gain (kg)') + theme(axis.title.x = element_text(size = 16)) + theme(axis.title.y = element_text(size = 16))


#plot2
library(ggplot2)
ggplot(data =RCS_data ,aes(x=age_week,y=mean2))+ geom_smooth(se=F,
                                                             method="loess",col = 'red')+
  geom_point(mapping=aes(x=age_week, y=weight_gain), data=mydata_new_constructed[mydata_new_constructed$BMI_group=='2',])+
  geom_smooth(mapping=aes(x=age_week, y=mean2+sd2), data=RCS_data,se=F,method="loess",linetype="dotted",col='green')+
  geom_smooth(mapping=aes(x=age_week, y=mean2-sd2), data=RCS_data,se=F,method="loess",linetype="dotted",col='green')+
  geom_smooth(mapping=aes(x=age_week, y=mean2+2*sd2), data=RCS_data,se=F,method="loess",linetype="dashed",col= 'blue')+
  geom_smooth(mapping=aes(x=age_week, y=mean2-2*sd2), data=RCS_data,se=F,method="loess",linetype="dashed",col='blue')+ 
  xlab('Gestational age (weeks)') + ylab('Gestational weight gain (kg)')+ theme(axis.title.x = element_text(size = 16)) + theme(axis.title.y = element_text(size = 16))


#plot3
library(ggplot2)
ggplot(data =RCS_data ,aes(x=age_week,y=mean3))+ geom_smooth(se=F,
                                                             method="loess",col = 'red')+
  geom_point(mapping=aes(x=age_week, y=weight_gain), data=mydata_new_constructed[mydata_new_constructed$BMI_group=='3',])+
  geom_smooth(mapping=aes(x=age_week, y=mean3+sd3), data=RCS_data,se=F,method="loess",linetype="dotted",col='green')+
  geom_smooth(mapping=aes(x=age_week, y=mean3-sd3), data=RCS_data,se=F,method="loess",linetype="dotted",col='green')+
  geom_smooth(mapping=aes(x=age_week, y=mean3+2*sd3), data=RCS_data,se=F,method="loess",linetype="dashed",col= 'blue')+
  geom_smooth(mapping=aes(x=age_week, y=mean3-2*sd3), data=RCS_data,se=F,method="loess",linetype="dashed",col='blue')+ 
  xlab('Gestational age (weeks)') + ylab('Gestational weight gain (kg)') + theme(axis.title.x = element_text(size = 16)) + theme(axis.title.y = element_text(size = 16))


鲜花

握手

雷人

路过

鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
上一篇:
R语言操作数据库发布时间:2022-07-18
下一篇:
《数据科学:R语言实战》一导读发布时间:2022-07-18
热门推荐
热门话题
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap