setwd("D:/Tool APP/R/R_Workspace/RCS")
library(haven)
library(flextable)
library(finalfit)
library(rms)
library(ggplot2)
library(tibble)


#library(sparseM)
#install.packages('ggplot2')

data <- read_dta('D:/Tool APP/stata/RCS/17variables(drop).dta')

data<-as.data.frame(data)
#Setting factor and Numeric Variables
data$anxiety <- as.factor(data$anxiety)
data$sex <- as.factor(data$sex)
data$edu <- as.factor(data$edu)
data$FHistory <- as.factor(data$FHistory)
data$midsleeptime <- as.factor(data$midsleeptime)
data$complication <- as.factor(data$complication)
data$exercise <- as.factor(data$exercise)
data$bingcheng <- as.factor(data$bingcheng)
data$sleeptime1 <- as.factor(data$sleeptime1)
data$worksth <- as.numeric(data$worksth)
data$BMI <- as.numeric(data$BMI)
data$age <- as.numeric(data$age)
data$sleeptime <- as.numeric(data$sleeptime)
data$startsleep <- as.factor(data$startsleep)
data$startsleep1 <- as.factor(data$startsleep1)
data$startsleep3 <- as.factor(data$startsleep3)
data$startsleep4 <- as.factor(data$startsleep4)
data$startsleep5 <- as.factor(data$startsleep5)
data$sleeptime2 <- as.factor(data$sleeptime2)
data$sleeptime3 <- as.factor(data$sleeptime3)
data$lianhe <- as.factor(data$lianhe)
data$lianhe1 <- as.factor(data$lianhe1)
data$lianhe2 <- as.factor(data$lianhe2)
data$age1 <- as.factor(data$age1)
data$sleeptime4 <- as.factor(data$sleeptime4)
data$anxiety3 <- as.factor(data$anxiety3)


x <- c("sleeptime", "FHistory", "bingcheng", "age","sex","midsleeptime","complication","edu","exercise")
#x <- c("sleeptime3", "FHistory", "bingcheng", "age","midsleeptime","complication","edu","exercise")
y <- "anxiety"

data %>%
  finalfit.glm(dependent = y, explanatory= x) %>%
  flextable(cwidth = 1)

#make the result table
colon_s %>%
  finalfit(dependent = y, explanatory= x) %>%
  flextable(cwidth = 1) %>%
  save_as_docx(path = "table.docx")
#make the forest plot
data %>%
  or_plot(dependent = y, explanatory= x)

#make RCS curve
#1.find the optimal number of nodes
for (i in 3:7) {
  fit1<-glm(anxiety~sleeptime+BMI+FHistory+bingcheng+age+sex+midsleeptime+complication+exercise+edu,data=data,family=binomial())
  tmp <- extractAIC(fit1)
  if(i == 3) {AIC = tmp[2]; nk = 3}
  if(tmp[2] < AIC) {AIC = tmp[2]; nk = i} 
}
nk
#2.the optimal number og nodes above is 3

dd <- datadist(data)
options(datadist='dd') 
#3.fit a model
fit<- lrm(anxiety3~ rcs(sleeptime,3) +FHistory+bingcheng+age+sex+complication+exercise+edu,data=data) 
fit<- lrm(anxiety~ rcs(sleeptime,3) +FHistory+bingcheng+age+midsleeptime+complication+exercise+edu,data=data) 
anova(fit)#check for linerity，P<0.05
an<-anova(fit)

dd$limits$sleeptime[3] <-7###set the standard
fit1=update(fit)###updated models
OR<-Predict(fit1, sleeptime,fun=exp,ref.zero = TRUE) ##generating predictive values


ggplot()+
  geom_line(data=OR, aes(sleeptime,yhat),
            linetype="solid",size=1,alpha = 0.7,colour="#0070b9")+
  geom_ribbon(data=OR, 
              aes(sleeptime,ymin = lower, ymax = upper),
              alpha = 0.1,fill="#0070b9")+
  theme_classic()+
  geom_hline(yintercept=1, linetype=2,size=1)+
  #geom_vline(xintercept=48.89330,size=1,color = '#d40e8c')+ 
  labs( x="Night sleep duration", y="OR (95%CI)")

#RCS in different gender
OR1<-Predict(fit, sleeptime,sex=c('1','0'),fun=exp,ref.zero = TRUE) 
ggplot()+
  geom_line(data=OR1, aes(sleeptime,yhat, color = sex),
            linetype="solid",size=1,alpha = 0.7)+
  geom_ribbon(data=OR1, 
              aes(sleeptime,ymin = lower, ymax = upper,fill = sex),
              alpha = 0.1)+
  scale_color_manual(values = c('#0070b9','#d40e8c'))+
  scale_fill_manual(values = c('#0070b9','#d40e8c'))+
  theme_classic()+
  geom_hline(yintercept=1, linetype=2,size=1)+
  labs(title = "Anxiety Risk", x="sleep duration", y="OR (95%CI)") 

ggsave("sex2.tiff", units="in", width=5, height=4, dpi=300, compression = 'lzw')

#RCS in different age
OR2<-Predict(fit, sleeptime,age1=c('0','1','2','3','4'),fun=exp,ref.zero = TRUE,digits=5) 
ggplot()+
  geom_line(data=OR1, aes(sleeptime,yhat, color = age1),
            linetype="solid",size=1,alpha = 0.7)+
  geom_ribbon(data=OR1, 
              aes(sleeptime,ymin = lower, ymax = upper,fill = age1),
              alpha = 0.1)+
  scale_color_manual(values = c('#0070b9','#d40e8c','#F7BA0B'))+
  scale_fill_manual(values = c('#0070b9','#d40e8c','#F7BA0B'))+
  theme_classic()+
  geom_hline(yintercept=1, linetype=2,size=1)+
  labs(title = "Anxiety Risk", x="sleep duration", y="OR (95%CI)")