1. analysis<-function(df,variable_list,type_list,group,group_list=0,t_test=FALSE)
    2. {
    3.   library(carData)
    4.   library(car)
    5.   library(nortest)
    7.   df <- df[which(!is.na(df[,group])),] #筛选掉没有分组的观测
    10.   if(group_list == 0)
    11.     group_list<-as.character(as.data.frame(table(df[,group]))[,1]) #group的类名
    13.   n_group<-length(unique(na.omit(df[,group]))) #group类数
    15.   #生成分析结果的dataframe
    16.   df_result<-data.frame(Variables=c('N',variable_list),Overall=NA)
    17.   each_group_num<-as.data.frame(table(df[,group]))
    18.   df_result[1,'Overall']<-count(df)
    20.   for(i in 1:n_group) #添加group的列
    21.   {
    22.     # print(group_list[i])
    23.     df_result[,group_list[i]]=NA
    24.     df_result[1,group_list[i]]=each_group_num[which(each_group_num[,1]==group_list[i]),2]
    25.   }
    27.   #添加总的四分位数列
    28.   df_result[,'Median']=NA
    29.   for(i in 1:n_group) #添加group的四分位数列
    30.   {
    31.     df_result[,paste(group_list[i],'Median')]=NA
    32.   }
    33.   df_result[,'是否正态']<-NA
    34.   df_result[,'是否方差齐性']<-NA
    35.   df_result[,'anova p-value']<-NA
    36.   df_result[,'kruskal-wallis p-value']<-NA
    37.   df_result[,'卡方检验 p-value']<-NA
    38.   if(t_test==TRUE)
    39.     df_result[,'t检验 p-value']<-NA
    40.   df_result[,'wilcox秩和检验 p-value']<-NA
    43.   #遍历字段
    44.   for(i in seq(1,dim(df_result)[1]-1))
    45.   {
    46.     mystr=variable_list[i] #字段名
    47.     df2<-df
    48.     df2<-df2[which(!is.na(df2[,group])),] #去group字段na的数据
    50.     #查找是否有对应的flag字段，若有，筛选出flag=0的数据
    51.     mystr2=ifelse(substring(mystr,1,3)=='val',paste('flag_',substring(mystr,5,nchar(mystr)),sep=""),paste('flag_',mystr,sep=""))
    52.     if(mystr2 %in% colnames(df))
    53.       df2<-df2[which(df2[,mystr2]==0),]
    55.     #连续和分类变量分别分析
    56.     if(type_list[i]=='连续')
    57.     {
    58.       df2<-df2[which(!is.na(df2[,mystr])),] #去自变量na的数据
    59.       #计算均值标准差，无法转成数值型的去掉
    60.       df_result[i+1,'Overall']=paste(as.character(formatC(mean(as.numeric(df2[,mystr]),na.rm = T),digits = 2, format = "f")),' ± ',as.character(formatC(sd(as.numeric(df2[,mystr]),na.rm = T),digits = 2, format = "f")), sep = "")
    61.       for(j in 1:n_group)
    62.       {
    63.         df_result[i+1,group_list[j]]=paste(as.character(formatC(mean(as.numeric(df2[which(df2[,group]==group_list[j]),mystr]),na.rm = T),digits = 2, format = "f")),' ± ',as.character(formatC(sd(as.numeric(df2[which(df2[,group]==group_list[j]),mystr]),na.rm = T),digits = 2, format = "f")),sep = "")
    64.         df_result[i+1,'Median']=paste(as.character(formatC(median(as.numeric(df2[,mystr]),na.rm = T),digits = 2, format = "f")),' (',as.character(formatC(quantile(as.numeric(df2[,mystr]),na.rm = T)[2],digits = 2, format = "f")),', ',as.character(formatC(quantile(as.numeric(df2[,mystr]),na.rm = T)[4],digits = 2, format = "f")),') ',sep = "")
    65.         df_result[i+1,paste(group_list[j],'Median')]=paste(as.character(formatC(median(as.numeric(df2[which(df2[,group]==group_list[j]),mystr]),na.rm = T),digits = 2, format = "f")),' (',as.character(formatC(quantile(as.numeric(df2[which(df2[,group]==group_list[j]),mystr]),na.rm = T)[2],digits = 2, format = "f")),', ',as.character(formatC(quantile(as.numeric(df2[which(df2[,group]==group_list[j]),mystr]),na.rm = T)[4],digits = 2, format = "f")),') ',sep = "")
    66.       }
    68.       #是否正态，样本量>5000时不能用shapiro 检验
    69.       df_result[i+1,'是否正态']=ifelse(ad.test(as.numeric(df2[,mystr]))$p.value>0.05,'是','否')
    71.       #是否方差齐性
    72.       df_result[i+1,'是否方差齐性']=ifelse(leveneTest(as.numeric(df2[,mystr])~df2[,group])[1,3]>0.05,'是','否')
    74.       #anova p-value
    75.       anova_p <- summary(aov(as.numeric(df2[,mystr])~df2[,group]))[[1]]["Pr(>F)"][1,1]
    76.       df_result[i+1,'anova p-value']<-ifelse(anova_p<0.000001,formatC(anova_p,format='e',digits=2),formatC(anova_p,digits = 6, format = "f"))
    78.       #kruskal-wallis p-value
    79.       kw_p <- kruskal.test(as.numeric(df2[,mystr])~df2[,group])$p.value
    80.       df_result[i+1,'kruskal-wallis p-value']<-ifelse(kw_p<0.000001,formatC(kw_p,format='e',digits=2),formatC(kw_p,digits = 6, format = "f"))
    82.       #t-test p-value
    83.       if(t_test==TRUE)
    84.       {
    85.         variance_equal <- ifelse(df_result[i+1,'是否方差齐性']=='是',TRUE,FALSE)  #根据方差是否齐性选择t检验函数参数
    86.         t_p <- t.test(as.numeric(df2[,mystr])~df2[,group],var.equal=variance_equal)$p.value
    87.         df_result[i+1,'t检验 p-value']<-ifelse(t_p<0.000001,formatC(t_p,format='e',digits=2),formatC(t_p,digits = 6, format = "f"))
    88.       }
    89.     }else
    90.     {
    91.       #计算freq和占比
    92.       t<-as.data.frame(table(df2[,mystr]))
    93.       df_result[i+1,'Overall']=paste(as.character(t[which(t[,1]==1),'Freq']),' (',as.character(formatC(t[which(t[,1]==1),'Freq']/dim(df)[1]*100,digits = 2, format = "f")),'%)',sep = "")
    94.       t<-as.data.frame(table(df2[,c(mystr,group)]))
    95.       for(j in 1:n_group)
    96.       {
    97.         df_result[i+1,group_list[j]]=paste(as.character(t[which(t[,1]=='1' & t[,2]==group_list[j]),'Freq']),' (',
    98.                                            as.character(formatC(t[which(t[,1]=='1' & t[,2]==group_list[j]),'Freq']/each_group_num[which(each_group_num[,1]==group_list[j]),2]*100,digits = 2, format = "f")),
    99.                                            '%)',sep='')
    100.       }
    102.       #卡方检验 p-value
    103.       chisq_p <- chisq.test(as.matrix(table(df2[,c(mystr,group)])))$p.value
    104.       df_result[i+1,'卡方检验 p-value']<-ifelse(chisq_p<0.000001,formatC(chisq_p,format='e',digits=2),formatC(chisq_p,digits =6, format = "f"))
    105.     }
    106.   }
    107.   print("0.00e+00代表该数小于1e-323")
    108.   return(df_result)
    109. }
    115. dev_yinshi <- dbGetQuery(con,"select * from dev_yinshi")
    116. data_v1_new<-dbGetQuery(con,"select * from data_v1_new")
    117. dev_yinshi<-dev_yinshi[which(!(is.na(dev_yinshi$id_980_v1) &
    118.                                  is.na(dev_yinshi$id_982_v1) &
    119.                                  is.na(dev_yinshi$id_981_v1) &
    120.                                  is.na(dev_yinshi$id_983_v1) &
    121.                                  is.na(dev_yinshi$id_984_v1) )),]
    122. dev_yinshi$group1<-NA
    123. dev_yinshi$group1[dev_yinshi$id_980_v1==0]<-0
    124. dev_yinshi$group1[dev_yinshi$id_980_v1==1]<-1
    126. dev_yinshi$group2<-NA
    127. dev_yinshi$group2[dev_yinshi$id_982_v1==98202]<-1
    128. dev_yinshi$group2[dev_yinshi$id_982_v1==98201 |
    129.                     dev_yinshi$id_982_v1==98203 ]<-0
    131. dev_yinshi$group3<-NA
    132. dev_yinshi$group3[dev_yinshi$id_981_v1==98101 |
    133.                     dev_yinshi$id_981_v1==98102 ]<-0
    134. dev_yinshi$group3[dev_yinshi$id_981_v1==98103 |
    135.                     dev_yinshi$id_981_v1==98104 |
    136.                     dev_yinshi$id_981_v1==98105 ]<-1
    138. dev_yinshi$group4<-NA
    139. dev_yinshi$group4[dev_yinshi$id_983_v1==98301]<-1
    140. dev_yinshi$group4[dev_yinshi$id_983_v1==98302 |
    141.                     dev_yinshi$id_983_v1==98303 |
    142.                     dev_yinshi$id_983_v1==98304 ]<-0
    144. dev_yinshi$group5<-NA
    145. dev_yinshi$group5[dev_yinshi$id_984_v1==98401 |
    146.                     dev_yinshi$id_984_v1==98402 |
    147.                     dev_yinshi$id_984_v1==98403 ]<-1
    148. dev_yinshi$group5[dev_yinshi$id_984_v1==98404]<-0
    149. dev_yinshi$score<-ifelse( ifelse(is.na(dev_yinshi$group1),0,dev_yinshi$group1)==1,1,0)+
    150.   ifelse( ifelse(is.na(dev_yinshi$group2),0,dev_yinshi$group2)==1,1,0)+
    151.   ifelse( ifelse(is.na(dev_yinshi$group3),0,dev_yinshi$group3)==1,1,0)+
    152.   ifelse( ifelse(is.na(dev_yinshi$group4),0,dev_yinshi$group4)==1,1,0)+
    153.   ifelse( ifelse(is.na(dev_yinshi$group5),0,dev_yinshi$group5)==1,1,0)
    154. dev_yinshi$score[is.na(dev_yinshi$group1) & is.na(dev_yinshi$group2) & is.na(dev_yinshi$group3) &
    155.                    is.na(dev_yinshi$group4) & is.na(dev_yinshi$group5) ]<-NA
    156. dev_yinshi$score_2<-ifelse( ifelse(is.na(dev_yinshi$group1),0,dev_yinshi$group1)==1,1,0)+
    157.   ifelse( ifelse(is.na(dev_yinshi$group2),0,dev_yinshi$group2)==1,1,0)+
    158.   ifelse( ifelse(is.na(dev_yinshi$group3),0,dev_yinshi$group3)==1,1,0)+
    159.   ifelse( ifelse(is.na(dev_yinshi$group4),0,dev_yinshi$group4)==1,1,0)+
    160.   ifelse( ifelse(is.na(dev_yinshi$group5),0,dev_yinshi$group5)==1,1,0)
    161. dev_yinshi$score_2[is.na(dev_yinshi$group1) | is.na(dev_yinshi$group2) | is.na(dev_yinshi$group3) |
    162.                      is.na(dev_yinshi$group4) | is.na(dev_yinshi$group5) ]<-NA
    163. dev_yinshi<-dev_yinshi[which(!is.na(dev_yinshi$score_2)),]
    164. ###id_979
    165. data_1<-merge(dev_yinshi,data_v1_new,by='user_id')
    166. data_1$db_duration_year<-as.numeric(data_1$diabetic_duration)/12
    167. data_1$ideal_smoking<-NA
    168. data_1$ideal_smoking[which(data_1$smoke_info==3)]<-as.character(1)
    169. data_1$ideal_smoking[which(data_1$smoke_info==1 | data_1$smoke_info==2)]<-as.character(0)
    170. data_1$current_drinker<-NA
    171. data_1$current_drinker[which(data_1$drink_info==3)]<-as.character(1)
    172. data_1$current_drinker[which(data_1$drink_info==1 | data_1$drink_info==2)]<-as.character(0)
    174. data_1$avg_pwv <- NA
    175. data_1$val_929003[which(data_1$flag_929003!=0)]<-NA
    176. data_1$val_929004[which(data_1$flag_929004!=0)]<-NA
    177. for (i in 1:nrow(data_1)) {
    178.   if(is.na(data_1$val_929003[i]) & !is.na(data_1$val_929004[i])) {
    179.     data_1$avg_pwv[i]<-as.numeric(data_1$val_929004[i])
    180.   }
    181.   if(!is.na(data_1$val_929003[i]) & is.na(data_1$val_929004[i])) {
    182.     data_1$avg_pwv[i]<-as.numeric(data_1$val_929003[i])
    183.   }
    184.   if (!is.na(data_1$val_929003[i]) & !is.na(data_1$val_929004[i])) {
    185.     data_1$avg_pwv[i]<-(as.numeric(data_1$val_929004[i])+as.numeric(data_1$val_929003[i]))*0.5
    186.   }
    187. }
    189. data_1$ABI <- NA
    190. data_1$val_929001[which(data_1$flag_929001!=0)]<-NA
    191. data_1$val_929002[which(data_1$flag_929002!=0)]<-NA
    192. for (i in 1:nrow(data_1)) {
    193.   if(is.na(data_1$val_929001[i]) & !is.na(data_1$val_929002[i])) {
    194.     data_1$ABI[i]<-as.numeric(data_1$val_929002[i])
    195.   }
    196.   if(!is.na(data_1$val_929001[i]) & is.na(data_1$val_929002[i])) {
    197.     data_1$ABI[i]<-as.numeric(data_1$val_929001[i])
    198.   }
    199.   if (!is.na(data_1$val_929001[i]) & !is.na(data_1$val_929002[i])) {
    200.     data_1$ABI[i]<-(as.numeric(data_1$val_929002[i])+as.numeric(data_1$val_929001[i]))*0.5
    201.   }
    202. }
    204. data_1$IMT <- NA
    205. data_1$val_925001[which(data_1$flag_925001!=0)]<-NA
    206. data_1$val_925009[which(data_1$flag_925009!=0)]<-NA
    207. for (i in 1:nrow(data_1)) {
    208.   if(is.na(data_1$val_925001[i]) & !is.na(data_1$val_925009[i])) {
    209.     data_1$IMT[i]<-as.numeric(data_1$val_925009[i])
    210.   }
    211.   if(!is.na(data_1$val_925001[i]) & is.na(data_1$val_925009[i])) {
    212.     data_1$IMT[i]<-as.numeric(data_1$val_925001[i])
    213.   }
    214.   if (!is.na(data_1$val_925001[i]) & !is.na(data_1$val_925009[i])) {
    215.     data_1$IMT[i]<-(as.numeric(data_1$val_925009[i])+as.numeric(data_1$val_925001[i]))*0.5
    216.   }
    217. }
    219. data_1$DR<-NA
    220. data_1[which(data_1$manual_read_res==0),'DR']<-0
    221. data_1[which(data_1$manual_read_res>=1),'DR']<-1
    223. group<-'id_979'
    224. variable_list<-c('age','sex','education_level',
    225.                  'diabetes_history','id_477','id_478','db_duration_year',
    226.                  'ideal_smoking','current_drinker','bmi','weight','vat',
    227.                  'waistline','hips','whr','sbp','dbp',
    228.                  'val_875001','val_875011','homa2','val_875004','val_875014',
    229.                  'val_876','val_878040','val_878041','val_878042','val_901013',
    230.                  'avg_pwv','ABI','IMT','DR','group2',
    231.                  'group3','group5','group4')
    232. type_list<-c('连续','离散','离散',
    233.              '离散','离散','离散','连续',
    234.              '离散','离散','连续','连续','连续',
    235.              '连续','连续','连续','连续','连续',
    236.              '连续','连续','连续','连续','连续',
    237.              '连续','连续','连续','连续','连续',
    238.              '连续','连续','连续','离散','离散',
    239.              '离散','离散','离散')
    245. df_result<-analysis(data_1,variable_list,type_list,group,group_list=0,t_test=FALSE)
    246. df_result$p<-ifelse(!is.na(df_result$`anova p-value`),df_result$`anova p-value`,df_result$`卡方检验 p-value`)
    247. tem1<-c('Variables','Overall','1','2','3','4','p')
    248. df_result_2<-df_result[,tem1]