## library packages
library(segmented) #分段回归
library(splines) 
library(ggplot2)#绘图
library(Hmisc)
library(rms) #立方样条
library(nnet) #的multinom函数构建无序多分类logistic回归模型
library(openxlsx)
library(tidyverse)
## load data

load('diagnosis.RData')


## add labels

data.wide$Status <- factor(data.wide$Status, 
                           levels=c(0,1),
                           labels=c("No Disease", "Disease"))

data.wide$Sex <- factor(data.wide$Sex, 
                        levels=c(0,1),
                        labels=c("Male", "Female"))



# attach data.wide
attach(data.wide)


# 0. Fit model1 TestC~Age

model1 <- lm(TestC ~ Age)
y_model1 <- predict(model1)


# 1. Polynomial Regression

# fit polymomial regression

model1.poly2 <- lm(TestC ~ Age + I(Age^2))
model1.poly3 <- lm(TestC ~ Age + I(Age^2) + I(Age^3))

# check models
summary(model1.poly2)
summary(model1.poly3)

# compare models
anova(model1,model1.poly2)
anova(model1.poly2,model1.poly3)

# make predictions
y_model1.poly2 <- predict(model1.poly2)
y_model1.poly3 <- predict(model1.poly3)


# 2. Segmented regression


model1.segmented1 <- segmented(model1)
# model1.segmented2 <- segmented(model1,psi=c(35,50))

summary(model1.segmented1)
intercept(model1.segmented1)
slope(model1.segmented1)



plot(Age,TestC, pch=1, cex=1.5)
abline(a=coef(model1)[1],b=coef(model1)[2],col="red",lwd= 2.5)
plot(model1.segmented1, col='blue', lwd= 2.5 ,add=T)
# plot(model1.segmented2, col='green', lwd= 2.5 ,add=T)


# 3. Spline

# restricted cubic spline

model.spline1 <- lm(TestC ~ rcs(Age, 3))
model.spline2 <- lm(TestC ~ rcs(Age, 5))

# compare models (to determine number of knots)
AIC(model.spline1,model.spline2)

# make predictions
y_model.spline1 <- predict(model.spline1)
y_model.spline2 <- predict(model.spline2)




# or use ggplot

# 3 knots
ggplot(data=data.wide, aes(x=Age, y=TestC))  +
  geom_point()+geom_smooth(aes(x=Age, y=TestC), method="lm", formula=y ~ rcs(x, 3), se=FALSE, colour="blue")

# 5 knots
ggplot(data=data.wide, aes(x=Age, y=TestC))  +
  geom_point()+geom_smooth(aes(x=Age, y=TestC), method="lm", formula=y ~ rcs(x, 5), se=FALSE, colour="red")



# 4. Lowess

plot(Age,TestC, pch=1, cex=1.5)
lines(lowess(Age, TestC), col=2,lwd= 2.5)


# binary outcome
plot(Age, Status, pch=1, cex=1.5)
lines(lowess(Age, Status), col=2,lwd= 2.5)

analysedata <- sjydata %>% select(overweight,underweight,contains("Slope_"))
# 5.二分类logisistic回归的限制性立方样条
?lrm
sjydata <- read.xlsx("C:/Users/99405/Desktop/增长轨迹/xlsx数据/11.9sjydata6.xlsx")
analysedata$overweight <- as.factor(analysedata$overweight)
d3 <- datadist(analysedata) #打包数据集
options(datadist='d3')

fit1<-lrm(overweight~rcs(Slope_WEI_120,3),data=analysedata)
rcs1 <- Predict(fit1,Slope_WEI_120)
fit2<-lrm(overweight~rcs(Slope_WEI_240,3),data=analysedata)
rcs2 <- Predict(fit2,Slope_WEI_240)
fit3<-lrm(overweight~rcs(Slope_WEI_360,3),data=analysedata)
rcs3 <- Predict(fit3,Slope_WEI_360)
fit4<-lrm(overweight~rcs(Slope_WEI_480,3),data=analysedata)
rcs4 <- Predict(fit4,Slope_WEI_480)
fit5<-lrm(overweight~rcs(Slope_WEI_720,3),data=analysedata)
rcs5 <- Predict(fit5,Slope_WEI_720)
fit6<-lrm(overweight~rcs(Slope_WEI_960,3),data=analysedata)
rcs6 <- Predict(fit6,Slope_WEI_960)

p1 <- ggplot()+geom_line(data=rcs1, aes(Slope_WEI_120,yhat),linetype=1,linewidth=1,alpha = 0.8,colour="red")
p2 <- ggplot()+geom_line(data=rcs2, aes(Slope_WEI_240,yhat),linetype=1,linewidth=1,alpha = 0.8,colour="red")
p3 <- ggplot()+geom_line(data=rcs3, aes(Slope_WEI_360,yhat),linetype=1,linewidth=1,alpha = 0.8,colour="red")
p4 <- ggplot()+geom_line(data=rcs4, aes(Slope_WEI_480,yhat),linetype=1,linewidth=1,alpha = 0.8,colour="red")
p5 <- ggplot()+geom_line(data=rcs5, aes(Slope_WEI_720,yhat),linetype=1,linewidth=1,alpha = 0.8,colour="red")
p6 <- ggplot()+geom_line(data=rcs6, aes(Slope_WEI_960,yhat),linetype=1,linewidth=1,alpha = 0.8,colour="red")



# 无序多分类logisistic回归的限制性立方样条

d3 <- datadist(xob) #打包数据集
options(datadist='d3') 
fit1<-ols(你的Y~rcs(你样条化的X,c(-1.54887,-0.23625,0.93533))#c(-1.54887,-0.23625,0.93533)是你定的节点，修改
          +treatment协变量没有就删了,data=xob)
summary(fit1)
an<-anova(fit1)
p <-round(an[,5],3)
OLS1<-Predict(fit1,你样条化的X)
ggplot()+geom_line(data=OLS1, aes(你样条化的X,yhat),linetype=1,linewidth=1,alpha = 0.8,colour="red")+
  geom_ribbon(data=OLS1, aes(你样条化的X,ymin = lower, ymax = upper),alpha = 0.7,fill="grey")+
  theme_classic()+labs(title = paste0("RCS-你的Y-3"), x="你样条化的X", y="你的Y")+
  geom_text(family = "Times New Roman")+
  annotate("text", label =paste0("P-overall = ",ifelse(round(p[1],3) < 0.001,"< 0.001",round(p[1],3)),
                                 "\nP-non-linear = ",ifelse(round(p[2],3) < 0.001,"< 0.001",round(p[2],3))),
           x = 1, y = 65, size = 3, colour = "black")
p<-as.data.frame(p)
table1<-p
ggsave(file=paste0("F:/Rwork/RCSPlots/","你的Y","~","你样条化的X",3,".jpg"),width=10,height=8)