library(rio)
library(dplyr)
library(compareGroups)
library(riskRegression)
library(cmprsk)
library(prodlim)
library(survminer)
library(survival)
library(export)
library(ggplot2)
library(ggsci)
library(scales)
library(tidyverse)
library(Hmisc)
library(tableone)
library(naniar)
library(lubridate)
library(simputation)
library(ggprism)
library(janitor)
library(mice)
library(lattice)
library(MASS)
library(nnet)
library(missRanger)
library(forestmodel)
df=import('df.xlsx')

setwd(paste0(PATH,'/result'))

for (i in c(2,4,8:10,14)+1) {df[,i]=factor(df[,i])}

df$age=as.numeric(df$age)

df$group=ifelse(df$gender=='male'&df$SMI<=28.8,1,ifelse(df$gender=='female'&df$SMI<=20.8,1,0))

#-------------------------------------------------------------------group
df$group<-factor(df$group,levels = c(0,1),labels = c('normal muscle mass group','low muscle mass group'))





var0=c('SMI','T12SMRA','age','gender')
var0

vary=c('CT imaging features before treatment','Number of cavities before treatment',	'Whole Lung Lesion Proportions Before Treatment',	'Percentage of cavity absorption')
vary


#----------------------------------------------------------1
Y=vary[1]
Y
df$Y=df[,Y]
FML<-as.formula(paste0("Y~",paste0(var0,collapse = "+")))
fit<-glm(FML,family =gaussian(),df)
summary(fit)




#----------------------------------------------------------2
Y=vary[2]
Y
df$Y=df[,Y]
FML<-as.formula(paste0("Y~",paste0(var0,collapse = "+")))
fit<-glm(FML,family =gaussian(),df)
summary(fit)





#----------------------------------------------------------3
Y=vary[3]
Y
df$Y=df[,Y]
FML<-as.formula(paste0("Y~",paste0(var0,collapse = "+")))
fit<-glm(FML,family =gaussian(),df)
summary(fit)


#----------------------------------------------------------4
Y=vary[4]
Y
df$Y=df[,Y]
FML<-as.formula(paste0("Y~",paste0(var0,collapse = "+")))
fit<-glm(FML,family =gaussian(),df)
summary(fit)





vary=c('prognosis')
vary

var0=c('Skeletal muscle index','T12SMRA','age','gender')
var0
#----------------------------------------------------------1
Y=vary[1]
Y
df$Y=df[,Y]
FML<-as.formula(paste0("Y~",paste0(var0,collapse = "+")))
fit<-glm(FML,family =binomial(),df)
summary(fit)





var0=c('group','T12SMRA','age','gender')
var0
#----------------------------------------------------------1
Y=vary[1]
Y
df$Y=df[,Y]
FML<-as.formula(paste0("Y~",paste0(var0,collapse = "+")))
fit<-glm(FML,family =binomial(),df)
summary(fit)






w=10;h=8
library(psych)

varx=c('Skeletal muscle index','T12SMRA','BMI')
varx



vary=c('Proportion of total lung lesion before treatment ',' total lung consolidation before treatment ',' total lung strip shadow before treatment ',
       'Total lung glazing before treatment ',' total pulmonary focus absorption ratio ', 'percentage of cavity absorption')
vary
#---------------------------------------------------------------------------------------
dt<-df
cor_soil <- corr.test(dt[,varx],dt[,vary], method = "spearman", adjust = "none")#Do correlation
cp <- as.data.frame(cor_soil$p)
cr <- as.data.frame(cor_soil$r)
cr <- round(cr,3)

cp$names <- rownames(cp)
longp <- melt(cp,idvar = "names",v.names = "abd",direction = "long")
colnames(longp)[c(3)] <- c('p')

cr$names <- rownames(cr)
longr <- melt(cr,idvar = "names",v.names = "abd",direction = "long")

