doit<-function(var,lab,outcome,type='C',...){ library(survey) library(scales) library(Cairo) library(openxlsx) outcome.f<-paste0("~",deparse(substitute(outcome))) %>% as.formula() var.f<-paste0("~",deparse(substitute(var))) %>% as.formula() y.str<-deparse(substitute(outcome)) x.str<-deparse(substitute(var)) var<-enquo(var) outcome<-enquo(outcome) t<-theme_set(theme_bw()) t<-theme_update( axis.title.x=element_text(hjust=0.5), #axis.line=element_line(colour="black",lineend="square"), panel.grid.minor=element_blank(), #panel.grid.major.y=element_line(colour="black",size=0.1), #panel.grid.major=element_blank(), axis.ticks=element_line(colour="black",size=0.25), #axis.ticks.x=element_blank(), #strip.text=element_blank(), legend.position = "bottom", legend.direction="horizontal", legend.box="vertical", #legend.title=element_blank(), legend.text=element_text(size=8), legend.spacing=unit(0,"cm"), legend.margin=margin(b=0,t=0,unit='cm'), #axis.text.x=element_text(angle=45, hjust=1, vjust=1,size=8), axis.text.y=element_text(size=8), plot.title=element_text(size=12,margin=margin(b=10),hjust=0), plot.subtitle=element_text(size=8,color="darkslategrey",margin=margin(b=25,l=-25)), plot.caption = element_text(size=8,margin=margin(t=10),color="darkslategrey",hjust=0) ) if (type=='C'){ temp<-foranalysis %>% ungroup() %>% mutate(id=ifelse(is.na(ori),row_number(),ori)) %>% arrange(!!var) %>% mutate( pctl=row_number()/sum(is.na(!!var)==FALSE) ,quintile=case_when( pctl<=.2~1 ,pctl<=.4~2 ,pctl<=.6~3 ,pctl<=.8~4 ,pctl<=1~5 ) ) %>% group_by(quintile) dens<-temp %>% summarise(n(),den=sum(wtcomb*population)) lemas.svy<-svydesign(ids=~id,weights=~wtcomb,strata=~strat,data=temp,nest=TRUE) est<-svyby(outcome.f,~quintile,design=lemas.svy,svytotal) %>% left_join(dens,by="quintile") %>% mutate(rate=!!outcome*1000000/den,rate_se=se*1000000/den,ul=rate+1.96*rate_se,ll=rate-1.96*rate_se) print(est) forplot1<-foranalysis %>% ungroup() %>% mutate( rate=wtcomb*!!outcome*n()*1000000/sum(wtcomb*population) ,drate=wtcomb*!!outcome*1000000/(wtcomb*population) ,wt=wtcomb*population ,type=str_glue("Fatal Encounters ({y.str}) per 1M") ) forplot2<-foranalysis %>% ungroup() %>% mutate( rate=wtcomb*viol*n()*1000/sum(wtcomb*population) ,drate=wtcomb*viol*1000/(wtcomb*population) ,wt=wtcomb*population ,type="Violent Crimes per 1K" ) forplot<-bind_rows(forplot1,forplot2) fr<-forplot %>% mutate(num=case_when( type==str_glue("Fatal Encounters ({y.str}) per 1M")~wtcomb*!!outcome*1000000 ,type=="Violent Crimes per 1K"~wtcomb*viol*1000 ) ) %>% group_by(type) %>% summarise(n(),fulcrum=mean(!!var,na.rm=TRUE),rate=sum(num)/sum(wtcomb*population)) print(fr) forplot %>% group_by(type) %>% arrange(desc(rate)) %>% slice_head(n=10) %>% select(type,Year,prelink1,rate,drate,wt) %>% print() p<-forplot %>% ggplot(.)+ stat_smooth(inherit.aes=FALSE,aes(x=!!var,y=rate),formula=y~x,method="lm",color="blue",linetype=2,se=TRUE)+ geom_point(data=fr,inherit.aes=FALSE,aes(x=fulcrum,y=rate),color="white",fill="red",shape=21,size=4,alpha=0.7)+ facet_wrap(~type,nrow=1)+ labs(x=paste0("\n",lab),y="Rate\n") plotulim<-layer_scales(p)$y$range$range[[2]] plotllim<-layer_scales(p)$y$range$range[[1]] range<-plotulim-plotllim breaks.raw<-plotllim:plotulim %>% pretty(10) p<-p+ scale_y_continuous(labels=comma,breaks=breaks.raw[which(breaks.raw>=plotllim)],expand=expansion(mult = c(0.03,0.05)))+ annotation_custom(grob=grid::rectGrob(),xmin=-Inf,xmax=Inf,ymin=-Inf,ymax=plotllim)+ geom_boxplot(data=forplot,aes(x=!!var,y=plotllim-.05*range),orientation="y",width=0.05*range,outlier.alpha=0.1) } else if (type=='D'){ temp<-foranalysis %>% ungroup() %>% mutate(id=ifelse(is.na(ori),row_number(),ori)) %>% group_by(!!var) dens<-temp %>% summarise(n(),den=sum(wtcomb*population)) extargs<-c(...) l<-length(extargs) forplot1<-foranalysis %>% group_by(!!var) %>% mutate( rate=wtcomb*!!outcome*n()*1000000/sum(wtcomb*population) ,type=str_glue("Fatal Encounters ({y.str}) per 1M") ) forplot2<-foranalysis %>% group_by(!!var) %>% mutate( rate=wtcomb*viol*n()*1000/sum(wtcomb*population) ,type="Violent Crimes per 1K" ) forplot<-bind_rows(forplot1,forplot2) if (l>0 & (l %% 2==0)){ levels<-extargs[1:(l/2)] labels<-extargs[(l/2+1):l] forplot<-forplot %>% mutate(!!var:=factor(!!var,levels=levels,labels=labels)) temp<-temp %>% mutate(!!var:=factor(!!var,levels=levels,labels=labels)) dens<-dens %>% mutate(!!var:=factor(!!var,levels=levels,labels=labels)) } lemas.svy<-svydesign(ids=~id,weights=~wtcomb,strata=~strat,data=temp,nest=TRUE) est<-svyby(outcome.f,var.f,design=lemas.svy,svytotal) %>% left_join(dens,by=str_glue("{x.str}")) %>% mutate(rate=!!outcome*1000000/den,rate_se=se*1000000/den,ul=rate+1.96*rate_se,ll=rate-1.96*rate_se) print(est) means<-forplot %>% ungroup() %>% group_by(type,!!var) %>% summarise(n(),mean=mean(rate)) print(means) if (l>6){ t<-theme_update( axis.text.x=element_text(angle=45, hjust=1, vjust=1,size=8) ) } counts<-forplot %>% group_by(!!var) %>% summarise(n=n()/2) %>% mutate(nf=paste0("n=",prettyNum(n,big.mark=",",scientific=FALSE))) p<-forplot %>% ggplot(.,aes(x=factor(!!var),y=rate))+ stat_summary(fun.data="mean_cl_boot",geom="crossbar",fill="#00334d",alpha=0.3,position="dodge")+ facet_wrap(~type,nrow=1)+ labs(x=paste0("\n",lab),y="Rate\n") plotulim<-ggplot_build(p)$layout$panel_scales_y[[1]]$range$range[[2]] plotllim<-ggplot_build(p)$layout$panel_scales_y[[1]]$range$range[[1]] range<-plotulim-plotllim breaks.raw<-plotllim:plotulim %>% pretty(10) p<-p+ scale_y_continuous(labels=comma,breaks=breaks.raw[which(breaks.raw>=plotllim)],expand=expansion(mult = c(0.05,0.05)))+ annotation_custom(grob=grid::rectGrob(),xmin=-Inf,xmax=Inf,ymin=-Inf,ymax=plotllim)+ geom_text(data=counts,aes(x=!!var,y=plotllim-.05*range,label=nf)) } print(p) CairoPNG(str_glue("//Rtpnfil02/0218273_dcra/Stats/Output/Plots/LEMAS_{y.str}_{x.str}.png"),width=7,height=5,units="in",res=500) print(p) dev.off() write.xlsx(est,str_glue("//Rtpnfil02/0218273_dcra/Stats/Output/Plots/LEMAS_{y.str}_{x.str}.xlsx")) }
vaporwave
andulkaphoto
discoholicmusic
jacobvanloon