library(openxlsx) #读写xlsx文件
library(tidyverse) #数据处理，绘图,批量化(map)
data <- read.xlsx("C:/Users/Administrator/Desktop/时间序列joinpoint.xlsx")
data3 <- read.xlsx("C:/Users/Administrator/Desktop/人口构成比.xlsx")
data1 <- data %>% filter(审核状态=="已终审卡")
data1$发病日期 <- as.Date(data1$发病日期,origin="1900-01-01")
data1 <- data1 %>% mutate(year=year(发病日期))
data1 <- data1 %>% mutate(occupation=case_when(人群分类=="农民"~1,
                                              TRUE~0))
data1 <- data1 %>% mutate(职业分类=case_when(人群分类=="家务及待业"~"其他",
                                         人群分类=="其他:居民"~"其他",
                                         人群分类=='其他:厨师'~"其他",
                                         TRUE~人群分类))

data1$职业分类 <- factor(data1$职业分类)
data1$职业分类 <- fct_infreq(data1$职业分类)
data1$人群分类 <- factor(data1$人群分类,level="农民","家务及待业",'其它')
data1$year <- factor(data1$year)
class(data1$year)
data1 <- data1 %>% 
data2 <- read.xlsx("C:/Users/Administrator/Desktop/澄城县人口数据.xlsx")
data2 <- data2[4:20,]
data2 <- data2 %>% rename("year"=年份)
data2 <- data2 %>% mutate(常住人口=as.numeric(常住人口))
result1 <- data1 %>% group_by(year) %>% summarise(n=n())
result2 <- data1 %>% group_by(性别,year) %>% summarise(n=n())
result3 <- data1 %>% group_by(occupation,year) %>% summarise(n=n())


result1 <- result1 %>% left_join(data2[,c(1,8)],by="year")
result1 <- result1 %>% mutate(incidence=n/常住人口/100)
result1 <- result1 %>% mutate(std=sqrt((1-incidence)*incidence/(常住人口*10000)))

result2 <- result2 %>% left_join(data2[,c(1,c(4,5))],by="year")
result2 <- result2 %>% replace_na(list(X4="151502",X5="152597"))
result2 <- result2 %>% mutate(X4=as.numeric(X4))
result2 <- result2 %>% mutate(X5=as.numeric(X5))
result2 <- result2 %>% mutate(incidence1=100*n/X4,incidence2=100*n/X5)
result2$incidence <- ifelse(result2$性别=="女",result2$incidence2,result2$incidence1)
result2 <- result2 %>% mutate(std1=sqrt((1-incidence)*incidence/X4))
result2 <- result2 %>% mutate(std2=sqrt((1-incidence)*incidence/X5))
result2$std <- ifelse(result2$性别=="女",result2$std2,result2$std1)

write.xlsx(result1,"C:/Users/Administrator/Desktop/joinpointresult.xlsx")
write.xlsx(result2,"C:/Users/Administrator/Desktop/joinpointresult2.xlsx")



ggplot(data1,aes(x=year,fill=factor(职业分类)))+geom_bar(position='fill')

plotdata <- data1 %>% group_by(year, 职业分类) %>% 
  summarise(count = n()) %>%
  mutate(percent = count*100/sum(count)) %>% mutate(label=paste0(sprintf("%.1f", percent), "%"))


ggplot(plotdata,aes(year,percent,fill=职业分类))+
  geom_bar(stat="identity",position = "fill")+
scale_y_continuous(labels = scales::percent_format(scale = 100))+ #百分比y轴
  labs(x="年份",y="人群占比",fill="")+
  scale_fill_brewer(palette="Set3")+
geom_text(aes(label=label),position=position_fill(vjust=0.5),size=2,color="black")

palette
RColorBrewer::display.brewer.all()