allnew <- as.data.frame(cbind(longr,longp$p))
colnames(allnew) <- c('names','variable','r','p')

allnew[which(allnew$p<0.001),'sig'] <- '***'
allnew[which(allnew$p<0.01 & allnew$p>0.001),'sig'] <- '**'
allnew[which(allnew$p<0.05 & allnew$p>0.01),'sig'] <- '*'


rt<-allnew
rt$r2<-sprintf('%0.3f',rt$r)
rt$pvalue<-sprintf('%0.3f',rt$p)
rt$pvalue<-ifelse(rt$pvalue=='0.000','<0.001',rt$pvalue)
rt


ggplot(rt, aes(variable,names)) +
    geom_tile(aes(fill=r)) +
    geom_text(aes(label=r2), color="black", size=6) + 
    geom_text(aes(label=sig), color="black", size=6,vjust = 1.8) + 
    scale_fill_gradient2(low='blue', high='red',mid = 'white', limit=c(-1,1),name=paste0("Correlation")) + 
    labs(x=NULL,y=NULL) + 
    theme_bw(base_size = 16)+
    theme(axis.text.x = element_text(size=16,angle = 45,hjust = 1,vjust = 1,color = "black"),
          axis.text.y = element_text(size=16,color = "black"),
          axis.ticks.y = element_blank(),
          panel.background=element_blank())+
    ggtitle('')





vary=c('Proportion of total lung lesion before treatment ',' total lung consolidation before treatment ',' total lung strip shadow before treatment ',
       'Whole-lung glass grinding before treatment',' whole-lung lesion absorption ratio')
vary
#---------------------------------------------------------------------------------------
dt<-df
cor_soil <- corr.test(dt[,varx],dt[,vary], method = "spearman", adjust = "none")#做相关性
cp <- as.data.frame(cor_soil$p)
cr <- as.data.frame(cor_soil$r)
cr <- round(cr,3)

cp$names <- rownames(cp)
longp <- melt(cp,idvar = "names",v.names = "abd",direction = "long")
colnames(longp)[c(3)] <- c('p')

cr$names <- rownames(cr)
longr <- melt(cr,idvar = "names",v.names = "abd",direction = "long")

allnew <- as.data.frame(cbind(longr,longp$p))
colnames(allnew) <- c('names','variable','r','p')

allnew[which(allnew$p<0.001),'sig'] <- '***'
allnew[which(allnew$p<0.01 & allnew$p>0.001),'sig'] <- '**'
allnew[which(allnew$p<0.05 & allnew$p>0.01),'sig'] <- '*'


rt<-allnew
rt$r2<-sprintf('%0.3f',rt$r)
rt$pvalue<-sprintf('%0.3f',rt$p)
rt$pvalue<-ifelse(rt$pvalue=='0.000','<0.001',rt$pvalue)
rt


ggplot(rt, aes(variable,names)) +
    geom_tile(aes(fill=r)) +
    geom_text(aes(label=r2), color="black", size=6) + 
    geom_text(aes(label=sig), color="black", size=6,vjust = 1.8) + 
    scale_fill_gradient2(low='blue', high='red',mid = 'white', limit=c(-1,1),name=paste0("Correlation")) + 
    labs(x=NULL,y=NULL) + 
    theme_bw(base_size = 16)+
    theme(axis.text.x = element_text(size=16,angle = 45,hjust = 1,vjust = 1,color = "black"),
          axis.text.y = element_text(size=16,color = "black"),
          axis.ticks.y = element_blank(),
          panel.background=element_blank())+
    ggtitle('')








varx=c('Skeletal muscle index','T12SMRA','BMI')
varx



vary=c('Proportion of total lung lesions before treatment, proportion of total lung lesions absorbed, percentage of cavity absorbed')
vary
#---------------------------------------------------------------------------------------
dt<-df
cor_soil <- corr.test(dt[,varx],dt[,vary], method = "spearman", adjust = "none")#做相关性
cp <- as.data.frame(cor_soil$p)
cr <- as.data.frame(cor_soil$r)
cr <- round(cr,3)

cp$names <- rownames(cp)
longp <- melt(cp,idvar = "names",v.names = "abd",direction = "long")
colnames(longp)[c(3)] <- c('p')

cr$names <- rownames(cr)
longr <- melt(cr,idvar = "names",v.names = "abd",direction = "long")

allnew <- as.data.frame(cbind(longr,longp$p))
colnames(allnew) <- c('names','variable','r','p')

allnew[which(allnew$p<0.001),'sig'] <- '***'
allnew[which(allnew$p<0.01 & allnew$p>0.001),'sig'] <- '**'
allnew[which(allnew$p<0.05 & allnew$p>0.01),'sig'] <- '*'


rt<-allnew
rt$r2<-sprintf('%0.3f',rt$r)
rt$pvalue<-sprintf('%0.3f',rt$p)
rt$pvalue<-ifelse(rt$pvalue=='0.000','<0.001',rt$pvalue)
rt


ggplot(rt, aes(variable,names)) +
    geom_tile(aes(fill=r)) +
    geom_text(aes(label=r2), color="black", size=6) + 
    geom_text(aes(label=sig), color="black", size=6,vjust = 1.8) + 
    scale_fill_gradient2(low='blue', high='red',mid = 'white', limit=c(-1,1),name=paste0("Correlation")) + 
    labs(x=NULL,y=NULL) + 
    theme_bw(base_size = 16)+
    theme(axis.text.x = element_text(size=16,angle = 45,hjust = 1,vjust = 1,color = "black"),
          axis.text.y = element_text(size=16,color = "black"),
          axis.ticks.y = element_blank(),
          panel.background=element_blank())+
    ggtitle('')




vary=c('Proportion of total lung lesion before treatment, proportion of total lung lesion absorption')
vary
#---------------------------------------------------------------------------------------
dt<-df
cor_soil <- corr.test(dt[,varx],dt[,vary], method = "spearman", adjust = "none")#做相关性
cp <- as.data.frame(cor_soil$p)
cr <- as.data.frame(cor_soil$r)
cr <- round(cr,3)

cp$names <- rownames(cp)
longp <- melt(cp,idvar = "names",v.names = "abd",direction = "long")
colnames(longp)[c(3)] <- c('p')

cr$names <- rownames(cr)
longr <- melt(cr,idvar = "names",v.names = "abd",direction = "long")

allnew <- as.data.frame(cbind(longr,longp$p))
colnames(allnew) <- c('names','variable','r','p')

allnew[which(allnew$p<0.001),'sig'] <- '***'
allnew[which(allnew$p<0.01 & allnew$p>0.001),'sig'] <- '**'
allnew[which(allnew$p<0.05 & allnew$p>0.01),'sig'] <- '*'


rt<-allnew
rt$r2<-sprintf('%0.3f',rt$r)
rt$pvalue<-sprintf('%0.3f',rt$p)
rt$pvalue<-ifelse(rt$pvalue=='0.000','<0.001',rt$pvalue)
rt


ggplot(rt, aes(variable,names)) +
    geom_tile(aes(fill=r)) +
    geom_text(aes(label=r2), color="black", size=6) + 
    geom_text(aes(label=sig), color="black", size=6,vjust = 1.8) + 
    scale_fill_gradient2(low='blue', high='red',mid = 'white', limit=c(-1,1),name=paste0("Correlation")) + 
    labs(x=NULL,y=NULL) + 
    theme_bw(base_size = 16)+
    theme(axis.text.x = element_text(size=16,angle = 45,hjust = 1,vjust = 1,color = "black"),
          axis.text.y = element_text(size=16,color = "black"),
          axis.ticks.y = element_blank(),
          panel.background=element_blank())+
    ggtitle('')