library(glmnet) library(data.table) options("width"=420) #C2020dir2="/mnt/vafs01/data/VA_MVP011_CDW/mcmahonb/" #C2020dir="/mnt/vafs01/data/VA_MVP011_CDW/mcmahonb/C2022/" C2020dir2="/lustre/valustre/VA_MVP011_CDW/mcmahonb/" C2020dir="/lustre/valustre/VA_MVP011_CDW/mcmahonb/C2022/" Mode=function(x){ux=unique(x);ux[which.max(tabulate(match(x,ux)))]} ndiA = readRDS(paste(c(C2020dir2,"/ndi.RDS"),collapse="")) ndiA$OCDate=as.numeric(julian(as.Date(as.character(ndiA$DOD_NDI),format="%m/%d/%Y"),origin=as.Date("2000-01-01")))/365.25 ndiA=ndiA[MatchFound==1,] CntF=1.0 set.seed(100) #Run="encode" #study="Rcc" #study="Pds-e" study="Pds-c" ;CntF=0.25 #study="nested" Run="PullData" #cohort=50 print(Run) if (Run == "PullData"){ #for (cohort in 86){ #for (cohort in c(36,39,41,44,46,49,51,54,56,59,61,64,66,69,71,74,76,79,81,84)){ for (cohort in seq(0,100,1)){ #cohort=24 #for (cohort in seq(1,2,3)){ print(cohort) # for (cohort in c(1,2,4,67)){ # for (cohort in c(3)){ #define cohort and outcome #predict suicide attempt after oct 1, 2015, but keep track of date, and keep only patients seen at least 2 years prior to outcome. Select controls with randomly chosen Dx code date for outcome. ndi=ndiA visits=readRDS(paste(c(C2020dir,"outpat/oFac.",cohort,".RDS"),collapse="")) # for stopcodes, 323 is primary care / medicine, 502 is mental health abbout 10x as many primary care visits as mental health if (study=="Pds-e"){ OCDate=10.5;tWin=3;tGap=0. visitT=visits[(VisitDateTime < OCDate-tGap & VisitDateTime > OCDate-tWin-tGap & StopCode==502) , .(PatientICN,VisitDateTime)];colnames(visitT)[2]="OCDate" pts=visitT[,min(OCDate),by="PatientICN"];colnames(pts)[2]="Tpr" rm(visits) } if (study=="Pds-c"){ timePr=17.; OCDate=timePr;tWin=0.25;tGap=0. visitT=visits[(VisitDateTime < OCDate-tGap & VisitDateTime > OCDate-tWin-tGap ) , .(PatientICN,VisitDateTime)];colnames(visitT)[2]="OCDate" pts=visitT[,min(OCDate),by="PatientICN"];colnames(pts)[2]="Tpr" pts$Tpr=timePr rm(visits) } #Read demo file for gender, race, ethnicity, and complete list of patients in cohort. read in inpatient and outpatient Dx codes, the select cases by Dx codes demo=readRDS(paste(c(C2020dir,"DemoU.",cohort,".RDS"),collapse="")) fee9=readRDS(paste(c(C2020dir,"fee/FeeService.Dx9.",cohort,".RDS"),collapse="")) fee10=readRDS(paste(c(C2020dir,"fee/FeeService.Dx10.",cohort,".RDS"),collapse="")) Fx9=fee9[,c("PatientICN","InitialTreatmentDateTime","ICDCode")]; colnames(Fx9)[2]="VisitDateTime" Fx10=fee10[,c("PatientICN","InitialTreatmentDateTime","ICD10Code")]; colnames(Fx10)[2:3]=c("VisitDateTime","ICDCode") out9=readRDS(paste(c(C2020dir,"outpat/oDx9.",cohort,".RDS"),collapse="")) out10=readRDS(paste(c(C2020dir,"outpat/oDx10.",cohort,".RDS"),collapse="")) in9=readRDS(paste(c(C2020dir,"inpat/iDx9.",cohort,".RDS"),collapse="")) in10=readRDS(paste(c(C2020dir,"inpat/iDx10.",cohort,".RDS"),collapse="")) in9d=readRDS(paste(c(C2020dir,"inpat/idDx9.",cohort,".RDS"),collapse="")) in10d=readRDS(paste(c(C2020dir,"inpat/idDx10.",cohort,".RDS"),collapse="")) Dx9a=out9[,c("PatientICN","VisitDateTime","ICDCode")] Dx10a=out10[,c("PatientICN","VisitDateTime","ICD10Code")];colnames(Dx10a)[3]="ICDCode" Dx9b=in9[,c("PatientICN","DischargeDateTime","ICDCode")];colnames(Dx9b)[2]="VisitDateTime" Dx10b=in10[,c("PatientICN","DischargeDateTime","ICD10Code")];colnames(Dx10b)[2]="VisitDateTime";colnames(Dx10b)[3]="ICDCode" Dx9c=in9d[,c("PatientICN","AdmitDateTime","ICDCode")];colnames(Dx9c)[2]="VisitDateTime" Dx10c=in10d[,c("PatientICN","AdmitDateTime","ICD10Code")];colnames(Dx10c)[2]="VisitDateTime";colnames(Dx10c)[3]="ICDCode" mDx=rbind(Fx9,Fx10,Dx9a,Dx9b,Dx9c,Dx10a,Dx10b,Dx10c) Dx=mDx[,c("PatientICN","VisitDateTime","ICDCode")] Dx[,"lower":=tolower(ICDCode)] Dx[VisitDateTime < 15.75,"ICDCode"]=Dx[VisitDateTime < 15.75,"lower"] #ndi = readRDS(paste(c(C2020dir2,"NDI/ndi.",cohort,".RDS"),collapse="")) #ndi$OCDate=as.numeric(julian(as.Date(as.character(ndi$DOD_NDI),format="%m/%d/%Y"),origin=as.Date("2000-01-01")))/365.25 #ndi=ndi[MatchFound==1,] if (study=="Pds-e" | study=="Pds-c"){ Dxp=merge(Dx,pts,all.x=T,all.y=F,by="PatientICN") Dx=Dxp[is.na(Tpr)==F,] if(study=="Pds-c"){Dx$Tpr=timePr} DxBf=Dx[VisitDateTime-Tpr < 0.12,] Dx=Dx[VisitDateTime-Tpr >= 0.12,] ndip=merge(ndi,pts,all.x=T,all.y=F,by="PatientICN") ndi=ndip[is.na(Tpr)==F,] if(study=="Pds-c"){ndi$Tpr=timePr} ndiBf=ndi[OCDate-Tpr < 0.12, ] ndi=ndi[OCDate-Tpr >= 0.12,] } #SuicideAttempt=Dx[c(grep("^e95|^e98",ICDCode), grep("^T14.91|^T3[6789]|^T[456]...X2|^X[78]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] #SuicideIdeation=Dx[c(grep("^V62.84",ICDCode), grep("^R45.851",ICDCode)),.(PatientICN,VisitDateTime)] tmp=Dx[c(grep("^e95",ICDCode), grep("^T14.91|^T3[6789]..X2|^T[456]...X2|^X[678]..XX|^Y87.0",ICDCode)),.(PatientICN,ICDCode,VisitDateTime)] #left off E98 SuicideAttempt=tmp[,min(VisitDateTime),by="PatientICN"] ; SAtmp=tmp tmp=Dx[c(grep("^e850|^e9[67]|^e935",ICDCode), grep("^T40",ICDCode)),.(PatientICN,ICDCode,VisitDateTime)] Overdose=tmp[,min(VisitDateTime),by="PatientICN"] ; ODtmp=tmp suicide=unique( ndi[c(grep("^X[678]|^T3[6789]..X2|^T[456]...X2|^Y87.0|^X[78]",UnderlyingCause_NDI),grep("^X[678]|^T3[6789]..X2|^T[456]...X2|^Y87.0|^X[78]",rec_cond1),grep("^X[678]|^T3[6789]..X2|^T[456]...X2|^Y87.0|^X[78]",rec_cond2),grep("^X[678]|^T3[6789]..X2|^T[456]...X2|^Y87.0|^X[78]",rec_cond3),grep("^X[678]|^T3[6789]..X2|^T[456]...X2|^Y87.0|^X[78]",rec_cond4),grep("^X[678]|^T3[6789]..X2|^T[456]...X2|^Y87.0|^X[78]",rec_cond5),grep("^X[678]|^T3[6789]..X2|^T[456]...X2|^Y87.0|^X[78]",rec_cond6)),.(PatientICN, DOD_NDI)]) #suicide=unique( ndi[c(grep("^X[678]",UnderlyingCause_NDI),grep("^X[678]",rec_cond1),grep("^X[678]",rec_cond2),grep("^X[678]",rec_cond3),grep("^X[678]",rec_cond4),grep("^X[678]",rec_cond5),grep("^X[678]",rec_cond6)),.(PatientICN, DOD_NDI)]) suicide = suicide[is.na(match(PatientICN, unique(Dx$PatientICN)))==F,] # remove some extra patients from ndi that are not in Dx colnames(suicide)[2]="OCDate";colnames(SuicideAttempt)[2]="OCDate";colnames(Overdose)[2]="OCDate" suicide$OCDate=as.numeric(julian(as.Date(as.character(suicide$OCDate),format="%m/%d/%Y"),origin=as.Date("2000-01-01")))/365.25 #save(SAtmp,ODtmp,suicide,file=paste(c("subOC.",cohort,".RData"),collapse='')) OC.S=suicide;OC.SA=SuicideAttempt;OC.OD=Overdose OC.S$suicide=1;colnames(OC.S)[2]="S.OCDate" OC.SA$SA=1;colnames(OC.SA)[2]="SA.OCDate" OC.OD$OD=1;colnames(OC.OD)[2]="OD.OCDate" tmp=merge(OC.S,OC.SA,by="PatientICN",all.x=T,all.y=T) subOC=merge(OC.OD,tmp,by="PatientICN",all.x=T,all.y=T) subOC$OD[is.na(subOC$OD)==T]=0 subOC$SA[is.na(subOC$SA)==T]=0 subOC$suicide[is.na(subOC$suicide)==T]=0 setcolorder(subOC,c(1,5,4,7,6,3,2)) #saveRDS(subOC,file=paste(c("subOC.",cohort,".RDS"),collapse="")) #for (cohort in seq(1,99)){load(paste(c("subOC.",cohort,".RData"),collapse='')) #suicide$ICDCode="NDI";suicide$OC=1;SAtmp$OC=2;ODtmp$OC=3;s=suicide[,c(1,3,2,4)];suicide=s;colnames(SAtmp)[3]="OCDate";colnames(ODtmp)[3]="OCDate";ocC=rbind(suicide,SAtmp,ODtmp);ocC$C=cohort #OCAll=rbind(OCAll,ocC) # outcome=stroke print("total Pts with Dx, suicide, SuicideAttempt, Overdose, ") print(dim(unique(Dx,by="PatientICN"))) print(dim(unique(suicide,by="PatientICN"))) print(dim(unique(SuicideAttempt,by="PatientICN"))) print(dim(unique(Overdose,by="PatientICN"))) suicide$OC=1;SuicideAttempt$OC=2;Overdose$OC=3 m=match(SuicideAttempt$PatientICN,suicide$PatientICN) w=which(is.na(m)==T) outcome=rbind(suicide,SuicideAttempt[w,]) m=match(Overdose$PatientICN,outcome$PatientICN) w=which(is.na(m)==T) outcome=rbind(outcome,Overdose[w,]) print(c(cohort,dim(outcome)[1])) #write.csv(outcome,file=paste(c("outcome.",cohort,".csv"),collapse=""),row.names=F) if (study=="Pds-e" | study == "Pds-c") { ptsDemo=merge(pts,demo,all.x=T,all.y=F,by="PatientICN") # ptsDemo[,"BirthDate":=as.numeric(julian(BirthDateTime,origin=as.Date("2000-01-01")))/365.25] xOC=merge(ptsDemo,outcome,all.x=T,all.y=F,by="PatientICN") xOC[is.na(OC)==T,"OC"]=0 xOC[OC==0,"OCDate"]=xOC[OC==0,"LastDate"] xOC[is.na(OCDate)==T & is.na(DeathDateTime)==F,"OCDate"]=xOC[is.na(OCDate)==T & is.na(DeathDateTime)==F,"DeathDateTime"] xOC[is.na(OCDate)==T & is.na(DeathDateTime)==T,"OCDate"]=22.2 xOC$Tpr=xOC$Tpr+.02 #ad 1 week to ensure observations related to visit are recorded in first time bin xOC[is.na(OC)==T,"OC"]=0 xOC$Cohort=cohort OC=xOC[,c("PatientICN","Cohort","BirthDate","Gender","Race","Ethnicity","Height","Tpr","OCDate","OC")] if (CntF < 1.0){ y1=OC$OC y1[y1>0]=1 wY=which(y1>0) ; wZ=sample(which(y1==0),CntF*(dim(OC)[1]-length(wY))) ; w=c(wY,wZ) OC=OC[w,]} } if (study=="Rcc"){ CaOC=outcome[OCDate > 7.62,] AllCases=unique(outcome$PatientICN) cases=unique(CaOC$PatientICN) excl=AllCases[is.na(match(AllCases,cases))==T] #create table of visits by patients who have no noted suicide attempts and choose randomly among visits (weights patients with more visits more highly). Define OCDate also for controls, using randomly selected visit. Odx=Dx[is.na(match(PatientICN,cases))==T,];colnames(Odx)[2]="VisitDateTime" controls=sample(unique(Odx$PatientICN,by="PatientICN"),5*length(cases)) # pay attention to set.seed() at start of script. for repeatable results, keep first call to sample() Cadx=Dx[is.na(match(PatientICN,cases))==F,];colnames(Cadx)[2]="VisitDateTime" Cndx=Dx[is.na(match(PatientICN,controls))==F,];colnames(Cndx)[2]="VisitDateTime" CnOC=Dx[is.na(match(PatientICN,controls))==F,.SD[sample(.N,1)],by=PatientICN] CnOC$ICDCode=NULL colnames(CnOC)[2]="OCDate" CnOC$OC=0 CnOC=CnOC[OCDate > 7.5,] Cdx=rbind(Cndx,Cadx) COC=rbind(CaOC,CnOC) #COC$OC=c(rep(1,dim(CaOC)[1]),rep(0,dim(CnOC)[1])) DmOC=demo[COC] #DmOC[,"BirthDate":=as.numeric(julian(BirthDateTime,origin=as.Date("2000-01-01")))/365.25] DmOC[,"Cohort":=cohort] DmOC$"Tpr"=DmOC$OCDate-0.12 OC=DmOC[,c("PatientICN","Cohort","BirthDate","Gender","Race","Ethnicity","Height","Tpr","OCDate","OC")] } pts=OC$PatientICN if (study == "Pds-e" | study == "Pds-c"){Dx=DxBf;Dx$Tpr=NULL} else {Dx=Cdx[is.na(match(PatientICN,pts))==F,]} #tGap=0.25 parameters below determine start and end of prediction windows, used in #[VisitDateTime < OCDate-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"], where tGap comes from gindow and tWin from window #gindow=c(0.12,0.37,0.62,1.12,1.62,2.62,3.62,5.62,7.62) gindow=c(0.,0.25,0.5,1.,1.5,2.5,3.5,5.5) window=c(0.25,0.25,0.5,0.5,1,1,2,2) ####### at this point, we have OC for outcomes (one row for patient) and Dx, with all dx codes for patients and controls. #Reachvet2 variables #BirthDateTime Gender #agegroup_ls30 agegroup_ge30 agegroup_ge40 agegroup_ge50 agegroup_ge60 agegroup_ge70 oDemo=readRDS(paste(c(C2020dir,"outpat/oDemo.",cohort,".RDS"),collapse="")) moDemo=merge(OC,oDemo,all.x=T,by="PatientICN") oDemo=moDemo for (count in c(1,2,3,4,5,6,7,8)){ tWin=window[count];tGap=gindow[count] MaritalStatus=oDemo[(VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap), Mode(PatientMaritalStatus),by="PatientICN"];colnames(MaritalStatus)[2]="MaritalStatus" PatientFIPS=oDemo[(VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap), Mode(PatientFIPS),by="PatientICN"];colnames(PatientFIPS)[2]="PatientFIPS" PatientZIP=oDemo[(VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap), Mode(PatientZIP),by="PatientICN"];colnames(PatientZIP)[2]="PatientZIP" rvDemo=list(OC,MaritalStatus,PatientFIPS,PatientZIP) dDm=Reduce(function(...) merge(...,all=T,by="PatientICN"),rvDemo) #for (i in seq(10,10)){dDm[which(is.na(dDm[,..i])),i]=0} #make sure columns 3-16 are the counts per condition columns. replace NA with 0 #for (i in seq(9,84)){d1[which(d1[,..i] >=1),i]=1} if (count==1){dDemo1=dDm} if (count==2){dDemo2=dDm} if (count==3){dDemo3=dDm} if (count==4){dDemo4=dDm} if (count==5){dDemo5=dDm} if (count==6){dDemo6=dDm} if (count==7){dDemo7=dDm} if (count==8){dDemo8=dDm} } #MaritalStatusCode White_ serviceconnectedgroup region orh mst_status #mar_Unknown mar_Uncopled mar_Never_married mar_Married mst_unknown mst_no mst_declined mst_yes #orh_unknown orh_Rural region_unknown region_south region_northeast region_midwest region_west #serviceconnection_ge50 serviceconnection_ls_50 serviceconnection_0 sex_female white_unknown white_non_white #Diagnosis #Dx variables, ICD9 | ICD10 #moFac=merge(OC,oFac,all.y=T,by="PatientICN") Dim.Dx9=readRDS(paste(c(C2020dir,"Dim.ICD9.RDS"),collapse="")) Dim.Dx10=readRDS(paste(c(C2020dir,"Dim.ICD10.RDS"),collapse=""));colnames(Dim.Dx10)=c("ICDCode","ICDDescription") Dim.Dx=rbind(Dim.Dx9,Dim.Dx10);setkey(Dim.Dx,ICDCode) #use, eg. tmp=Dim.Dx[SuicideAttempt] mDx=merge(OC,Dx,all.x=T,by="PatientICN") Dx=mDx Dx[,"lower":=tolower(ICDCode)] Dx[VisitDateTime < 15.75,"ICDCode"]=Dx[VisitDateTime < 15.75,"lower"] SuicideAttempt=Dx[c(grep("^e95|^e98",ICDCode), grep("^T14.91|^T3[6789]..X2|^T[456]...X2|^X[678]..XX|^Y87.0|^X[78]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] SuicideIdeation=Dx[c(grep("^v62.84",ICDCode), grep("^R45.851",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] alcdx_poss=Dx[c(grep("^291|^303|^305.0",ICDCode), grep("^F10.",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] #alcdx_poss 291%' '303%' '305.0% Nicdx_poss=Dx[c(grep("^305.1",ICDCode), grep("^F17.2|^Z72.0",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] AmphetamineUseDisorder=Dx[c(grep("^304.4|^305.7",ICDCode), grep("^F15.",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Cannabis=Dx[c(grep("^304.3|^305.2",ICDCode), grep("^F12",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] COCNdx=Dx[c(grep("^304.2|^305.6",ICDCode), grep("^F14",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] DrugInducedSymptoms=Dx[c(grep("^292",ICDCode), grep("^F19",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] OpioidOverdose=Dx[c(grep("^965.0|e850.0|e850.1|e850.2|e980.0|e935.0|e935.1|e935.2",ICDCode), grep("^T40.[01234]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] OUD=Dx[c(grep("^304.0|305.5|304.7",ICDCode), grep("^F11",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] SedativeUseDisorder=Dx[c(grep("^304.1|^305.4",ICDCode), grep("^F13",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] SAE_sed=Dx[c(grep("^967.[08]|^968.0|^969.[12345]|^e851.|^e852.[01234589]|^e853.[01289]|^e937.[08]|^e938.0|^e939.[1245]|^e980.[123]",ICDCode), grep("^T42",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] SUD_CatchAll=Dx[c(grep("^304.[5689|^305.[389]",ICDCode), grep("^F1[1235689]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] SAE_Acet=Dx[c(grep("e850.4|e935.4|^967.0|^967.8|^968.0|^968.1|^968.2|^968.3|^968.4|^968.5|e851.|e852.[01234589]|e853.[01289]|e937.[089]|e939.[1245]|^e980.[13]",ICDCode), grep("^T39",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] SAE_OtherDrug=Dx[c(grep("^e855.1|969.05|e855.0|e854.1|e935.6|e855.6|e939.7|e980.5|e935.3|e855.5|969.02|e850.3|e940.1|e980.4|e855.9|e939.0|965.1|e855.4|e854.0|^969.6|^965.69|^970.1|e855.3|e939.6|^969|^965.6|^969.72|^969.03|e855.8|e855.2|^965.61|e854.3|^969.04",ICDCode), grep("^T40.[5789]|^T4[13]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] MDD=Dx[c(grep("^296.[23]",ICDCode), grep("^F33|^F32.[012345679]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Schiz=Dx[c(grep("295.0|295.1|295.2|295.3|295.4|295.6|295.7|295.8|295.9",ICDCode), grep("^F20|^F21|^F25",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Sleep=Dx[c(grep("307.4|291.82|327.0|327.1|327.2|780.50|780.52|780.54|780.55|780.56|780.58|780.59",ICDCode), grep("^G47|^G51",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] BIPOLAR=Dx[c(grep("^296.[0145678]",ICDCode), grep("^F30|F31",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] PTSD=Dx[c(grep("309.81",ICDCode), grep("^F43.1",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] MH_CatchAll=Dx[c(grep("^290.[89]|^293.[01234]|^295.00|^295.5|^296.82|^296.9|^300.1[1234569]|^300.[3456789]|^301.[13]|^301.21|^301.5[19]|^301.8[49]|^307.1|^307.5[12349]|^308|^309.[012349]|^309.8[239]|^31[124]|^780.1",ICDCode), grep("^F3[234]|^4[45]|^F6[38]|^F9[01]|^R45.7",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] ODEPRdx_poss=Dx[c(grep("^301.12|300.4|^293.83|^298.0|^301.1|^311|^296.9|^309.0|^309.1|^296.82",ICDCode), grep("^F32|^F33",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] OtherPsychosis=Dx[c(grep("^29[78]",ICDCode), grep("^F2[23489]|^F06.[02]|^F53.[01]|^F53|^F06.3.",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Thyroid=Dx[c(grep("^24[0123456]",ICDCode), grep("^E0[0123]|^E89.0",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] reach_othanxdis=Dx[c(grep("^300.[02]|^300.10",ICDCode), grep("^F40|^F41",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] reach_persond=Dx[c(grep("^301.0|^301.20|^301.22|^301.7|^301.83|^301.50|^301.81|^301.82|^301.6|^301.4|^301.9",ICDCode), grep("^F60",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] TBI_Dx=Dx[c(grep("^85[123456789]",ICDCode), grep("^S06.[23456]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Concuss=Dx[c(grep("^850",ICDCode), grep("^S06.[019]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Headache=Dx[c(grep("^784.0|^307.81|^339|^346",ICDCode), grep("^G4[34]|^M54.81|^R51",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] SAE_Falls=Dx[c(grep("^e888.|^e880.0|^e987.9|^e888.9|^e987.2|^e882.|^e883.9|^e881.1|^e885.1|^e885.2|^e888.0|^e884.9|^e884.2|^e884.0|^e987.0|^e885.9|^e884.3|^e887.|^e883.0|^e885.4|^e886.9|^e885.|^e883.2|^e880.1|^e880.9|^e888.8|^e885.3|^e884.5|^e886.0|^e888.1|^e884.1|^e987.1|^e884.6|^e884.4|^e929.3|^e881.0|^e885.0|^e883.1",ICDCode), grep("^W[01]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] SAE_OtherAccident=Dx[c(grep("e920.4|^e920.9|^e922.5|^e922.8|^e985.6|^e922.9|^e919.5|^e910.4|^e985.7|^e920.2|^e919.0|^e919.2|^e922.3|^e910.0|^e920.3|^e919.3|^e920.5|^e985.4|^e922.0|^e910.1|^e920.8|^e985.5|^e985.3|^e920.0|^e910.2|^e919.8|^e922.1|^e910.9|^e985.2|^e919.1|^e922.4|^e919.4|^e985.1|^e919.6|^e910.3|^e985.0|^e920.1|^e919.7|^e986.|^e919.9|^e922.2|^e910.8",ICDCode), grep("^W[23467]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] SAE_Vehicle=Dx[c(grep("^e836.9|^e835.6|^e842.6|^e827.0|^e833.7|^e811.3|^e820.3|^e819.6|^e825.4|^e841.4|^e830.2|^e835.4|^e815.4|^e816.1|^e818.6|^e846.|^e840.4|^e827.8|^e814.8|^e929.0|^e810.7|^e817.2|^e838.9|^e823.5|^e826.1|^e833.3|^e821.1|^e816.6|^e813.0|^e844.7|^e834.8|^e817.8|^e832.5|^e843.8|^e838.5|^e833.1|^e821.3|^e811.0|^e822.8|^e824.5|^e848.|^e818.0|^e805.3|^e837.5|^e806.8|^e812.7|^e831.2|^e821.7|^e813.6|^e824.3|^e844.5|^e828.2|^e814.3|^e821.8|^e838.6|^e802.8|^e836.1|^e818.7|^e836.6|^e835.1|^e828.9|^e834.1|^e802.2|^e840.7|^e817.5|^e825.3|^e841.5|^e830.3|^e825.9|^e800.0|^e829.4|^e807.1|^e840.1|^e824.8|^e818.1|^e837.4|^e810.4|^e833.0|^e838.2|^e821.4|^e817.9|^e826.4|^e812.3|^e843.3|^e811.5|^e831.9|^e832.0|^e805.0|^e830.7|^e820.4|^e843.5|^e804.8|^e820.6|^e824.6|^e836.4|^e805.9|^e820.1|^e823.9|^e834.0|^e836.2|^e800.9|^e802.9|^e814.2|^e814.1|^e833.9|^e825.8|^e836.7|^e816.3|^e828.8|^e823.7|^e826.9|^e815.0|^e842.9|^e832.7|^e817.4|^e840.6|^e800.1|^e818.8|^e830.0|^e830.6|^e824.9|^e810.2|^e825.2|^e810.5|^e807.0|^e819.9|^e838.1|^e833.5|^e811.4|^e843.6|^e804.3|^e819.7|^e806.0|^e845.0|^e844.3|^e834.2|^e813.8|^e837.9|^e821.9|^e804.9|^e812.9|^e837.3|^e818.9|^e841.7|^e828.4|^e803.9|^e801.9|^e820.2|^e834.3|^e811.2|^e833.8|^e822.4|^e829.8|^e801.3|^e827.3|^e829.0|^e801.1|^e825.7|^e814.7|^e817.7|^e816.4|^e816.0|^e823.4|^e815.7|^e823.32|^e802.0|^e840.3|^e810.8|^e817.1|^e814.5|^e830.1|^e845.9|^e826.2|^e845.8|^e825.1|^e816.9|^e819.4|^e831.7|^e832.2|^e834.5|^e813.1|^e844.2|^e838.4|^e807.9|^e833.2|^e843.1|^e804.0|^e811.7|^e805.2|^e813.5|^e813.7|^e844.6|^e844.8|^e837.8|^e837.2|^e824.0|^e806.9|^e831.5|^e842.7|^e803.3|^e814.4|^e835.8|^e836.0|^e801.8|^e803.8|^e811.1|^e820.5|^e834.4|^e829.9|^e822.3|^e810.0|^e817.6|^e816.7|^e827.2|^e833.4|^e818.2|^e815.8|^e823.1|^e807.2|^e802.1|^e825.0|^e832.9|^e841.6|^e830.4|^e840.0|^e804.1|^e819.3|^e837.1|^e816.8|^e812.0|^e807.8|^e821.5|^e838.3|^e843.4|^e831.8|^e802.3|^e806.2|^e819.5|^e824.7|^e830.8|^e805.1|^e831.0|^e803.2|^e805.8|^e820.9|^e844.9|^e813.2|^e831.6|^e836.3|^e814.9|^e820.0|^e814.0|^e835.5|^e800.2|^e818.3|^e816.2|^e841.1|^e835.3|^e815.5|^e838.8|^e810.6|^e825.5|^e810.3|^e823.6|^e811.9|^e840.8|^e819.2|^e806.1|^e815.2|^e801.0|^e832.4|^e841.0|^e834.7|^e844.0|^e807.3|^e803.1|^e822.7|^e822.1|^e837.0|^e819.0|^e812.4|^e821.6|^e843.9|^e831.3|^e826.0|^e840.9|^e813.9|^e837.6|^e822.5|^e812.6|^e813.4|^e800.3|^e836.5|^e816.5|^e814.6|^e825.6|^e841.2|^e810.9|^e815.6|^e827.4|^e835.2|^e840.2|^e838.7|^e841.8|^e835.0|^e817.0|^e818.4|^e843.2|^e812.2|^e834.6|^e820.7|^e819.1|^e826.3|^e811.8|^e832.1|^e844.1|^e843.0|^e801.2|^e811.6|^e837.7|^e822.6|^e832.3|^e824.1|^e831.4|^e812.5|^e803.0|^e842.8|^e830.5|^e836.8|^e835.7|^e841.9|^e823.8|^e822.0|^e822.2|^e828.0|^e800.8|^e840.5|^e815.3|^e815.1|^e827.9|^e817.3|^e810.1|^e818.5|^e823.0|^e832.8|^e815.9|^e841.3|^e806.3|^e804.2|^e834.9|^e823.3|^e838.0|^e833.6|^e822.9|^e821.2|^e843.7|^e826.8|^e832.6|^e812.1|^e835.9|^e819.8|^e812.8|^e831.1|^e824.4|^e821.0|^e847.|^e830.9|^e813.3|^e820.8|^e824.2|^e844.4",ICDCode), grep("^V[2345678]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Amputation=Dx[c(grep("^885|^886|^887|^895|^896|^897|^V49.6|^V49.7",ICDCode), grep("^S68|^S88",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] SpinalCordInj=Dx[c(grep("^806.[0123456789]|^952.[0123489]",ICDCode), grep("^S14|^S24|^S34|^G82",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Backpain=Dx[c(grep("^720|^721|^722|^723|^724",ICDCode), grep("^G54.[134]|^M43.[013589]|^M45.[0456789]|^M46.[01489]|^M47|^M48|^M49|^M51|^M53|^M54|^Q76|^S23|^S33",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_PARALYSIS=Dx[c(grep("^342.[019]|^34[34]",ICDCode), grep("^G04.1|^G11.4|^G80.[12]|^G8[123]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] DEMENTIA=Dx[c(grep("^046.[13]|^046.79|^290.[01234]|^291.2|^292.82|^294.[18]|^331.[01279]|^331.8[23]|^333.[04]|^294.2",ICDCode), grep("^A81.[0289]|^F0[123]|^F10.29]|^F13.[29]|^F18.[129]|^F19.[129]|^G23.1|^G30.[0189]|G31.09|G31.83|G90.3",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Parkinsons=Dx[c(grep("^332",ICDCode), grep("^G20",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] DeliriumTremens=Dx[c(grep("^291",ICDCode), grep("^F10.2[3567]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Huntington=Dx[c(grep("^333.4",ICDCode), grep("^G10",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_Lymphoma=Dx[c(grep("^238.6|^273.3|^202.[0123]|^202.[56789]|^203.01|^203.8[01]|^v10.7",ICDCode), grep("^C8[1234568]|^C90.[02]|^96",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_METCANCR=Dx[c(grep("^19[678]|^199.[01]",ICDCode), grep("^C80|^C7[6789]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_NMETTUMR=Dx[c(grep("^17[459]|^1[456]|^17[0129]|^18|^19[012345]",ICDCode), grep("^C[01234567]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] HIV=Dx[c(grep("^04[234]|^v08",ICDCode), grep("^B20|^Z21",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Tourette=Dx[c(grep("^307.23",ICDCode), grep("^F95.2",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Osteoporosis=Dx[c(grep("^733.0",ICDCode), grep("^M8[13]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Fibromyalgia=Dx[c(grep("^729.1",ICDCode), grep("^D86.87|^M60.8|^M79.[17]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Hearing=Dx[c(grep("^389",ICDCode), grep("^H90",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Visual=Dx[c(grep("^369",ICDCode), grep("^H5[34]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Neuro=Dx[c(grep("^256.60|^355.[09]|^357.2|^357.9|^356",ICDCode), grep("^B02.2[23]|E08.42|^E10.4[0123]|E11.4|E13.4|^G50.[089]|^G52.1|^G54.[89]|^G56.[89]|^G57.[89]|^G58.[789]|^G60.[2389]|^G61.9|^G62.[089]|^G90.09|^G99.0|^H46.[2389]|^H47.0|^M79.2",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] reach_chronic=Dx[c(grep("^338.0|^338.2|^338.4",ICDCode), grep("^G89.2",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] GastricUlcer=Dx[c(grep("^53[1234].[5679]",ICDCode), grep("^K2[5678].[79]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] SleepApnea=Dx[c(grep("^327.2|^780.5",ICDCode), grep("^G47.3",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] Homeless=Dx[c(grep("^v60.0",ICDCode), grep("^Z59.0",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] AMI=Dx[c(grep("^41[02]",ICDCode), grep("^I2[12]|I25.2",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] CAD=Dx[c(grep("^411|^413|^414",ICDCode), grep("^I25",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] ChestPain_Dx=Dx[c(grep("^786.50",ICDCode), grep("^R07.[189]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_ARRHYTH=Dx[c(grep("^426.10|^426.11|^426.13|^427.0|^427.2|^427.31|^427.60|^427.9|^785.0|^v45.0|^v53.3|^426.[234678]|^426.5[0123]",ICDCode), grep("^I4[89]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_BLANEMIA=Dx[c(grep("^280.0",ICDCode), grep("^D5[89]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_CHRNPULM=Dx[c(grep("^49|^506.4|^50[012345]",ICDCode), grep("^J4[01234]|^J6",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_COAG=Dx[c(grep("^286|^287.[1345]",ICDCode), grep("^D6[56789]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_COMDIAB=Dx[c(grep("^250.[45679]",ICDCode), grep("^E1[01].[012345]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_DefANEMIA=Dx[c(grep("^280.[189]|^281",ICDCode), grep("^D50.[89]|^D51.[012389]|^D52.[0189]|6D53.[01289]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_ELECTRLYTE=Dx[c(grep("^276",ICDCode), grep("^E22.2|^E86.[019]|^E87.[012345678]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_HEART=Dx[c(grep("^398.91|^402.11|^402.91|^404.1[13]|^404.9[13]|^428",ICDCode), grep("^I50",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_HYPERTENS=Dx[c(grep("^401.[19]|^402.[19]0|^404.[19]0|^405.[19]",ICDCode), grep("^I1[01]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_LIVER=Dx[c(grep("^070.3[23]|^070.54|^456.[012]|^571.[0235689]|^571.4[09]|^572.[38]|^v42.7",ICDCode), grep("^K7[0123456]|^B18.|^I85",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_OBESITY=Dx[c(grep("^278.0",ICDCode), grep("^E66",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_OTHNEURO=Dx[c(grep("^331.9|^332.0|^333.[45]|^340|^33[45]|^341.[1289]|^345.[014589]",ICDCode), grep("^G1[0123]|^G2[015]|^G31.[289]|^G32.[089[|^G35|^G3.[189|^G37.[01234589]|^G40.01234589AB]|^G.[14]|^R47|^R56.[019]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_PERIVASC=Dx[c(grep("^441.[2479]|^447.1|^557.[19]|^v43.4|^440|^443.[123456789]",ICDCode), grep("^I70.[12345678]|^I71.|^I72.|^I73.9",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_PULMCIRC=Dx[c(grep("^416|^417.9",ICDCode), grep("^I27",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_RENAL=Dx[c(grep("^403.11|^403.91|^404.[19]2|^58[56]|^v42.0|^v45.1|^v56.[48]",ICDCode), grep("^I12.0|^I13.1|^N18.[1234569]|^N19|^N25.0|^Z49.0|^Z94.0|^Z99.2",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_RHEUMART=Dx[c(grep("^701.0|^71[04]|^72[05]",ICDCode), grep("^M0[23456]|^M1[23]|^M45|^M32",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_UNCDIAB=Dx[c(grep("^250.0",ICDCode), grep("^E1[013].[39]|^E08.[89]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_VALVDIS=Dx[c(grep("^v42.4|^v43.3|^093.2|^39[456]|^397.[01]|^424.[0-8]|^424.9[01]|^746.[3456]",ICDCode), grep("^A52.0|^I0[56789]|^I3[456789]|^Q23.[0123]|^Z95.[234]",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] EH_WEIGHTLS=Dx[c(grep("^26[0123]|^783.21",ICDCode), grep("^E4[0123456]|^R63.4|^R64",ICDCode)),.(PatientICN,VisitDateTime,Tpr,ICDCode)] #pdf("conditions.pdf") #cond(SuicideAttempt,"SuicideAttempt");cond(SuicideIdeation,"SuicideIdeation");cond(alcdx_poss,"alcdx_poss");cond(Nicdx_poss,"Nicdx_poss");cond(AmphetamineUseDisorder,"AmphetamineUseDisorder");cond(Cannabis,"Cannabis");cond(COCNdx,"COCNdx");cond(DrugInducedSymptoms,"DrugInducedSymptoms");cond(OpioidOverdose,"OpioidOverdose");cond(OUD,"OUD");cond(SedativeUseDisorder,"SedativeUseDisorder");cond(SAE_sed,"SAE_sed");cond(SUD_CatchAll,"SUD_CatchAll");cond(SAE_Acet,"SAE_Acet");cond(SAE_OtherDrug,"SAE_OtherDrug");cond(MDD,"MDD");cond(Schiz,"Schiz");cond(Sleep,"Sleep");cond(BIPOLAR,"BIPOLAR");cond(PTSD,"PTSD");cond(MH_CatchAll,"MH_CatchAll");cond(ODEPRdx_poss,"ODEPRdx_poss");cond(OtherPsychosis,"OtherPsychosis");cond(Thyroid,"Thyroid");cond(reach_othanxdis,"reach_othanxdis");cond(reach_persond,"reach_persond");cond(TBI_Dx,"TBI_Dx");cond(Concuss,"Concuss");cond(Headache,"Headache");cond(SAE_Falls,"SAE_Falls");cond(SAE_OtherAccident,"SAE_OtherAccident");cond(SAE_Vehicle,"SAE_Vehicle");cond(Amputation,"Amputation");cond(SpinalCordInj,"SpinalCordInj");cond(Backpain,"Backpain");cond(EH_PARALYSIS,"EH_PARALYSIS");cond(DEMENTIA,"DEMENTIA");cond(Parkinsons,"Parkinsons");cond(DeliriumTremens,"DeliriumTremens");cond(Huntington,"Huntington");cond(EH_Lymphoma,"EH_Lymphoma");cond(EH_METCANCR,"EH_METCANCR");cond(EH_NMETTUMR,"EH_NMETTUMR");cond(HIV,"HIV");cond(Tourette,"Tourette");cond(Osteoporosis,"Osteoporosis");cond(Fibromyalgia,"Fibromyalgia");cond(Hearing,"Hearing");cond(Visual,"Visual");cond(Neuro,"Neuro");cond(reach_chronic,"reach_chronic");cond(GastricUlcer,"GastricUlcer");cond(SleepApnea,"SleepApnea");cond(Homeless,"Homeless");cond(AMI,"AMI");cond(CAD,"CAD");cond(ChestPain_Dx,"ChestPain_Dx");cond(EH_ARRHYTH,"EH_ARRHYTH");cond(EH_BLANEMIA,"EH_BLANEMIA");cond(EH_CHRNPULM,"EH_CHRNPULM");cond(EH_COAG,"EH_COAG");cond(EH_COMDIAB,"EH_COMDIAB");cond(EH_DefANEMIA,"EH_DefANEMIA");cond(EH_ELECTRLYTE,"EH_ELECTRLYTE");cond(EH_HEART,"EH_HEART");cond(EH_HYPERTENS,"EH_HYPERTENS");cond(EH_LIVER,"EH_LIVER");cond(EH_OBESITY,"EH_OBESITY");cond(EH_OTHNEURO,"EH_OTHNEURO");cond(EH_PERIVASC,"EH_PERIVASC");cond(EH_PULMCIRC,"EH_PULMCIRC");cond(EH_RENAL,"EH_RENAL");cond(EH_RHEUMART,"EH_RHEUMART");cond(EH_UNCDIAB,"EH_UNCDIAB");cond(EH_VALVDIS,"EH_VALVDIS");cond(EH_WEIGHTLS,"EH_WEIGHTLS") #dev.off() # #cond=function(cond,title){ # a=Dim.Dx[match(toupper(names(head(sort(table(cond$ICDCode),decreasing=T),10))),Dim.Dx$ICDCode),] ; b=head(sort(table(cond$ICDCode),decreasing=T),10) # provide list of top codes # df=data.table(b,a); df$V1=NULL # h1=hist(cond$VisitDateTime,breaks=240,main=title,xlim=c(0,22.5)) #for (i in seq(1,10)){text(0,max(h1$counts)*(0.95-i*.05),paste(df[i,],collapse=": "),adj=0,cex=.7)} #} ############# count instances for (count in c(1,2,3,4,5,6,7,8)){ tWin=window[count];tGap=gindow[count] SuicideAttempt.N=SuicideAttempt[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(SuicideAttempt.N)[2]="SuicideAttempt" SuicideIdeation.N=SuicideIdeation[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(SuicideIdeation.N)[2]="SuicideIdeation" alcdx_poss.N=alcdx_poss[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(alcdx_poss.N)[2]="alcdx_poss" Nicdx_poss.N=Nicdx_poss[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Nicdx_poss.N)[2]="Nicdx_poss" AmphetamineUseDisorder.N=AmphetamineUseDisorder[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(AmphetamineUseDisorder.N)[2]="AmphetamineUseDisorder" Cannabis.N=Cannabis[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Cannabis.N)[2]="Cannabis" COCNdx.N=COCNdx[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(COCNdx.N)[2]="COCNdx" DrugInducedSymptoms.N=DrugInducedSymptoms[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(DrugInducedSymptoms.N)[2]="DrugInducedSymptoms" OpioidOverdose.N=OpioidOverdose[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(OpioidOverdose.N)[2]="OpioidOverdose" OUD.N=OUD[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(OUD.N)[2]="OUD" SedativeUseDisorder.N=SedativeUseDisorder[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(SedativeUseDisorder.N)[2]="SedativeUseDisorder" SAE_sed.N=SAE_sed[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(SAE_sed.N)[2]="SAE_sed" SUD_CatchAll.N=SUD_CatchAll[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(SUD_CatchAll.N)[2]="SUD_CatchAll" SAE_Acet.N=SAE_Acet[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(SAE_Acet.N)[2]="SAE_Acet" SAE_OtherDrug.N=SAE_OtherDrug[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(SAE_OtherDrug.N)[2]="SAE_OtherDrug" MDD.N=MDD[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(MDD.N)[2]="MDD" Schiz.N=Schiz[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Schiz.N)[2]="Schiz" Sleep.N=Sleep[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Sleep.N)[2]="Sleep" BIPOLAR.N=BIPOLAR[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(BIPOLAR.N)[2]="BIPOLAR" PTSD.N=PTSD[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(PTSD.N)[2]="PTSD" MH_CatchAll.N=MH_CatchAll[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(MH_CatchAll.N)[2]="MH_CatchAll" ODEPRdx_poss.N=ODEPRdx_poss[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(ODEPRdx_poss.N)[2]="ODEPRdx_poss" OtherPsychosis.N=OtherPsychosis[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(OtherPsychosis.N)[2]="OtherPsychosis" Thyroid.N=Thyroid[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Thyroid.N)[2]="Thyroid" reach_othanxdis.N=reach_othanxdis[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(reach_othanxdis.N)[2]="reach_othanxdis" reach_persond.N=reach_persond[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(reach_persond.N)[2]="reach_persond" TBI_Dx.N=TBI_Dx[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(TBI_Dx.N)[2]="TBI_Dx" Concuss.N=Concuss[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Concuss.N)[2]="Concuss" Headache.N=Headache[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Headache.N)[2]="Headache" SAE_Falls.N=SAE_Falls[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(SAE_Falls.N)[2]="SAE_Falls" SAE_OtherAccident.N=SAE_OtherAccident[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(SAE_OtherAccident.N)[2]="SAE_OtherAccident" SAE_Vehicle.N=SAE_Vehicle[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(SAE_Vehicle.N)[2]="SAE_Vehicle" Amputation.N=Amputation[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Amputation.N)[2]="Amputation" SpinalCordInj.N=SpinalCordInj[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(SpinalCordInj.N)[2]="SpinalCordInj" Backpain.N=Backpain[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Backpain.N)[2]="Backpain" EH_PARALYSIS.N=EH_PARALYSIS[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_PARALYSIS.N)[2]="EH_PARALYSIS" DEMENTIA.N=DEMENTIA[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(DEMENTIA.N)[2]="DEMENTIA" Parkinsons.N=Parkinsons[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Parkinsons.N)[2]="Parkinsons" DeliriumTremens.N=DeliriumTremens[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(DeliriumTremens.N)[2]="DeliriumTremens" Huntington.N=Huntington[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Huntington.N)[2]="Huntington" EH_Lymphoma.N=EH_Lymphoma[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_Lymphoma.N)[2]="EH_Lymphoma" EH_METCANCR.N=EH_METCANCR[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_METCANCR.N)[2]="EH_METCANCR" EH_NMETTUMR.N=EH_NMETTUMR[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_NMETTUMR.N)[2]="EH_NMETTUMR" HIV.N=HIV[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(HIV.N)[2]="HIV" Tourette.N=Tourette[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Tourette.N)[2]="Tourette" Osteoporosis.N=Osteoporosis[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Osteoporosis.N)[2]="Osteoporosis" Fibromyalgia.N=Fibromyalgia[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Fibromyalgia.N)[2]="Fibromyalgia" Hearing.N=Hearing[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Hearing.N)[2]="Hearing" Visual.N=Visual[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Visual.N)[2]="Visual" Neuro.N=Neuro[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Neuro.N)[2]="Neuro" reach_chronic.N=reach_chronic[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(reach_chronic.N)[2]="reach_chronic" GastricUlcer.N=GastricUlcer[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(GastricUlcer.N)[2]="GastricUlcer" SleepApnea.N=SleepApnea[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(SleepApnea.N)[2]="SleepApnea" Homeless.N=Homeless[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Homeless.N)[2]="Homeless" AMI.N=AMI[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(AMI.N)[2]="AMI" CAD.N=CAD[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(CAD.N)[2]="CAD" ChestPain_Dx.N=ChestPain_Dx[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(ChestPain_Dx.N)[2]="ChestPain_Dx" EH_ARRHYTH.N=EH_ARRHYTH[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_ARRHYTH.N)[2]="EH_ARRHYTH" EH_BLANEMIA.N=EH_BLANEMIA[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_BLANEMIA.N)[2]="EH_BLANEMIA" EH_CHRNPULM.N=EH_CHRNPULM[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_CHRNPULM.N)[2]="EH_CHRNPULM" EH_COAG.N=EH_COAG[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_COAG.N)[2]="EH_COAG" EH_COMDIAB.N=EH_COMDIAB[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_COMDIAB.N)[2]="EH_COMDIAB" EH_DefANEMIA.N=EH_DefANEMIA[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_DefANEMIA.N)[2]="EH_DefANEMIA" EH_ELECTRLYTE.N=EH_ELECTRLYTE[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_ELECTRLYTE.N)[2]="EH_ELECTRLYTE" EH_HEART.N=EH_HEART[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_HEART.N)[2]="EH_HEART" EH_HYPERTENS.N=EH_HYPERTENS[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_HYPERTENS.N)[2]="EH_HYPERTENS" EH_LIVER.N=EH_LIVER[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_LIVER.N)[2]="EH_LIVER" EH_OBESITY.N=EH_OBESITY[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_OBESITY.N)[2]="EH_OBESITY" EH_OTHNEURO.N=EH_OTHNEURO[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_OTHNEURO.N)[2]="EH_OTHNEURO" EH_PERIVASC.N=EH_PERIVASC[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_PERIVASC.N)[2]="EH_PERIVASC" EH_PULMCIRC.N=EH_PULMCIRC[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_PULMCIRC.N)[2]="EH_PULMCIRC" EH_RENAL.N=EH_RENAL[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_RENAL.N)[2]="EH_RENAL" EH_RHEUMART.N=EH_RHEUMART[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_RHEUMART.N)[2]="EH_RHEUMART" EH_UNCDIAB.N=EH_UNCDIAB[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_UNCDIAB.N)[2]="EH_UNCDIAB" EH_VALVDIS.N=EH_VALVDIS[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_VALVDIS.N)[2]="EH_VALVDIS" EH_WEIGHTLS.N=EH_WEIGHTLS[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EH_WEIGHTLS.N)[2]="EH_WEIGHTLS" RVDx=list(OC,SuicideAttempt.N,SuicideIdeation.N,alcdx_poss.N,Nicdx_poss.N,AmphetamineUseDisorder.N,Cannabis.N,COCNdx.N,DrugInducedSymptoms.N,OpioidOverdose.N,OUD.N,SedativeUseDisorder.N,SAE_sed.N,SUD_CatchAll.N,SAE_Acet.N,SAE_OtherDrug.N,MDD.N,Schiz.N,Sleep.N,BIPOLAR.N,PTSD.N,MH_CatchAll.N,ODEPRdx_poss.N,OtherPsychosis.N,Thyroid.N,reach_othanxdis.N,reach_persond.N,TBI_Dx.N,Concuss.N,Headache.N,SAE_Falls.N,SAE_OtherAccident.N,SAE_Vehicle.N,Amputation.N,SpinalCordInj.N,Backpain.N,EH_PARALYSIS.N,DEMENTIA.N,Parkinsons.N,DeliriumTremens.N,Huntington.N,EH_Lymphoma.N,EH_METCANCR.N,EH_NMETTUMR.N,HIV.N,Tourette.N,Osteoporosis.N,Fibromyalgia.N,Hearing.N,Visual.N,Neuro.N,reach_chronic.N,GastricUlcer.N,SleepApnea.N,Homeless.N,AMI.N,CAD.N,ChestPain_Dx.N,EH_ARRHYTH.N,EH_BLANEMIA.N,EH_CHRNPULM.N,EH_COAG.N,EH_COMDIAB.N,EH_DefANEMIA.N,EH_ELECTRLYTE.N,EH_HEART.N,EH_HYPERTENS.N,EH_LIVER.N,EH_OBESITY.N,EH_OTHNEURO.N,EH_PERIVASC.N,EH_PULMCIRC.N,EH_RENAL.N,EH_RHEUMART.N,EH_UNCDIAB.N,EH_VALVDIS.N,EH_WEIGHTLS.N) d1=Reduce(function(...) merge(...,all=T,by="PatientICN"),RVDx) for (i in seq(11,86)){d1[which(is.na(d1[,..i])),i]=0} #make sure columns 3-16 are the counts per condition columns. replace NA with 0 #for (i in seq(9,84)){d1[which(d1[,..i] >=1),i]=1} if (count==1){dDx1=d1} if (count==2){dDx2=d1} if (count==3){dDx3=d1} if (count==4){dDx4=d1} if (count==5){dDx5=d1} if (count==6){dDx6=d1} if (count==7){dDx7=d1} if (count==8){dDx8=d1} } #fit=glm(OC ~ SuicideIdeation+alcdx_poss+Nicdx_poss+AmphetamineUseDisorder+Cannabis+COCNdx+DrugInducedSymptoms+OpioidOverdose+OUD+SedativeUseDisorder+SAE_sed+SUD_CatchAll+SAE_Acet+SAE_OtherDrug +MDD+Schiz+Sleep+BIPOLAR+PTSD+MH_CatchAll+ODEPRdx_poss+OtherPsychosis+Thyroid+reach_othanxdis+reach_persond+TBI_Dx+Concuss+Headache+SAE_Falls+SAE_OtherAccident+SAE_Vehicle+Amputation+SpinalCordInj+Backpain+EH_PARALYSIS+DEMENTIA+Parkinsons+DeliriumTremens+Huntington+EH_Lymphoma+EH_METCANCR+EH_NMETTUMR +HIV+Tourette+Osteoporosis+Fibromyalgia+Hearing+Visual+Neuro+reach_chronic+GastricUlcer+SleepApnea+Homeless+AMI+CAD+ChestPain_Dx+EH_ARRHYTH+EH_BLANEMIA+EH_CHRNPULM+EH_COAG+EH_COMDIAB+EH_DefANEMIA+EH_ELECTRLYTE+EH_HEART+EH_HYPERTENS+EH_LIVER+EH_OBESITY+EH_OTHNEURO+EH_PERIVALV+EH_PULMCIRC+EH_RENAL+EH_RHEUMART+EH_UNCDIAB+EH_VALVDIS+EH_WEIGHTLS, data=d1, family='binomial') #medications #Reachvet2 variables #Paton_Antipsychotic_Rx Antipsychotic_Rx_6 Antipsychotic_Rx_12 Paton_MoodStabilizer_GE3_Rx MoodStabilizer_GE3_Rx_6 MoodStabilizer_GE3_Rx_12 #Paton_Antidepressant_Rx Antidepressant_Rx_12 Paton_Stimulant_Rx Stimulant_Rx_12 #Paton_Anxiolytics_Rx Anxiolytics_Rx_12 Paton_Barbiturate_Rx Barbiturate_Rx_12 #Paton_Insomnia_Rx Insomnia_Rx_12 Paton_SUD_Rx SUD_Rx_12 Paton_PainAdjAnticonvulsant_Rx #PainAdjAnticonvulsant_Rx_12 #Paton_Combined_Analgesics_rx Combined_Analgesics_rx_12 #Paton_Reach_statin_RX Reach_statin_RX_12 #Paton_OpioidForPain_Rx OpioidForPain_Rx_6 OpioidForPain_Rx_12 #MaxMPR_Antipsychotic_Rx_LOW MaxMPR_MoodStabilizer_GE3_Rx_LOW MaxMPR_Antidepressant_Rx_LOW MaxMPR_Stimulant_Rx_LOW MaxMPR_Anxiolytics_Rx_LOW MaxMPR_Barbiturate_Rx_LOW #MaxMPR_Insomnia_Rx_LOW MaxMPR_SUD_Rx_LOW MaxMPR_PainAdjAnticonvulsant_Rx_LOW MaxMPR_Combined_Analgesics_rx_LOW MaxMPR_Reach_statin_RX_LOW MaxMPR_OpioidForPain_Rx_LOW CN=readRDS(paste(c(C2020dir,"meds/Rx.CN.",cohort,".RDS"),collapse="")) CN[,"LocalDrugNameWithDose":=NULL] CN[,"DFactor":=factor(paste(VAClassification,DrugNameWithoutDose,sep='.'))] CN$VAClassification=factor(CN$VAClassification) CV=readRDS(paste(c(C2020dir,"meds/Rx.CV.",cohort,".RDS"),collapse="")) CV[,"LocalDrugNameWithDose":=NULL] CV[,"DFactor":=factor(paste(VAClassification,DrugNameWithoutDose,sep='.'))] CV$VAClassification=factor(CV$VAClassification) A=readRDS(paste(c(C2020dir,"meds/Rx.A.",cohort,".RDS"),collapse="")) A[,"LocalDrugNameWithDose":=NULL] A[,"DFactor":=factor(paste(VAClassification,DrugNameWithoutDose,sep='.'))] A$VAClassification=factor(A$VAClassification) HS=readRDS(paste(c(C2020dir,"meds/Rx.HS.",cohort,".RDS"),collapse="")) HS[,"LocalDrugNameWithDose":=NULL] HS[,"DFactor":=factor(paste(VAClassification,DrugNameWithoutDose,sep='.'))] HS$VAClassification=factor(HS$VAClassification) oth=readRDS(paste(c(C2020dir,"meds/Rx.Oth.",cohort,".RDS"),collapse="")) oth[,"LocalDrugNameWithDose":=NULL] oth[,"DFactor":=factor(paste(VAClassification,DrugNameWithoutDose,sep='.'))] oth$VAClassification=factor(oth$VAClassification) Rx=rbind(CN,CV,A,HS,oth) mRx=merge(OC,Rx,all.x=T,by="PatientICN") Rx=mRx Insomnia_Rx=Rx[grep("SUVOREXANT|TASIMELTEON|RAMELTEON",DrugNameWithoutDose),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose)] Barbiturate_Rx=Rx[grep("AMOBARBITAL|butaBARBITAL|MEPHOBARBITAL|PENTOBARB|PHENOBARB|SECOBARB|PRIMIDONE|BUTALBITAL",DrugNameWithoutDose),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose)] Combined_Analgesics_Rx=Rx[grep("CN10[1235]|MS10[12]|DE650|DE[27]00",VAClassification),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose)] AlcoholPharmacotherapy_notop_Rx=Rx[grep("NALTREXONE|DISULFIRAM|ACAMPROSATE",DrugNameWithoutDose),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose)] MoodStabilizer_GE3_Rx=Rx[grep("CARBAMAZEPINE|DIVALPROEX|FELBAMATE|GABAPENTIN|LAMOTRIGINE|OXCARBAZEPINE|TOPIRAMATE|VALPROIC",DrugNameWithoutDose),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose)] PainAdjAnticonvulsant_Rx=Rx[grep("CARBAMAZEPINE|LEVETIRACETAM|ZONISAMIDE|TIAGABINE|GABAPENTIN|PREGABALIN|OXCARBAZEPINE|TOPIRAMATE|VALPROIC",DrugNameWithoutDose),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose)] Antidepressant_Rx=Rx[grep("CHLORPROMAZINE|FLUPHENAZINE|PERPHENAZINE|THIORIDAZINE|THIOTHIXENE|TRIFLUOPERAZINE|ARIPIPRAZOLE|ASENAPINE|BREXPIPRAZOLE|CARIPRAZINE|CLOZAPINE|HALOPERIDOL|ILOPERIDONE|LOXAPINE|LURASIDONE|OLANZAPINE|PALIPERIDONE|PIMAVANSERIN|QUETIAPINE|RISPERIDONE|ZIPRASIDONE|LITHIUM|ACETOPHENAZINE|ARIPIPRAZOLE|ASENAPINE|BREXPIPRAZOLE|CARIPRAZINE|CHLORPROTHIXENE|DROPERIDOL|FLUPHENAZINE|HALOPERIDOL|ILOPERIDONE|LURASIDONE|MESORIDAZINE|MOLINDONE|PIMOZIDE|PROCHLORPERAZINE|PROMETHAZINE",DrugNameWithoutDose),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose)] Stimulant_Rx=Rx[grep("AMPHETAMINE|METHYLPHENIDATE|MODAFINIL",DrugNameWithoutDose),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose)] Thiamine_Rx=Rx[grep("THIAMINE",DrugNameWithoutDose),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose)] Anxiolytics_Rx=Rx[grep("ALPRAZOLAM|CHLORDIAZEPOXIDE|CLONAZEPAM|CLORAZEPATE|DIAZEPAM|FLURAZEPAM|LORAZEPAM|OXAZEPAM|TEMAZEPAM|TRIAZOLAM|HALAZEPAM|ESTAZOLAM|ZOLPIDEM|QUAZEPAM|ZALEPLON|ESZOPICLONE|BUSPIRONE|CHLORAL HYDRATE|PRAZEPAM|CHLORDIAZEPOXIDE",DrugNameWithoutDose),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose)] OAtmp=Rx[grep("BUPRENORPHINE",DrugNameWithoutDose),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose,NationalDrugNameWithDose)] OpioidAgonist_Rx=OAtmp[!grep("PATCH|INJ",NationalDrugNameWithDose),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose)] Ninj1=Rx[grep("NALTREXONE",NationalDrugNameWithDose),] Ninj2=Ninj1[grep("INJ",NationalDrugNameWithDose),] NaltrexoneINJ_Rx=Ninj2[!grep("METHYL",NationalDrugNameWithDose),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose)] Op1=Rx[c(grep("CN101",VAClassification),grep("TRAMADOL",DrugNameWithoutDose))] OpioidForPain_Rx=Op1[!grep("DOVERIN|BUPRENOR|NALOX|IPECAC|BUPRE",NationalDrugNameWithDose),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose)] Reach_Statin_Rx=Rx[grep("ATORVASTATIN|FLUVASTATIN|LOVASTATIN|PITAVASTATIN|PRAVASTATIN|ROSUVASTATIN|SIMVASTATIN",DrugNameWithoutDose),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose)] CVD_Rx=Rx[grep("^CV",VAClassification),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose)] Antimicrobial_Rx=Rx[grep("^AM",VAClassification),.(PatientICN,FillDateTime,Tpr,DrugNameWithoutDose)] #CN60[019] = ACETOPHENAZINE|ARIPIPRAZOLE|ASENAPINE|BREXPIPRAZOLE|CARIPRAZINE|CHLORPROTHIXENE|DROPERIDOL|FLUPHENAZINE|HALOPERIDOL|ILOPERIDONE|LURASIDONE|MESORIDAZINE|MOLINDONE|PIMOZIDE|PROCHLORPERAZINE|PROMETHAZINE #cn302 MINUS MIDAZOLAM = ALPRAZOLAM CHLORDIAZEPOXIDE CLONAZEPAM CLORAZEPATE DIAZEPAM FLURAZEPAM LORAZEPAM OXAZEPAM TEMAZEPAM TRIAZOLAM #COUNT INSTANCES for (count in c(1,2,3,4,5,6,7,8)){ tWin=window[count];tGap=gindow[count] Insomnia_Rx.N=Insomnia_Rx[FillDateTime < Tpr-tGap & FillDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Insomnia_Rx.N)[2]="Insomnia_Rx" Barbiturate_Rx.N=Barbiturate_Rx[FillDateTime < Tpr-tGap & FillDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Barbiturate_Rx.N)[2]="Barbiturate_Rx" Combined_Analgesics_Rx.N=Combined_Analgesics_Rx[FillDateTime < Tpr-tGap & FillDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Combined_Analgesics_Rx.N)[2]="Combined_Analgesics_Rx" AlcoholPharmacotherapy_notop_Rx.N=AlcoholPharmacotherapy_notop_Rx[FillDateTime < Tpr-tGap & FillDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(AlcoholPharmacotherapy_notop_Rx.N)[2]="AlcoholPharmacotherapy_notop_Rx" MoodStabilizer_GE3_Rx.N=MoodStabilizer_GE3_Rx[FillDateTime < Tpr-tGap & FillDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(MoodStabilizer_GE3_Rx.N)[2]="MoodStabilizer_GE3_Rx" PainAdjAnticonvulsant_Rx.N=PainAdjAnticonvulsant_Rx[FillDateTime < Tpr-tGap & FillDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(PainAdjAnticonvulsant_Rx.N)[2]="PainAdjAnticonvulsant_Rx" Antidepressant_Rx.N=Antidepressant_Rx[FillDateTime < Tpr-tGap & FillDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Antidepressant_Rx.N)[2]="Antidepressant_Rx" Stimulant_Rx.N=Stimulant_Rx[FillDateTime < Tpr-tGap & FillDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Stimulant_Rx.N)[2]="Stimulant_Rx" Thiamine_Rx.N=Thiamine_Rx[FillDateTime < Tpr-tGap & FillDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Thiamine_Rx.N)[2]="Thiamine_Rx" Anxiolytics_Rx.N=Anxiolytics_Rx[FillDateTime < Tpr-tGap & FillDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Anxiolytics_Rx.N)[2]="Anxiolytics_Rx" OpioidAgonist_Rx.N=OpioidAgonist_Rx[FillDateTime < Tpr-tGap & FillDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(OpioidAgonist_Rx.N)[2]="OpioidAgonist_Rx" NaltrexoneINJ_Rx.N=NaltrexoneINJ_Rx[FillDateTime < Tpr-tGap & FillDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(NaltrexoneINJ_Rx.N)[2]="NaltrexoneINJ_Rx" OpioidForPain_Rx.N=OpioidForPain_Rx[FillDateTime < Tpr-tGap & FillDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(OpioidForPain_Rx.N)[2]="OpioidForPain_Rx" Reach_Statin_Rx.N=Reach_Statin_Rx[FillDateTime < Tpr-tGap & FillDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Reach_Statin_Rx.N)[2]="Reach_Statin_Rx" CVD_Rx.N=CVD_Rx[FillDateTime < Tpr-tGap & FillDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(CVD_Rx.N)[2]="CVD_Rx" Antimicrobial_Rx.N=Antimicrobial_Rx[FillDateTime < Tpr-tGap & FillDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Antimicrobial_Rx.N)[2]="Antimicrobial_Rx" RVRx=list(OC,Insomnia_Rx.N,Barbiturate_Rx.N,Combined_Analgesics_Rx.N,AlcoholPharmacotherapy_notop_Rx.N,MoodStabilizer_GE3_Rx.N,PainAdjAnticonvulsant_Rx.N,Antidepressant_Rx.N,Stimulant_Rx.N,Thiamine_Rx.N,Anxiolytics_Rx.N,OpioidAgonist_Rx.N,NaltrexoneINJ_Rx.N,OpioidForPain_Rx.N,Reach_Statin_Rx.N,CVD_Rx.N,Antimicrobial_Rx.N) dRx=Reduce(function(...) merge(...,all=T,by="PatientICN"),RVRx) for (i in seq(11,26)){dRx[which(is.na(dRx[,..i])),i]=0} #make sure columns 3-16 are the counts per condition columns. replace NA with 0 #for (i in seq(9,21)){dRx[which(dRx[,..i] >=1),i]=1} if (count==1){dRx1=dRx} if (count==2){dRx2=dRx} if (count==3){dRx3=dRx} if (count==4){dRx4=dRx} if (count==5){dRx5=dRx} if (count==6){dRx6=dRx} if (count==7){dRx7=dRx} if (count==8){dRx8=dRx}} #saveRDS(dRx,file=paste(c("Rx.",cohort,".RDS"),collapse="")) #fit=glm(OC ~ Insomnia_Rx+Barbiturate_Rx+Combined_Analgesics_Rx+AlcoholPharmacotherapy_notop_Rx+MoodStabilizer_GE3_Rx+PainAdjAnticonvulsant_Rx+Antidepressant_Rx+Stimulant_Rx+Anxiolytics_Rx+OpioidAgonist_Rx+NaltrexoneINJ_Rx+OpioidForPain_Rx+Reach_Statin_Rx, data=dRx, family='binomial') #CPT iCPT=readRDS(paste(c(C2020dir,"inpat/iCPT.",cohort,".RDS"),collapse="")) colnames(iCPT)=c("PatientICN","VisitSID","CPTCode","CPTDateTime") oCPT.EM=readRDS(paste(c(C2020dir,"outpat/oCPT.EM.",cohort,".RDS"),collapse="")) oCPT.EM$sta3n=NULL;oCPT.EM$VisitDateTime=NULL colnames(oCPT.EM)=c("PatientICN","VisitSID","CPTCode","CPTDateTime") oCPT.Med=readRDS(paste(c(C2020dir,"outpat/oCPT.Med.",cohort,".RDS"),collapse="")) oCPT.Med$sta3n=NULL;oCPT.Med$VisitDateTime=NULL colnames(oCPT.Med)=c("PatientICN","VisitSID","CPTCode","CPTDateTime") oCPT.Oth=readRDS(paste(c(C2020dir,"outpat/oCPT.Other.",cohort,".RDS"),collapse="")) oCPT.Oth$sta3n=NULL;oCPT.Oth$VisitDateTime=NULL colnames(oCPT.Oth)=c("PatientICN","VisitSID","CPTCode","CPTDateTime") CPT=rbind(iCPT,oCPT.EM,oCPT.Med,oCPT.Oth) mCPT=merge(OC,CPT,all.x=T,by="PatientICN") CPT=mCPT rm(oCPT.EM); rm(oCPT.Med); rm(oCPT.Oth) GroupTx=CPT[grep("9085[37]|96153",CPTCode),.(PatientICN,CPTDateTime,Tpr,CPTCode)] MedMgmnt=CPT[grep("90862|H0031|96152",CPTCode),.(PatientICN,CPTDateTime,Tpr,CPTCode)] PsyTx25min=CPT[grep("9080[45]",CPTCode),.(PatientICN,CPTDateTime,Tpr,CPTCode)] PsyTx45min=CPT[grep("9080[67]",CPTCode),.(PatientICN,CPTDateTime,Tpr,CPTCode)] PsyTx75min=CPT[grep("9080[89]",CPTCode),.(PatientICN,CPTDateTime,Tpr,CPTCode)] PsyInterview.Test=CPT[grep("90801|96100|9611[57]|96150",CPTCode),.(PatientICN,CPTDateTime,Tpr,CPTCode)] FamilyPsy=CPT[grep("90847|90887|9084[69]",CPTCode),.(PatientICN,CPTDateTime,Tpr,CPTCode)] AlcTx=CPT[grep("H000[123456]|H00[12]5|H0024",CPTCode),.(PatientICN,CPTDateTime,Tpr,CPTCode)] Breath.EtOH=CPT[grep("82075",CPTCode),.(PatientICN,CPTDateTime,Tpr,CPTCode)] MethadoneClinic=CPT[grep("H0020",CPTCode),.(PatientICN,CPTDateTime,Tpr,CPTCode)] PsyReport=CPT[grep("9088[59]",CPTCode),.(PatientICN,CPTDateTime,Tpr,CPTCode)] EnvironManipulation=CPT[grep("90882",CPTCode),.(PatientICN,CPTDateTime,Tpr,CPTCode)] InpatPsy20min=CPT[grep("9081[67]",CPTCode),.(PatientICN,CPTDateTime,Tpr,CPTCode)] InpatPsy45min=CPT[grep("90818",CPTCode),.(PatientICN,CPTDateTime,Tpr,CPTCode)] InteractivePsy=CPT[grep("9081[12]",CPTCode),.(PatientICN,CPTDateTime,Tpr,CPTCode)] for (count in c(1,2,3,4,5,6,7,8)){ tWin=window[count];tGap=gindow[count] GroupTx.N=GroupTx[CPTDateTime < Tpr-tGap & CPTDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(GroupTx.N)[2]="GroupTx" MedMgmnt.N=MedMgmnt[CPTDateTime < Tpr-tGap & CPTDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(MedMgmnt.N)[2]="MedMgmnt" PsyTx25min.N=PsyTx25min[CPTDateTime < Tpr-tGap & CPTDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(PsyTx25min.N)[2]="PsyTx25min" PsyTx45min.N=PsyTx45min[CPTDateTime < Tpr-tGap & CPTDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(PsyTx45min.N)[2]="PsyTx45min" PsyTx75min.N=PsyTx75min[CPTDateTime < Tpr-tGap & CPTDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(PsyTx75min.N)[2]="PsyTx75min" PsyInterview.Test.N=PsyInterview.Test[CPTDateTime < Tpr-tGap & CPTDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(PsyInterview.Test.N)[2]="PsyInterview.Test" FamilyPsy.N=FamilyPsy[CPTDateTime < Tpr-tGap & CPTDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(FamilyPsy.N)[2]="FamilyPsy" AlcTx.N=AlcTx[CPTDateTime < Tpr-tGap & CPTDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(AlcTx.N)[2]="AlcTx" Breath.EtOH.N=Breath.EtOH[CPTDateTime < Tpr-tGap & CPTDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(Breath.EtOH.N)[2]="Breath.EtOH" MethadoneClinic.N=MethadoneClinic[CPTDateTime < Tpr-tGap & CPTDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(MethadoneClinic.N)[2]="MethadoneClinic" PsyReport.N=PsyReport[CPTDateTime < Tpr-tGap & CPTDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(PsyReport.N)[2]="PsyReport" EnvironManipulation.N=EnvironManipulation[CPTDateTime < Tpr-tGap & CPTDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(EnvironManipulation.N)[2]="EnvironManipulation" InpatPsy20min.N=InpatPsy20min[CPTDateTime < Tpr-tGap & CPTDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(InpatPsy20min.N)[2]="InpatPsy20min" InpatPsy45min.N=InpatPsy45min[CPTDateTime < Tpr-tGap & CPTDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(InpatPsy45min.N)[2]="InpatPsy45min" InteractivePsy.N=InteractivePsy[CPTDateTime < Tpr-tGap & CPTDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(InteractivePsy.N)[2]="InteractivePsy" lCPT=list(OC,GroupTx.N,MedMgmnt.N,PsyTx25min.N,PsyTx45min.N,PsyTx75min.N,PsyInterview.Test.N,FamilyPsy.N,AlcTx.N,Breath.EtOH.N,MethadoneClinic.N,PsyReport.N,EnvironManipulation.N,InpatPsy20min.N,InteractivePsy.N) dCPT=Reduce(function(...) merge(...,all=T,by="PatientICN"),lCPT) for (i in seq(11,24)){dCPT[which(is.na(dCPT[,..i])),i]=0} #make sure columns 3-16 are the counts per condition columns. replace NA with 0 #for (i in seq(9,21)){dRx[which(dRx[,..i] >=1),i]=1} if (count==1){dCPT1=dCPT} if (count==2){dCPT2=dCPT} if (count==3){dCPT3=dCPT} if (count==4){dCPT4=dCPT} if (count==5){dCPT5=dCPT} if (count==6){dCPT6=dCPT} if (count==7){dCPT7=dCPT} if (count==8){dCPT8=dCPT} } #appointments #Reachvet2 variabl #HighMORate_91days HighMORate_182days HighMORate_365days HighMORate_730days #TotalNumOfAppt_91days_0 TotalNumOfAppt_91days_1up TotalNumOfAppt_91days_3up TotalNumOfAppt_91days_6up TotalNumOfAppt_91days_9up #TotalNumOfAppt_182days_0 TotalNumOfAppt_182days_1up TotalNumOfAppt_182days_6up TotalNumOfAppt_182days_12up TotalNumOfAppt_182days_18up #TotalNumOfAppt_365days_0 TotalNumOfAppt_365days_1up TotalNumOfAppt_365days_12up TotalNumOfAppt_365days_24up TotalNumOfAppt_365days_36up #TotalNumOfAppt_730days_0 TotalNumOfAppt_730days_1up TotalNumOfAppt_730days_24up TotalNumOfAppt_730days_48up TotalNumOfAppt_730days_96up # Mental health #Reachvet2 variables #PHQ_RawScore PHQ_TimeSinceLastSurvey #PHQ_NumberOfSurveyCount PHQ_Minimal_Depression PHQ_Mild_Depression PHQ_Moderate_Depression PHQ_Moderately_Severe_Depression PHQ_Severe_Depression #PHQ_MostRecent_Question9_Score1up PHQ_MostRecent_Question9_Score2up PHQ_MostRecent_Question9_Score3up PHQ_Hist_Question9_Score1up PHQ_Hist_Question9_Score2up PHQ_Hist_Question9_Score3up PHQ_MDD PHQ_ODD MH=readRDS(paste(c(C2020dir,"survey/SResult.",cohort,".RDS"),collapse="")) MHQ=readRDS(paste(c(C2020dir,"survey/SQuestions.",cohort,".RDS"),collapse="")) MH2=merge(OC,MH[SurveyName=="PHQ-2",],all.x=T,by="PatientICN") ;setkeyv(MH2,c("PatientICN","SurveyGivenDateTime")) MHQ2=merge(OC,MHQ[SurveyName=="PHQ-2",],all.x=T,by="PatientICN") ;setkeyv(MHQ2,c("PatientICN","SurveyGivenDateTime")) MH9=merge(OC,MH[SurveyName=="PHQ9",],all.x=T,by="PatientICN") ;setkeyv(MH9,c("PatientICN","SurveyGivenDateTime")) MHQ9=merge(OC,MHQ[SurveyName=="PHQ9",],all.x=T,by="PatientICN") ;setkeyv(MHQ9,c("PatientICN","SurveyGivenDateTime")) #COUNT INSTANCES for (count in c(1,2,3,4,5,6,7,8)){ tWin=window[count];tGap=gindow[count] MH2k=MH2[SurveyGivenDateTime < Tpr-tGap & SurveyGivenDateTime > Tpr-tWin-tGap,] MHQ2k=MHQ2[SurveyGivenDateTime < Tpr-tGap & SurveyGivenDateTime > Tpr-tWin-tGap,] MH9k=MH9[SurveyGivenDateTime < Tpr-tGap & SurveyGivenDateTime > Tpr-tWin-tGap,] MHQ9k=MHQ9[SurveyGivenDateTime < Tpr-tGap & SurveyGivenDateTime > Tpr-tWin-tGap,] PHQ2_Count=MH2k[,.N,by="PatientICN"];colnames(PHQ2_Count)[2]="PHQ2_Count" PHQ9_Count=MH9k[,.N,by="PatientICN"];colnames(PHQ9_Count)[2]="PHQ9_Count" PHQ2_Raw.max=MH2k[,max(RawScore),by="PatientICN"];colnames(PHQ2_Raw.max)[2]="PHQ2_Raw.max" PHQ9_Raw.max=MH9k[,max(RawScore),by="PatientICN"];colnames(PHQ9_Raw.max)[2]="PHQ9_Raw.max" PHQ2_Raw.recent=MH2k[,last(RawScore),by="PatientICN"];colnames(PHQ2_Raw.recent)[2]="PHQ2_Raw.recent" PHQ9_Raw.recent=MH9k[,last(RawScore),by="PatientICN"];colnames(PHQ9_Raw.recent)[2]="PHQ9_Raw.recent" Pmax=MHQ9k[grep("dead or of hurting",SurveyQuestionText),max(SurveyChoiceText),by="PatientICN"] Plast=MHQ9k[grep("dead or of hurting",SurveyQuestionText),last(SurveyChoiceText),by="PatientICN"] PHQ9_Q9.max1=data.table(c(0,Pmax[V1=="Not at all",PatientICN]),1) ;colnames(PHQ9_Q9.max1)=c("PatientICN","PHQ9_Q9.max1") PHQ9_Q9.recent1=data.table(c(0,Plast[V1=="Not at all",PatientICN]),1) ;colnames(PHQ9_Q9.recent1)=c("PatientICN","PHQ9_Q9.recent1") PHQ9_Q9.max2= data.table(c(0,Pmax[V1=="Several days",PatientICN]),1) ;colnames(PHQ9_Q9.max2)=c("PatientICN","PHQ9_Q9.max2") PHQ9_Q9.recent2= data.table(c(0,Plast[V1=="Several days",PatientICN]),1) ;colnames(PHQ9_Q9.recent2)=c("PatientICN","PHQ9_Q9.recent2") PHQ9_Q9.max3= data.table(c(0,Pmax[V1=="More than half the days",PatientICN]),1) ;colnames(PHQ9_Q9.max3)=c("PatientICN","PHQ9_Q9.max3") PHQ9_Q9.recent3= data.table(c(0,Plast[V1=="More than half the days",PatientICN]),1) ;colnames(PHQ9_Q9.recent3)=c("PatientICN","PHQ9_Q9.recent3") PHQ9_Q9.max4= data.table(c(0,Pmax[V1=="Nearly every day",PatientICN]),1) ;colnames(PHQ9_Q9.max4)=c("PatientICN","PHQ9_Q9.max4") PHQ9_Q9.recent4= data.table(c(0,Plast[V1=="Nearly every day",PatientICN]),1) ;colnames(PHQ9_Q9.recent4)=c("PatientICN","PHQ9_Q9.recent4") #Patient Record flag RVphq=list(OC,PHQ2_Count,PHQ9_Count,PHQ2_Raw.max,PHQ9_Raw.max,PHQ2_Raw.recent,PHQ9_Raw.recent,PHQ9_Q9.max1,PHQ9_Q9.max2,PHQ9_Q9.max3,PHQ9_Q9.max4,PHQ9_Q9.recent1,PHQ9_Q9.recent2,PHQ9_Q9.recent3,PHQ9_Q9.recent4) dphq=Reduce(function(...) merge(...,all.x=T,by="PatientICN"),RVphq) for (i in seq(11,24)){dphq[which(is.na(dphq[,..i])),i]=0} #make sure columns 3-16 are the counts per condition columns. replace NA with 0 if (count==1){dPHQ1=dphq} if (count==2){dPHQ2=dphq} if (count==3){dPHQ3=dphq} if (count==4){dPHQ4=dphq} if (count==5){dPHQ5=dphq} if (count==6){dPHQ6=dphq} if (count==7){dPHQ7=dphq} if (count==8){dPHQ8=dphq} } #saveRDS(dphq,file=paste(c("PHQ.",cohort,".RDS"),collapse="")) #fit=glm(OC ~ PHQ2_Count+PHQ9_Count+PHQ2_Raw.max+PHQ9_Raw.max+PHQ2_Raw.recent+PHQ9_Raw.recent+PHQ9_Q9.max1+PHQ9_Q9.max2+PHQ9_Q9.max3+PHQ9_Q9.max4+PHQ9_Q9.recent1+PHQ9_Q9.recent2+PHQ9_Q9.recent3+PHQ9_Q9.recent4, data=dphq, family='binomial') #summary(fit) #utilization #PoC_ED_Stop PoC_UC_Stop PoC_TobaccoCessation_Stop MHOC_HBPC_Stop MHOC_GeneralMentalHealth_Stop MHOC_PTSD_Stop MHOC_SUD_Stop # , MHOC_TSES_Stop MHOC_PRRC_Stop MHOC_MHICM_Stop MHOC_PCMHI_Stop MHOC_Residential_Stop MHOC_Homeless_Stop MHOC_OtherMentalHealth_Stop oFac=readRDS(paste(c(C2020dir,"outpat/oFac.",cohort,".RDS"),collapse="")) oFac=oFac[is.na(match(PatientICN,pts))==F,] moFac=merge(OC,oFac,all.y=T,by="PatientICN") oFac=moFac[,c("PatientICN","VisitDateTime","StopCode","StopCode.1","VisitDateTime","Tpr")] PoC_ED=oFac[is.na(match(StopCode,c(130,297)))==F | is.na(match(StopCode.1,c(130,297)))==F,.(PatientICN,VisitDateTime,Tpr)] PoC_UC=oFac[is.na(match(StopCode,c(131)))==F | is.na(match(StopCode.1,c(131)))==F,.(PatientICN,VisitDateTime,Tpr)] PoC_TobaccoCessation=oFac[is.na(match(StopCode,c(566,707,708)))==F | is.na(match(StopCode.1,c(566,707,708)))==F,.(PatientICN,VisitDateTime,Tpr)] MHOC_HBPC=oFac[is.na(match(StopCode,c(156,157)))==F | is.na(match(StopCode.1,c(156,157)))==F,.(PatientICN,VisitDateTime,Tpr)] MHOC_GeneralMentalHealth=HBPC=oFac[is.na(match(StopCode,c(502,509,510,550,557,558)))==F | is.na(match(StopCode.1,c(502,509,510,550,557,558)))==F,.(PatientICN,VisitDateTime,Tpr)] MHOC_PTSD=oFac[is.na(match(StopCode,c(516,562,519,525,540,542,561,580,581)))==F | is.na(match(StopCode.1,c(516,562,519,525,540,542,561,580,581)))==F,.(PatientICN,VisitDateTime,Tpr)] MHOC_SUD=oFac[is.na(match(StopCode,c(513,514,523,545,547,548,560)))==F | is.na(match(StopCode.1,c(513,514,523,545,547,548,560)))==F,.(PatientICN,VisitDateTime,Tpr)] MHOC_TSES=oFac[is.na(match(StopCode,c(535,568,573,574,575,536)))==F | is.na(match(StopCode.1,c(535,568,573,574,575,536)))==F,.(PatientICN,VisitDateTime,Tpr)] MHOC_PRRC=oFac[is.na(match(StopCode,c(582,583,584)))==F | is.na(match(StopCode.1,c(582,583,584)))==F,.(PatientICN,VisitDateTime,Tpr)] MHOC_MHICM=oFac[is.na(match(StopCode,c(546,552,567)))==F | is.na(match(StopCode.1,c(546,552,567)))==F,.(PatientICN,VisitDateTime,Tpr)] MHOC_PCMHI=oFac[is.na(match(StopCode,c(534,539)))==F | is.na(match(StopCode.1,c(534,539)))==F,.(PatientICN,VisitDateTime,Tpr)] MHOC_Residential=oFac[is.na(match(StopCode,c(586,587,588,593,594,595,596,597,598,599)))==F | is.na(match(StopCode.1,c(586,587,588,593,594,595,596,597,598,599)))==F,.(PatientICN,VisitDateTime,Tpr)] MHOC_Homeless=oFac[is.na(match(StopCode,c(504,507,508,511,522,528,529,530,555,556,590,591,592)))==F | is.na(match(StopCode.1,c(504,507,508,511,522,528,529,530,555,556,590,591,592)))==F,.(PatientICN,VisitDateTime,Tpr)] MHOC_OtherMH=oFac[is.na(match(StopCode,c(292,503,505,506,512,524,527,531,532,533,537,538,553,554,559,563,564,565,566,571,572,576,577,578,579,589,713)))==F | is.na(match(StopCode.1,c(292,503,505,506,512,524,527,531,532,533,537,538,553,554,559,563,564,565,566,571,572,576,577,578,579,589,713)))==F,.(PatientICN,VisitDateTime,Tpr)] #COUNT INSTANCES for (count in c(1,2,3,4,5,6,7,8)){ tWin=window[count];tGap=gindow[count] PoC_ED.N=PoC_ED[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(PoC_ED.N)[2]="PoC_ED" PoC_UC.N=PoC_UC[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(PoC_UC.N)[2]="PoC_UC" PoC_TobaccoCessation.N=PoC_TobaccoCessation[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(PoC_TobaccoCessation.N)[2]="PoC_TobaccoCessation" MHOC_HBPC.N=MHOC_HBPC[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(MHOC_HBPC.N)[2]="MHOC_HBPC" MHOC_GeneralMentalHealth.N=MHOC_GeneralMentalHealth[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(MHOC_GeneralMentalHealth.N)[2]="MHOC_GeneralMentalHealth" MHOC_PTSD.N=MHOC_PTSD[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(MHOC_PTSD.N)[2]="MHOC_PTSD" MHOC_SUD.N=MHOC_SUD[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(MHOC_SUD.N)[2]="MHOC_SUD" MHOC_TSES.N=MHOC_TSES[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(MHOC_TSES.N)[2]="MHOC_TSES" MHOC_PRRC.N=MHOC_PRRC[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(MHOC_PRRC.N)[2]="MHOC_PRRC" MHOC_MHICM.N=MHOC_MHICM[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(MHOC_MHICM.N)[2]="MHOC_MHICM" MHOC_PCMHI.N=MHOC_PCMHI[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(MHOC_PCMHI.N)[2]="MHOC_PCMHI" MHOC_Residential.N=MHOC_Residential[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(MHOC_Residential.N)[2]="MHOC_Residential" MHOC_Homeless.N=MHOC_Homeless[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(MHOC_Homeless.N)[2]="MHOC_Homeless" MHOC_OtherMH.N=MHOC_OtherMH[VisitDateTime < Tpr-tGap & VisitDateTime > Tpr-tWin-tGap,.N,by="PatientICN"];colnames(MHOC_OtherMH.N)[2]="MHOC_OtherMH" Util=list(OC,PoC_ED.N,PoC_UC.N,PoC_TobaccoCessation.N,MHOC_HBPC.N,MHOC_GeneralMentalHealth.N,MHOC_PTSD.N,MHOC_SUD.N,MHOC_TSES.N,MHOC_PRRC.N,MHOC_MHICM.N,MHOC_PCMHI.N,MHOC_Residential.N,MHOC_Homeless.N,MHOC_OtherMH.N) d2=Reduce(function(...) merge(...,all=T,by="PatientICN"),Util) for (i in seq(11,24)){d2[which(is.na(d2[,..i])),i]=0} #make sure columns 3-16 are the counts per condition columns. replace NA with 0 #for (i in seq(9,22)){d2[which(d2[,..i] >=1),i]=1} if (count==1){dUtil1=d2} if (count==2){dUtil2=d2} if (count==3){dUtil3=d2} if (count==4){dUtil4=d2} if (count==5){dUtil5=d2} if (count==6){dUtil6=d2} if (count==7){dUtil7=d2} if (count==8){dUtil8=d2} } # VitalSigns Vbp=readRDS(paste(c(C2020dir,"vital/Vital.bp.",cohort,".RDS"),collapse="")) mVbp=merge(OC,Vbp,all.x=T,by="PatientICN") Vbp=mVbp Vo=readRDS(paste(c(C2020dir,"vital/Vital.other.",cohort,".RDS"),collapse="")) mVo=merge(OC,Vo,all.x=T,by="PatientICN") Vo=mVo for (count in c(1,2,3,4,5,6,7,8)){ tWin=window[count];tGap=gindow[count] Systolic=Vbp[VitalSignTakenDateTime < Tpr-tGap & VitalSignTakenDateTime > Tpr-tWin-tGap,median(Systolic),by="PatientICN"];colnames(Systolic)[2]="Systolic" Diastolic=Vbp[VitalSignTakenDateTime < Tpr-tGap & VitalSignTakenDateTime > Tpr-tWin-tGap,median(Diastolic),by="PatientICN"];colnames(Diastolic)[2]="Diastolic" Weight=Vo[(VitalSignTakenDateTime < Tpr-tGap & VitalSignTakenDateTime > Tpr-tWin-tGap) & (VitalTypeAbbreviation=="WT"), median(VitalResultNumeric),by="PatientICN"];colnames(Weight)[2]="Weight" Pulse=Vo[(VitalSignTakenDateTime < Tpr-tGap & VitalSignTakenDateTime > Tpr-tWin-tGap) & (VitalTypeAbbreviation=="P"), median(VitalResultNumeric),by="PatientICN"];colnames(Pulse)[2]="Pulse" Respiration=Vo[(VitalSignTakenDateTime < Tpr-tGap & VitalSignTakenDateTime > Tpr-tWin-tGap) & (VitalTypeAbbreviation=="R"), median(VitalResultNumeric),by="PatientICN"];colnames(Respiration)[2]="Respiration" PO2=Vo[(VitalSignTakenDateTime < Tpr-tGap & VitalSignTakenDateTime > Tpr-tWin-tGap) & (VitalTypeAbbreviation=="PO2"), median(VitalResultNumeric),by="PatientICN"];colnames(PO2)[2]="PO2" Pain=Vo[(VitalSignTakenDateTime < Tpr-tGap & VitalSignTakenDateTime > Tpr-tWin-tGap) & (VitalTypeAbbreviation=="PN"), median(VitalResultNumeric),by="PatientICN"];colnames(Pain)[2]="Pain" lV=list(OC,Systolic,Diastolic,Weight,Pulse,Respiration,PO2,Pain) dV=Reduce(function(...) merge(...,all.x=T,by="PatientICN"),lV) for (i in seq(11,16)){dV[which(is.na(dV[,..i])),i]=0} #make sure columns 3-16 are the counts per condition columns. replace NA with 0 if (count==1){dV1=dV} if (count==2){dV2=dV} if (count==3){dV3=dV} if (count==4){dV4=dV} if (count==5){dV5=dV} if (count==6){dV6=dV} if (count==7){dV7=dV} if (count==8){dV8=dV} } # Labs LabL=readRDS(paste(c(C2020dir,"labs/Lipid.Sugars.",cohort,".RDS"),collapse="")) mLabL=merge(OC,LabL,all.x=T,by="PatientICN") LabL=mLabL LabK=readRDS(paste(c(C2020dir,"labs/M.kidney.protein.",cohort,".RDS"),collapse="")) mLabK=merge(OC,LabK,all.x=T,by="PatientICN") LabK=mLabK LabRBC=readRDS(paste(c(C2020dir,"labs/H.RBC.iron.",cohort,".RDS"),collapse="")) mLabRBC=merge(OC,LabRBC,all.x=T,by="PatientICN") LabRBC=mLabRBC for (count in c(1,2,3,4,5,6,7,8)){ tWin=window[count];tGap=gindow[count] A1C=LabL[(LabChemSpecimenDateTime < Tpr-tGap & LabChemSpecimenDateTime > Tpr-tWin-tGap) & (ShortName=="A1C"), median(LabChemResultNumericValue),by="PatientICN"];colnames(A1C)[2]="A1C" Glucose=LabL[(LabChemSpecimenDateTime < Tpr-tGap & LabChemSpecimenDateTime > Tpr-tWin-tGap) & (ShortName=="Glucose"), median(LabChemResultNumericValue),by="PatientICN"];colnames(Glucose)[2]="Glucose" HDLC=LabL[(LabChemSpecimenDateTime < Tpr-tGap & LabChemSpecimenDateTime > Tpr-tWin-tGap) & (ShortName=="HDLC"), median(LabChemResultNumericValue),by="PatientICN"];colnames(HDLC)[2]="HDLC" Hemoglobin=LabRBC[(LabChemSpecimenDateTime < Tpr-tGap & LabChemSpecimenDateTime > Tpr-tWin-tGap) & (ShortName=="Hemoglobin"), median(LabChemResultNumericValue),by="PatientICN"];colnames(Hemoglobin)[2]="Hemoglobin" LDLC=LabL[(LabChemSpecimenDateTime < Tpr-tGap & LabChemSpecimenDateTime > Tpr-tWin-tGap) & (ShortName=="LDLC"), median(LabChemResultNumericValue),by="PatientICN"];colnames(LDLC)[2]="LDLC" TotChol=LabL[(LabChemSpecimenDateTime < Tpr-tGap & LabChemSpecimenDateTime > Tpr-tWin-tGap) & (ShortName=="TotChol"), median(LabChemResultNumericValue),by="PatientICN"];colnames(TotChol)[2]="TotChol" Trig=LabL[(LabChemSpecimenDateTime < Tpr-tGap & LabChemSpecimenDateTime > Tpr-tWin-tGap) & (ShortName=="Trig"), median(LabChemResultNumericValue),by="PatientICN"];colnames(Trig)[2]="Trig" Creat=LabK[(LabChemSpecimenDateTime < Tpr-tGap & LabChemSpecimenDateTime > Tpr-tWin-tGap) & (ShortName=="Creat - BSP"), median(LabChemResultNumericValue),by="PatientICN"];colnames(Creat)[2]="Creat" eGFR=LabK[(LabChemSpecimenDateTime < Tpr-tGap & LabChemSpecimenDateTime > Tpr-tWin-tGap) & (ShortName=="eGFR"), median(LabChemResultNumericValue),by="PatientICN"];colnames(eGFR)[2]="eGFR" lLab=list(OC,A1C,Glucose,HDLC,Hemoglobin,LDLC,TotChol,Trig,Creat,eGFR) dLab=Reduce(function(...) merge(...,all.x=T,by="PatientICN"),lLab) for (i in seq(11,19)){dLab[which(is.na(dLab[,..i])),i]=0} #make sure columns 3-16 are the counts per condition columns. replace NA with 0 if (count==1){dLab1=dLab} if (count==2){dLab2=dLab} if (count==3){dLab3=dLab} if (count==4){dLab4=dLab} if (count==5){dLab5=dLab} if (count==6){dLab6=dLab} if (count==7){dLab7=dLab} if (count==8){dLab8=dLab} } save(dDemo1,dDemo2,dDemo3,dDemo4,dDemo5,dDemo6,dDemo7,dDemo8,dDx1,dDx2,dDx3,dDx4,dDx5,dDx6,dDx7,dDx8,dRx1,dRx2,dRx3,dRx4,dRx5,dRx6,dRx7,dRx8,dCPT1,dCPT2,dCPT3,dCPT4,dCPT5,dCPT6,dCPT7,dCPT8,dPHQ1,dPHQ2,dPHQ3,dPHQ4,dPHQ5,dPHQ6,dPHQ7,dPHQ8,dUtil1,dUtil2,dUtil3,dUtil4,dUtil5,dUtil6,dUtil7,dUtil8,dV1,dV2,dV3,dV4,dV5,dV6,dV7,dV8,dLab1,dLab2,dLab3,dLab4,dLab5,dLab6,dLab7,dLab8,file=paste(c("rv10.",cohort,".RData"),collapse="")) } # for loop over cohort } # end of pullData if (Run == "encode"){ ACS=readRDS("../ACS1.all.RDS") census=data.table(ACS$zip,ACS$Black/ACS$Pop.y,ACS$AmIndian/ACS$Pop.y,ACS$gt.75k.Income/ACS$Pop.ge25,ACS$Under0.5.x/ACS$Families,ACS$lt.HighSchoolGrad/ACS$Workers.ge16.x,ACS$Bachelors/ACS$Workers.ge16.x,ACS$Spanish/ACS$Workers.ge16,ACS$WithSSI/ACS$Pop.ge3,ACS$InSchool/ACS$Pop.Determined,ACS$mILDiff.18.34/ACS$mNoILDiff.18.34,ACS$mILDiff.35.64/ACS$mNoILDiff.35.64,ACS$mILDiff.65.74/ACS$mNoILDiff.65.74,ACS$"18.64.Vet.Pov"/ACS$"18.64.TotVet",ACS$"18.64.Vet.Pov.Disabil"/ACS$"18.64.Vet.Pov.NoDisabil") colnames(census)=c("PatientZIP","Frac_Black","Frac_AmInd","Frac_gt_75k","Frac_lt_HalfPovLine","Frac_NoHighSchool","Frac_Bachelors","Frac_Spanish","Frac_SSI","Frac_InSchool","Ratio_18_34_IndLivDiff-not","Ratio_35_64_IndLivDiff-not","Ratio_65_74_IndLivDiff-not","Frac_18-64-VetPov","Ratio_18-64-Pov-Disable_not") first=T for (cohort in seq(1,101,2)){ print(cohort) load(paste(c("rv10.",cohort,".RData"),collapse="")) dRx1$Combined_Analgesics_5=floor(dRx1$Combined_Analgesics_Rx/5) ;dRx1$CVD_5=floor(dRx1$CVD_Rx/5) dRx2$Combined_Analgesics_5=floor(dRx2$Combined_Analgesics_Rx/5) ;dRx2$CVD_5=floor(dRx2$CVD_Rx/5) dRx3$Combined_Analgesics_5=floor(dRx3$Combined_Analgesics_Rx/5) ;dRx3$CVD_5=floor(dRx3$CVD_Rx/5) dRx4$Combined_Analgesics_5=floor(dRx4$Combined_Analgesics_Rx/5) ;dRx4$CVD_5=floor(dRx4$CVD_Rx/5) dRx5$Combined_Analgesics_5=floor(dRx5$Combined_Analgesics_Rx/5) ;dRx5$CVD_5=floor(dRx5$CVD_Rx/5) dRx6$Combined_Analgesics_5=floor(dRx6$Combined_Analgesics_Rx/5) ;dRx6$CVD_5=floor(dRx6$CVD_Rx/5) dRx7$Combined_Analgesics_5=floor(dRx7$Combined_Analgesics_Rx/5) ;dRx7$CVD_5=floor(dRx7$CVD_Rx/5) dRx8$Combined_Analgesics_5=floor(dRx8$Combined_Analgesics_Rx/5) ;dRx8$CVD_5=floor(dRx8$CVD_Rx/5) Nt=8; Npt=dim(dDemo1)[1];NCat=1+5+9;NCatL=69 # 1 demo, 3 vital, 9 labs, 3 categories each NDx=76;NRx=18;NCPT=14;NUtil=14;NCod=NDx+NRx+NCPT+NUtil # Dx + Rx + CPT #dUtil1$SuicideAttempt=NULL;dUtil2$SuicideAttempt=NULL;dUtil3$SuicideAttempt=NULL Cat=array(rep(0, Npt*NCat*Nt),dim=c(Npt,NCat,Nt)) CatL=array(rep(0, Npt*NCatL*Nt),dim=c(Npt,NCatL,Nt)) ; dimnames(CatL)[[2]]=rep("Unassigned",NCatL) #presence / absence structured data Cod=array(rep(0, Npt*NCod*Nt),dim=c(Npt,NCod,Nt)); dimnames(Cod)[[2]]=c(colnames(dDx1)[11:(NDx+10)],colnames(dRx1)[11:(NRx+10)],colnames(dCPT1)[11:(NCPT+10)],colnames(dUtil1)[11:(NUtil+10)]) Cod[,,1]=as.matrix(cbind(dDx1[,11:(NDx+10)],dRx1[,11:(NRx+10)],dCPT1[,11:(NCPT+10)],dUtil1[,11:(NUtil+10)])) ; Cod[,,2]=as.matrix(cbind(dDx2[,11:(NDx+10)],dRx2[,11:(NRx+10)],dCPT2[,11:(NCPT+10)],dUtil2[,11:(NUtil+10)])) Cod[,,3]=as.matrix(cbind(dDx3[,11:(NDx+10)],dRx3[,11:(NRx+10)],dCPT3[,11:(NCPT+10)],dUtil3[,11:(NUtil+10)])) ; Cod[,,4]=as.matrix(cbind(dDx4[,11:(NDx+10)],dRx4[,11:(NRx+10)],dCPT4[,11:(NCPT+10)],dUtil4[,11:(NUtil+10)])) Cod[,,5]=as.matrix(cbind(dDx5[,11:(NDx+10)],dRx5[,11:(NRx+10)],dCPT5[,11:(NCPT+10)],dUtil5[,11:(NUtil+10)])) ; Cod[,,6]=as.matrix(cbind(dDx6[,11:(NDx+10)],dRx6[,11:(NRx+10)],dCPT6[,11:(NCPT+10)],dUtil6[,11:(NUtil+10)])) Cod[,,7]=as.matrix(cbind(dDx7[,11:(NDx+10)],dRx7[,11:(NRx+10)],dCPT7[,11:(NCPT+10)],dUtil7[,11:(NUtil+10)])) ; Cod[,,8]=as.matrix(cbind(dDx8[,11:(NDx+10)],dRx8[,11:(NRx+10)],dCPT8[,11:(NCPT+10)],dUtil8[,11:(NUtil+10)])) for (j in seq(1,NCod)){ for (k in seq(1,Nt)){Cod[Cod[,j,k]>1,j,k]=1} } # time-independent demographic variables dDemo=dDemo1[,c(2,4,5,6)] dDemo$Female=0;dDemo$Female[which(dDemo$Gender=="F")]=1 dDemo$Rac=0;dDemo$Rac[which(dDemo$Race=='BLACK OR AFRICAN AMERICAN')]=1 dDemo$Eth=0;dDemo$Eth[which(dDemo$Ethnicity=='HISPANIC OR LATINO')]=1 dDemo$Young=0;dDemo$Young[which((dDemo1$Tpr-dDemo1$BirthDate)< 55)]=1 dDemo$Old=0;dDemo$Old[which((dDemo1$Tpr-dDemo1$BirthDate)> 65)]=1 ################# Census data - zip code tmp=as.matrix(cbind(dDemo1[,13],dDemo2[,13],dDemo3[,13],dDemo4[,13],dDemo5[,13],dDemo6[,13],dDemo7[,13],dDemo8[,13])) #zip codes for (i in seq(1,8)){tmp[is.na(tmp[,i]),i]=0} for (i in seq(1,8)){Cat[tmp[,i]>0,13,i]=tmp[tmp[,i]>0,i]} for (i in seq(8,2)){tmp[tmp[,(i-1)]==0,(i-1)]=tmp[tmp[,(i-1)]==0,i]} for (i in seq(1,7)){tmp[tmp[,(i+1)]==0,(i+1)]=tmp[tmp[,(i+1)]==0,i]} for (i in seq(1,8)){tmp[tmp[,i]==0,i]=NA} z=data.table(tmp[,1]);colnames(z)="PatientZIP" ACSpt=merge(z,census,all.x=T) dDemo$Black_ACS=0;dDemo$Black_ACS[which(ACSpt$Frac_Black > 0.2)]=1 dDemo$AmInd_ACS=0;dDemo$AmInd_ACS[which(ACSpt$Frac_AmInd > 0.02)]=1 dDemo$gt75k_ACS=0;dDemo$gt75k_ACS[which(ACSpt$Frac_gt_75k > 0.2)]=1 dDemo$ltHalfPov_ACS =0;dDemo$ltHalfPov_ACS[which(ACSpt$Frac_lt_HalfPovLine > 0.4)]=1 dDemo$NoHighSchool_ACS =0;dDemo$NoHighSchool_ACS[which(ACSpt$Frac_NoHighSchool > 0.4)]=1 dDemo$Bachelors_ACS =0;dDemo$Bachelors_ACS[which(ACSpt$Frac_Bachelors > 0.4)]=1 dDemo$Spanish_ACS =0;dDemo$Spanish_ACS[which(ACSpt$Frac_Spanish > 0.07)]=1 dDemo$SSI_ACS =0;dDemo$SSI_ACS[which(ACSpt$Frac_SSI > 0.1)]=1 dDemo$NotInSchool_ACS =0;dDemo$NotInSchool_ACS[which(ACSpt$Frac_InSchool < 0.18)]=1 dDemo$DifficIndivLiv18_34_ACS =0;dDemo$DifficIndivLiv18_34_ACS[which(ACSpt$"Ratio_18_34_IndLivDiff-not" > 0.05)]=1 dDemo$DifficIndivLiv35_64_ACS =0;dDemo$DifficIndivLiv35_64_ACS[which(ACSpt$"Ratio_35_64_IndLivDiff-not" > 0.1)]=1 dDemo$DifficIndivLiv65_74_ACS =0;dDemo$DifficIndivLiv65_74_ACS[which(ACSpt$"Ratio_65_74_IndLivDiff-not" >0.15)]=1 dDemo$Vet_Pov_18_64_ACS =0;dDemo$Vet_Pov_18_64_ACS[which(ACSpt$"Frac_18-64-VetPov" > 0.2)]=1 dDemo$Vet_Pov_Dis_18_64_ACS =0;dDemo$Vet_Pov_Dis_18_64_ACS[which(ACSpt$"Ratio_18-64-Pov-Disable_not" > 0.4)]=1 #column="Troponin" #tmp=as.matrix(cbind(dLab1[,..column],dLab2[,..column],dLab3[,..column],dLab4[,..column],dLab5[,..column],dLab6[,..column],dLab7[,..column],dLab8[,..column])) #Troponin # time-dependent categorical and numerical variables. column="MaritalStatus" Married=as.matrix(cbind(dDemo1[,..column],dDemo2[,..column],dDemo3[,..column],dDemo4[,..column],dDemo5[,..column],dDemo6[,..column],dDemo7[,..column],dDemo8[,..column])) for (i in seq(8,2)){Married[is.na(match(Married[,(i-1)],c("D","M","N","S","U","W"))),(i-1)]= Married[is.na(match(Married[,(i-1)],c("D","M","N","S","U","W"))),(i)] } #carry entry forward in time for (i in seq(1,7)){Married[is.na(match(Married[,(i+1)],c("D","M","N","S","U","W"))),(i+1)]= Married[is.na(match(Married[,(i+1)],c("D","M","N","S","U","W"))),(i)] } #carry entry backwards in time for (i in seq(1,8)){CatL[Married[,i]=="M",1,i]=1} ; dimnames(CatL)[[2]][[1]]="Married" # marital status ... 'M' = 1 otherwise = 0 for (i in seq(1,7)){CatL[(Married[,i]=="M") & (Married[,(i+1)]!="M"),2,i]=1}# change in marital status, getting married for (i in seq(1,7)){CatL[(Married[,i]!="M") & (Married[,(i+1)]=="M"),2,i]=1} ; dimnames(CatL)[[2]][[2]]="MaritalChange" # change in marital status, getting unmarried column1="Weight";column2="Height" BMI=as.matrix(cbind(703*dV1[,..column1]/dV1[,..column2]^2,703*dV2[,..column1]/dV2[,..column2]^2,703*dV3[,..column1]/dV3[,..column2]^2,703*dV4[,..column1]/dV4[,..column2]^2,703*dV5[,..column1]/dV5[,..column2]^2,703*dV6[,..column1]/dV6[,..column2]^2,703*dV7[,..column1]/dV7[,..column2]^2,703*dV8[,..column1]/dV8[,..column2]^2)) for (i in seq(1,8)){BMI[BMI[,i] < 10 | BMI[,i] > 100,i] = 0} for (i in seq(1,8)){BMI[is.na(BMI[,i]),i]=0} for (i in seq(8,2)){BMI[BMI[,(i-1)]==0,(i-1)]=BMI[BMI[,(i-1)]==0,i]} for (i in seq(1,7)){BMI[BMI[,(i+1)]==0,(i+1)]=BMI[BMI[,(i+1)]==0,i]} for (i in seq(1,8)){BMI[BMI[,i]==0,i]=NA} Cat[,2,]=BMI for (i in seq(1,8)){CatL[BMI[,i] < 18.5,3,i]=1} ; dimnames(CatL)[[2]][[3]]="Underweight" for (i in seq(1,8)){CatL[BMI[,i] > 25 & BMI[,i] < 30,4,i]=1} ; dimnames(CatL)[[2]][[4]]="Overweight" for (i in seq(1,8)){CatL[BMI[,i] > 30,5,i]=1} ; dimnames(CatL)[[2]][[5]]="Obese" column1="Systolic" tmp=as.matrix(cbind(dV1[,..column1],dV2[,..column1],dV3[,..column1],dV4[,..column1],dV5[,..column1],dV6[,..column1],dV7[,..column1],dV8[,..column1])) #Systolic for (i in seq(1,8)){tmp[is.na(tmp[,i]),i]=0} for (i in seq(8,2)){tmp[tmp[,(i-1)]==0,(i-1)]=tmp[tmp[,(i-1)]==0,i]} for (i in seq(1,7)){tmp[tmp[,(i+1)]==0,(i+1)]=tmp[tmp[,(i+1)]==0,i]} for (i in seq(1,8)){tmp[tmp[,i]==0,i]=NA} Cat[,3,]=tmp # Systolic for (i in seq(1,8)){CatL[tmp[,i]<90,6,i]=1} ; dimnames(CatL)[[2]][[6]]="Hypotension" for (i in seq(1,8)){CatL[(tmp[,i]>120 & tmp[,i]<130),7,i]=1} ; dimnames(CatL)[[2]][[7]]="Elevated BP" for (i in seq(1,8)){CatL[(tmp[,i]>130 & tmp[,i]<140),8,i]=1} ; dimnames(CatL)[[2]][[8]]="Hypertension_1" for (i in seq(1,8)){CatL[tmp[,i]>140,9,i]=1} ; dimnames(CatL)[[2]][[9]]="Hypertension_2" column2="Diastolic" tmp=as.matrix(cbind(dV1[,..column1]-dV1[,..column2],dV2[,..column1]-dV2[,..column2],dV3[,..column1]-dV3[,..column2],dV4[,..column1]-dV4[,..column2],dV5[,..column1]-dV5[,..column2],dV6[,..column1]-dV6[,..column2],dV7[,..column1]-dV7[,..column2],dV8[,..column1]-dV8[,..column2])) # pulse pressure for (i in seq(1,8)){tmp[is.na(tmp[,i]),i]=0} for (i in seq(8,2)){tmp[tmp[,(i-1)]==0,(i-1)]=tmp[tmp[,(i-1)]==0,i]} for (i in seq(1,7)){tmp[tmp[,(i+1)]==0,(i+1)]=tmp[tmp[,(i+1)]==0,i]} for (i in seq(1,8)){tmp[tmp[,i]==0,i]=NA} Cat[,4,]=tmp # PulsePressure for (i in seq(1,8)){CatL[tmp[,i]<30,10,i]=1} ; dimnames(CatL)[[2]][[10]]="PulsePr.lt.30" for (i in seq(1,8)){CatL[tmp[,i]>50 & tmp[,i]<70, 11,i]=1} ; dimnames(CatL)[[2]][[11]]="50.lt.PulsePr.lt.70" for (i in seq(1,8)){CatL[tmp[,i]>70 & tmp[,i]<90, 12,i]=1} ; dimnames(CatL)[[2]][[12]]="70.lt.PulsePr.lt.90" for (i in seq(1,8)){CatL[tmp[,i]>90,13,i]=1} ; dimnames(CatL)[[2]][[13]]="90.lt.PulsePr" column="Pulse" tmp=as.matrix(cbind(dV1[,..column],dV2[,..column],dV3[,..column],dV4[,..column],dV5[,..column],dV6[,..column],dV7[,..column],dV8[,..column])) #Pulse for (i in seq(1,8)){tmp[is.na(tmp[,i]),i]=0} for (i in seq(8,2)){tmp[tmp[,(i-1)]==0,(i-1)]=tmp[tmp[,(i-1)]==0,i]} for (i in seq(1,7)){tmp[tmp[,(i+1)]==0,(i+1)]=tmp[tmp[,(i+1)]==0,i]} for (i in seq(1,8)){tmp[tmp[,i]==0,i]=NA} Cat[,5,]=tmp # PulseRate Gender=dDemo1$Gender for (i in seq(1,8)){CatL[Gender=="M" & tmp[,i]<66,14,i]=1} ; dimnames(CatL)[[2]][[14]]="LowPulseRate" for (i in seq(1,8)){CatL[Gender=="F" & tmp[,i]<70,14,i]=1} for (i in seq(1,8)){CatL[Gender=="M" & tmp[,i]>82,15,i]=1} ; dimnames(CatL)[[2]][[15]]="HighPulseRate" for (i in seq(1,8)){CatL[Gender=="F" & tmp[,i]>85,15,i]=1} column="PO2" tmp=as.matrix(cbind(dV1[,..column],dV2[,..column],dV3[,..column],dV4[,..column],dV5[,..column],dV6[,..column],dV7[,..column],dV8[,..column])) #PO2 #tmp=as.matrix(cbind(dV1[,14],dV2[,14],dV3[,14],dV4[,14],dV5[,14],dV6[,14],dV7[,14],dV8[,14])) #PO2 for (i in seq(1,8)){tmp[is.na(tmp[,i]),i]=0} for (i in seq(8,2)){tmp[tmp[,(i-1)]==0,(i-1)]=tmp[tmp[,(i-1)]==0,i]} for (i in seq(1,7)){tmp[tmp[,(i+1)]==0,(i+1)]=tmp[tmp[,(i+1)]==0,i]} for (i in seq(1,8)){tmp[tmp[,i]==0,i]=NA} Cat[,6,]=tmp # PO2 for (i in seq(1,8)){CatL[is.na(tmp[,i]),16,i]=1} ; dimnames(CatL)[[2]][[16]]="PO2.NA" for (i in seq(1,8)){CatL[tmp[,i]<88,17,i]=1} ; dimnames(CatL)[[2]][[17]]="PO2.lt.88" for (i in seq(1,8)){CatL[tmp[,i]>=88 & tmp[,i]<91,18,i]=1} ; dimnames(CatL)[[2]][[18]]="88.le.PO2.lt.91" for (i in seq(1,8)){CatL[tmp[,i]>=91 & tmp[,i]<95,19,i]=1} ; dimnames(CatL)[[2]][[19]]="91.le.PO2.lt.95" #reference is ge 95 column="Pain" tmp=as.matrix(cbind(dV1[,..column],dV2[,..column],dV3[,..column],dV4[,..column],dV5[,..column],dV6[,..column],dV7[,..column],dV8[,..column])) #Pain for (i in seq(1,8)){tmp[tmp[,i]==0,i]=NA} Cat[,7,]=tmp # Pain for (i in seq(1,8)){CatL[is.na(tmp[,i]),20,i]=1} ; dimnames(CatL)[[2]][[20]]="Pain.NA" #for (i in seq(1,8)){CatL[tmp[,i]==0,14,i]=1} ; dimnames(CatL)[[2]][[14]]="No Pain" #reference for (i in seq(1,8)){CatL[tmp[,i]>=1 & tmp[,i]<5,21,i]=1} ; dimnames(CatL)[[2]][[21]]="1.le.Pain.lt.5" for (i in seq(1,8)){CatL[tmp[,i]>=5,22,i]=1} ; dimnames(CatL)[[2]][[22]]="Pain.ge.5" Ndv=22 ####### Begin Labs ######## column="A1C" tmp=as.matrix(cbind(dLab1[,..column],dLab2[,..column],dLab3[,..column],dLab4[,..column],dLab5[,..column],dLab6[,..column],dLab7[,..column],dLab8[,..column])) #A1C for (i in seq(1,8)){tmp[is.na(tmp[,i]),i]=0} for (i in seq(8,2)){tmp[tmp[,(i-1)]==0,(i-1)]=tmp[tmp[,(i-1)]==0,i]} for (i in seq(1,7)){tmp[tmp[,(i+1)]==0,(i+1)]=tmp[tmp[,(i+1)]==0,i]} for (i in seq(1,8)){tmp[tmp[,i]==0,i]=NA} Cat[,8,]=tmp # A1C #for (i in seq(1,8)){CatL[tmp[,i]<6,16,i]=1} ; dimnames(CatL)[[2]][[16]]="HbA1c.lt.6" #reference for (i in seq(1,8)){CatL[is.na(tmp[,i]),(Ndv+1),i]=1} ; dimnames(CatL)[[2]][[(Ndv+1)]]="HbA1c.NA" for (i in seq(1,8)){CatL[(tmp[,i]>6 & tmp[,i]<7),(Ndv+2),i]=1} ; dimnames(CatL)[[2]][[(Ndv+2)]]="6.lt.HbA1c.lt.7" for (i in seq(1,8)){CatL[(tmp[,i]>=7 & tmp[,i]<8),(Ndv+3),i]=1} ; dimnames(CatL)[[2]][[(Ndv+3)]]="7.le.HbA1c.lt.8" for (i in seq(1,8)){CatL[(tmp[,i]>=8 & tmp[,i]<9),(Ndv+4),i]=1} ; dimnames(CatL)[[2]][[(Ndv+4)]]="8.le.HbA1c.lt.9" for (i in seq(1,8)){CatL[tmp[,i]>9,(Ndv+5),i]=1} ; dimnames(CatL)[[2]][[(Ndv+5)]]="9.lt.HbA1c" ############ column="Glucose" tmp=as.matrix(cbind(dLab1[,..column],dLab2[,..column],dLab3[,..column],dLab4[,..column],dLab5[,..column],dLab6[,..column],dLab7[,..column],dLab8[,..column])) #Glucose for (i in seq(1,8)){tmp[is.na(tmp[,i]),i]=0} for (i in seq(8,2)){tmp[tmp[,(i-1)]==0,(i-1)]=tmp[tmp[,(i-1)]==0,i]} for (i in seq(1,7)){tmp[tmp[,(i+1)]==0,(i+1)]=tmp[tmp[,(i+1)]==0,i]} for (i in seq(1,8)){tmp[tmp[,i]==0,i]=NA} Cat[,9,]=tmp # Glucose for (i in seq(1,8)){CatL[is.na(tmp[,i]),(Ndv+6),i]=1} ; dimnames(CatL)[[2]][[(Ndv+6)]]="Glucose.NA" for (i in seq(1,8)){CatL[tmp[,i]<70,(Ndv+7),i]=1} ; dimnames(CatL)[[2]][[Ndv+7]]="Glucose.lt.70" #for (i in seq(1,8)){CatL[tmp[,i]>6,23,i]=1} ; dimnames(CatL)[[2]][[23]]="6.lt.HbA1c.lt.7" reference for (i in seq(1,8)){CatL[tmp[,i]>=100 & tmp[,i]<125,(Ndv+8),i]=1} ; dimnames(CatL)[[2]][[(Ndv+8)]]="100.le.Glucose.lt.125" for (i in seq(1,8)){CatL[tmp[,i]>=125 & tmp[,i]<150,(Ndv+9),i]=1} ; dimnames(CatL)[[2]][[(Ndv+9)]]="125.le.Glucose.lt.150" for (i in seq(1,8)){CatL[tmp[,i]>=150,(Ndv+10),i]=1} ; dimnames(CatL)[[2]][[(Ndv+10)]]="150.le.Glucose" ############### column="Hemoglobin" tmp=as.matrix(cbind(dLab1[,..column],dLab2[,..column],dLab3[,..column],dLab4[,..column],dLab5[,..column],dLab6[,..column],dLab7[,..column],dLab8[,..column])) #Hemoglobin for (i in seq(1,8)){tmp[is.na(tmp[,i]),i]=0} for (i in seq(8,2)){tmp[tmp[,(i-1)]==0,(i-1)]=tmp[tmp[,(i-1)]==0,i]} for (i in seq(1,7)){tmp[tmp[,(i+1)]==0,(i+1)]=tmp[tmp[,(i+1)]==0,i]} for (i in seq(1,8)){tmp[tmp[,i]==0,i]=NA} Cat[,10,]=tmp # Hemoglobin for (i in seq(1,8)){CatL[is.na(tmp[,i]),(Ndv+11),i]=1} ; dimnames(CatL)[[2]][[(Ndv+11)]]="Hemoglobin.NA" for (i in seq(1,8)){CatL[tmp[,i]<8,(Ndv+12),i]=1} ; dimnames(CatL)[[2]][[(Ndv+12)]]="Hemoglobin.lt.8" for (i in seq(1,8)){CatL[tmp[,i]>=8 & tmp[,i]<10,(Ndv+13),i]=1} ; dimnames(CatL)[[2]][[(Ndv+13)]]="8.le.Hemoglobin.lt.10" for (i in seq(1,8)){CatL[tmp[,i]>=10 & tmp[,i]<11,(Ndv+14),i]=1} ; dimnames(CatL)[[2]][[(Ndv+14)]]="10.le.Hemoglobin.lt.11" #for (i in seq(1,8)){CatL[tmp[,i]>11,24,i]=1} ; dimnames(CatL)[[2]][[24]]="Hemoglobin.gt.11" reference ############## column="HDLC" tmp=as.matrix(cbind(dLab1[,..column],dLab2[,..column],dLab3[,..column],dLab4[,..column],dLab5[,..column],dLab6[,..column],dLab7[,..column],dLab8[,..column])) #HDLC for (i in seq(1,8)){tmp[is.na(tmp[,i]),i]=0} for (i in seq(8,2)){tmp[tmp[,(i-1)]==0,(i-1)]=tmp[tmp[,(i-1)]==0,i]} for (i in seq(1,7)){tmp[tmp[,(i+1)]==0,(i+1)]=tmp[tmp[,(i+1)]==0,i]} for (i in seq(1,8)){tmp[tmp[,i]==0,i]=NA} Cat[,11,]=tmp # HDL for (i in seq(1,8)){CatL[is.na(tmp[,i]),(Ndv+15),i]=1} ; dimnames(CatL)[[2]][[(Ndv+15)]]="HDL.NA" #for (i in seq(1,8)){CatL[is.na(tmp[,i]),26,i]=1} ; dimnames(CatL)[[2]][[26]]="HDL.lt.100" #reference for (i in seq(1,8)){CatL[tmp[,i]<40,(Ndv+16),i]=1} ; dimnames(CatL)[[2]][[(Ndv+16)]]="HDL.lt.40" for (i in seq(1,8)){CatL[tmp[,i]>=60 & tmp[,i]<90,(Ndv+17),i]=1} ; dimnames(CatL)[[2]][[(Ndv+17)]]="40.le.HDL.lt.60" for (i in seq(1,8)){CatL[tmp[,i]>=90,(Ndv+18),i]=1} ; dimnames(CatL)[[2]][[(Ndv+18)]]="HDL.ge.90" #################### column="LDLC" tmp=as.matrix(cbind(dLab1[,..column],dLab2[,..column],dLab3[,..column],dLab4[,..column],dLab5[,..column],dLab6[,..column],dLab7[,..column],dLab8[,..column])) #LDLC for (i in seq(1,8)){tmp[is.na(tmp[,i]),i]=0} for (i in seq(8,2)){tmp[tmp[,(i-1)]==0,(i-1)]=tmp[tmp[,(i-1)]==0,i]} for (i in seq(1,7)){tmp[tmp[,(i+1)]==0,(i+1)]=tmp[tmp[,(i+1)]==0,i]} for (i in seq(1,8)){tmp[tmp[,i]==0,i]=NA} Cat[,12,]=tmp # LDL for (i in seq(1,8)){CatL[is.na(tmp[,i]),(Ndv+19),i]=1} ; dimnames(CatL)[[2]][[(Ndv+19)]]="LDL.NA" #for (i in seq(1,8)){CatL[is.na(tmp[,i]),26,i]=1} ; dimnames(CatL)[[2]][[26]]="LDL.lt.100" #reference for (i in seq(1,8)){CatL[tmp[,i]>100 & tmp[,i]<130,(Ndv+20),i]=1} ; dimnames(CatL)[[2]][[(Ndv+20)]]="LDL.lt.130" for (i in seq(1,8)){CatL[tmp[,i]>=130 & tmp[,i]<160,(Ndv+21),i]=1} ; dimnames(CatL)[[2]][[(Ndv+21)]]="130.le.LDL.lt.160" for (i in seq(1,8)){CatL[tmp[,i]>=160 & tmp[,i]<190,(Ndv+22),i]=1} ; dimnames(CatL)[[2]][[(Ndv+22)]]="160.le.LDL.lt.190" for (i in seq(1,8)){CatL[tmp[,i]>=190,(Ndv+23),i]=1} ; dimnames(CatL)[[2]][[(Ndv+23)]]="190.le.LDL" ################### column="TotChol" tmp=as.matrix(cbind(dLab1[,..column],dLab2[,..column],dLab3[,..column],dLab4[,..column],dLab5[,..column],dLab6[,..column],dLab7[,..column],dLab8[,..column])) #TotChol for (i in seq(1,8)){tmp[is.na(tmp[,i]),i]=0} for (i in seq(8,2)){tmp[tmp[,(i-1)]==0,(i-1)]=tmp[tmp[,(i-1)]==0,i]} for (i in seq(1,7)){tmp[tmp[,(i+1)]==0,(i+1)]=tmp[tmp[,(i+1)]==0,i]} for (i in seq(1,8)){tmp[tmp[,i]==0,i]=NA} Cat[,13,]=tmp # Tot Chol for (i in seq(1,8)){CatL[is.na(tmp[,i]),(Ndv+24),i]=1} ; dimnames(CatL)[[2]][[(Ndv+24)]]="TotChol.NA" #for (i in seq(1,8)){CatL[is.na(tmp[,i]),26,i]=1} ; dimnames(CatL)[[2]][[26]]="TotChol.lt.200" #reference for (i in seq(1,8)){CatL[tmp[,i]>200 & tmp[,i]<220,(Ndv+25),i]=1} ; dimnames(CatL)[[2]][[(Ndv+25)]]="200.lt.TotChol.lt.220" for (i in seq(1,8)){CatL[tmp[,i]>=220,(Ndv+26),i]=1} ; dimnames(CatL)[[2]][[(Ndv+26)]]="220.le.TotChol" ##################### column="Trig" tmp=as.matrix(cbind(dLab1[,..column],dLab2[,..column],dLab3[,..column],dLab4[,..column],dLab5[,..column],dLab6[,..column],dLab7[,..column],dLab8[,..column])) #Trig for (i in seq(1,8)){tmp[is.na(tmp[,i]),i]=0} for (i in seq(8,2)){tmp[tmp[,(i-1)]==0,(i-1)]=tmp[tmp[,(i-1)]==0,i]} for (i in seq(1,7)){tmp[tmp[,(i+1)]==0,(i+1)]=tmp[tmp[,(i+1)]==0,i]} for (i in seq(1,8)){tmp[tmp[,i]==0,i]=NA} Cat[,14,]=tmp # Trig for (i in seq(1,8)){CatL[is.na(tmp[,i]),(Ndv+27),i]=1} ; dimnames(CatL)[[2]][[(Ndv+27)]]="Trig.NA" #for (i in seq(1,8)){CatL[is.na(tmp[,i]),26,i]=1} ; dimnames(CatL)[[2]][[26]]="Trig.lt.150" #reference for (i in seq(1,8)){CatL[tmp[,i]>150 & tmp[,i]<200,(Ndv+28),i]=1} ; dimnames(CatL)[[2]][[Ndv+28]]="150.lt.Trig.lt.200" for (i in seq(1,8)){CatL[tmp[,i]>=200,(Ndv+29),i]=1} ; dimnames(CatL)[[2]][[(Ndv+29)]]="200.le.Trig" #################### column="Creat" tmp=as.matrix(cbind(dLab1[,..column],dLab2[,..column],dLab3[,..column],dLab4[,..column],dLab5[,..column],dLab6[,..column],dLab7[,..column],dLab8[,..column])) #Creat #tmp=as.matrix(cbind(dLab1[,17],dLab2[,17],dLab3[,17],dLab4[,17],dLab5[,17],dLab6[,17],dLab7[,17],dLab8[,17])) #Creat for (i in seq(1,8)){tmp[is.na(tmp[,i]),i]=0} for (i in seq(1,8)){tmp[tmp[,i] > 10,i]=tmp[tmp[,i] > 10,i]/88.42} #convert mmol/l to g/dL, based on ranges. need to do properly when pulling data for (i in seq(1,8)){tmp[tmp[,i] > 3 & tmp[,i] < 10,i]=0} for (i in seq(1,8)){tmp[,i]=186*tmp[,i]^(-1.154)*(dDemo1$Tpr-dDemo1$BirthDate)^(-0.203)} #formula for eGFR for non-black males for (i in seq(1,8)){tmp[dDemo$Female==1,i]=0.742*tmp[dDemo$Female==1,i]} for (i in seq(1,8)){tmp[dDemo$Rac==1,i]=1.210*tmp[dDemo$Rac==1,i]} for (i in seq(1,8)){tmp[tmp[,i]==Inf,i]=0} Cat[,15,]=tmp # Creat converted to eGFR --- merge with reported eGFR, below column="eGFR" tmp=as.matrix(cbind(dLab1[,..column],dLab2[,..column],dLab3[,..column],dLab4[,..column],dLab5[,..column],dLab6[,..column],dLab7[,..column],dLab8[,..column])) #eGFR for (i in seq(1,8)){tmp[is.na(tmp[,i]),i]=0} for (i in seq(1,8)){Cat[tmp[,i]>0,13,i]=tmp[tmp[,i]>0,i]} for (i in seq(8,2)){tmp[tmp[,(i-1)]==0,(i-1)]=tmp[tmp[,(i-1)]==0,i]} for (i in seq(1,7)){tmp[tmp[,(i+1)]==0,(i+1)]=tmp[tmp[,(i+1)]==0,i]} for (i in seq(1,8)){tmp[tmp[,i]==0,i]=NA} Cat[,15,]=tmp # eGFR combined with creatinine converted to eGFR for (i in seq(1,8)){CatL[is.na(tmp[,i]),45,i]=1} ; dimnames(CatL)[[2]][[Ndv+30]]="eGFR.NA" #for (i in seq(1,8)){CatL[tmp[,i]>90,38,i]=1} ; dimnames(CatL)[[2]][[38]]= #reference for (i in seq(1,8)){CatL[tmp[,i]>=60 & tmp[,i]<90,(Ndv+30),i]=1} ; dimnames(CatL)[[2]][[(Ndv+31)]]="60.lt.eGFR.lt.90" for (i in seq(1,8)){CatL[tmp[,i]>=45 & tmp[,i]<60,(Ndv+31),i]=1} ; dimnames(CatL)[[2]][[(Ndv+31)]]="45.lt.eGFR.lt.60" for (i in seq(1,8)){CatL[tmp[,i]>=30 & tmp[,i]<45,(Ndv+32),i]=1} ; dimnames(CatL)[[2]][[(Ndv+32)]]="30.lt.eGFR.lt.45" for (i in seq(1,8)){CatL[tmp[,i]>=15 & tmp[,i]<30,(Ndv+33),i]=1} ; dimnames(CatL)[[2]][[(Ndv+33)]]="15.lt.eGFR.lt.30" for (i in seq(1,8)){CatL[tmp[,i]<15,(Ndv+34),i]=1} ; dimnames(CatL)[[2]][[(Ndv+34)]]="15.gt.eGFR" ################### Ndvl=Ndv+34 ################### column="PHQ2_Count" tmp=as.matrix(cbind(dPHQ1[,..column],dPHQ2[,..column],dPHQ3[,..column],dPHQ4[,..column],dPHQ5[,..column],dPHQ6[,..column],dPHQ7[,..column],dPHQ8[,..column])) #PHQ2 administered for (i in seq(1,8)){CatL[tmp[,i]==0 ,(Ndvl+1),i]=1} ; dimnames(CatL)[[2]][[(Ndvl+1)]]="PHQ2.not.administered" column="PHQ2_Raw.max" tmp=as.matrix(cbind(dPHQ1[,..column],dPHQ2[,..column],dPHQ3[,..column],dPHQ4[,..column],dPHQ5[,..column],dPHQ6[,..column],dPHQ7[,..column],dPHQ8[,..column])) #PHQ2 max score for (i in seq(1,8)){CatL[tmp[,i]>=1 & tmp[,i]<3,(Ndvl+2),i]=1} ; dimnames(CatL)[[2]][[(Ndvl+2)]]="PHQ2.eq.1.2" #reference is PHQ2 = 0 for (i in seq(1,8)){CatL[tmp[,i]>=3 & tmp[,i]<5,(Ndvl+3),i]=1} ; dimnames(CatL)[[2]][[(Ndvl+3)]]="PHQ2.eq.3.4" for (i in seq(1,8)){CatL[tmp[,i]>=5 ,(Ndvl+4),i]=1} ; dimnames(CatL)[[2]][[(Ndvl+4)]]="PHQ2.eq.5.6" ################### column="PHQ9_Count" tmp=as.matrix(cbind(dPHQ1[,..column],dPHQ2[,..column],dPHQ3[,..column],dPHQ4[,..column],dPHQ5[,..column],dPHQ6[,..column],dPHQ7[,..column],dPHQ8[,..column])) #PHQ9 administered for (i in seq(1,8)){CatL[tmp[,i]==0 ,(Ndvl+5),i]=1} ; dimnames(CatL)[[2]][[(Ndvl+5)]]="PHQ9.not.administered" column="PHQ9_Raw.max" tmp=as.matrix(cbind(dPHQ1[,..column],dPHQ2[,..column],dPHQ3[,..column],dPHQ4[,..column],dPHQ5[,..column],dPHQ6[,..column],dPHQ7[,..column],dPHQ8[,..column])) #PHQ9 max score for (i in seq(1,8)){CatL[tmp[,i]>=5 & tmp[,i]<10,(Ndvl+6),i]=1} ; dimnames(CatL)[[2]][[(Ndvl+6)]]="PHQ9.score.5-9" for (i in seq(1,8)){CatL[tmp[,i]>=10 & tmp[,i]<15,(Ndvl+7),i]=1} ; dimnames(CatL)[[2]][[(Ndvl+7)]]="PHQ9.score.10-14" for (i in seq(1,8)){CatL[tmp[,i]>=15 & tmp[,i]<20,(Ndvl+8),i]=1} ; dimnames(CatL)[[2]][[(Ndvl+8)]]="PHQ9.score.15.19" for (i in seq(1,8)){CatL[tmp[,i]>=20,(Ndvl+9),i]=1} ; dimnames(CatL)[[2]][[(Ndvl+9)]]="PHQ9.20-27" ################# column="PHQ9_Q9.recent1" tmp=as.matrix(cbind(dPHQ1[,..column],dPHQ2[,..column],dPHQ3[,..column],dPHQ4[,..column],dPHQ5[,..column],dPHQ6[,..column],dPHQ7[,..column],dPHQ8[,..column])) #PHQ9 Q9=1 for (i in seq(1,8)){CatL[tmp[,i]==1,(Ndvl+10),i]=1} ; dimnames(CatL)[[2]][[(Ndvl+10)]]="PHQ9.Q9.score=1" column="PHQ9_Q9.recent2" tmp=as.matrix(cbind(dPHQ1[,..column],dPHQ2[,..column],dPHQ3[,..column],dPHQ4[,..column],dPHQ5[,..column],dPHQ6[,..column],dPHQ7[,..column],dPHQ8[,..column])) #PHQ9 Q9=2 for (i in seq(1,8)){CatL[tmp[,i]==1,(Ndvl+11),i]=1} ; dimnames(CatL)[[2]][[(Ndvl+11)]]="PHQ9.Q9.score=2" column="PHQ9_Q9.recent3" tmp=as.matrix(cbind(dPHQ1[,..column],dPHQ2[,..column],dPHQ3[,..column],dPHQ4[,..column],dPHQ5[,..column],dPHQ6[,..column],dPHQ7[,..column],dPHQ8[,..column])) #PHQ9 Q9=3 for (i in seq(1,8)){CatL[tmp[,i]==1,(Ndvl+12),i]=1} ; dimnames(CatL)[[2]][[(Ndvl+12)]]="PHQ9.Q9.score=3" column="PHQ9_Q9.recent4" tmp=as.matrix(cbind(dPHQ1[,..column],dPHQ2[,..column],dPHQ3[,..column],dPHQ4[,..column],dPHQ5[,..column],dPHQ6[,..column],dPHQ7[,..column],dPHQ8[,..column])) #PHQ9 Q9=4 for (i in seq(1,8)){CatL[tmp[,i]==1,(Ndvl+13),i]=1} ; dimnames(CatL)[[2]][[(Ndvl+13)]]="PHQ9.Q9.score=4" names=c(dimnames(CatL)[[2]],dimnames(Cod)[[2]]) Nn=length(names) m1=cbind(as.matrix(dDemo[,5:23]),CatL[,,1],Cod[,,1],CatL[,,2],Cod[,,2],CatL[,,3],Cod[,,3],CatL[,,4],Cod[,,4],CatL[,,5],Cod[,,5],CatL[,,6],Cod[,,6],CatL[,,7],Cod[,,7],CatL[,,8],Cod[,,8]) rownames(m1)=dRx1$PatientICN ext=c(".a",".b",".c",".d",".e",".f",".g",".h") for (i in seq(1,8)){colnames(m1)[(20+(i-1)*Nn):(19+i*Nn)]=paste(names,ext[i],sep="")} y1=dDemo1[,c("Tpr","OCDate","OC")] if (first == T){m=m1;y=y1; first=F} if (first == F){m=rbind(m,m1);y=rbind(y,y1)} } } yss=y;save(yss,m,file="rv10.2.100.RData") ##################################### filter out controls, and filter according to time gap bewteen Tpr and OCDate, then run glm / cox model with model selection if (Run == "glm"){ library(survival) library(data.table) library(glmnet) options("width"=420) w=which(y$OCDate<23 & y$OCDate > y$Tpr) #for rv8m this should be fixed in input file generation... km=m[w,];ky=yss[w] #for survival, skip next 7 lines y1=unlist(ky$OC) # basic fix, for glm y1[y1==2]=1; y1[y1==3]=1 y=as.factor(y1) y=as.vector(y) #fit.cv=cv.glmnet(km,ky,family="binomial") #fit.cv=cv.glmnet(m,y,family="binomial") w=which(yss$OCDate > 19 & yss$OC > 0) #move 1/4 of cases after 2 years to controls, and ignore rest. nw=sample(w,.75*length(w)) yss[w,"OC"]=0 ky=yss[-nw,] km=m[-nw,] y=ky$OC ; y[y>1]=0 #fit.cv=cv.glmnet(km,y,family="binomial") #for cox model, strip times earlier than 0.12 years after Tpr, then downselect controls to managable calculation size (< 1M patients) #load("rv8.1.99.RData") y=yss y$time=y$OCDate-y$Tpr # time, for survival analysis w=which(ky$time > 0.12) km=m[w,] ky=y[w,] y1=y$OC;y1[y1 > 1]=1 #wY=which(y1==1) ; wZ=sample(which(y1==0),4*length(wY)) ; w=c(wY,wZ) #km=m[w,];ky=y[w,] #y1=ky$OC;y1[y1 > 1]=1 ys=ky;ys$OC[ys$OC>1]=1 #m=km;y=ky yss=Surv(ys[,time],ys[,OC]) save(ky,km,yss,file="data.rv10m.even.RData") rm(m) cv.fit=cv.glmnet(km,yss,family="cox",nfolds=5) #fit=speedglm.wfit(ky,km,family="binomial") } ##################################### make attributes table, with coefficients from Rcc logistic regression ######### if (Run == "attributes") { #create table of OR, # casese, etc.\ library(data.table) library(glmnet) library(survival) options("width"=420) fit.cv.odd=readRDS("../rv10c.15/fit.cv.rv10c.15.1.99.RDS") #for rv8m #fit.cv.odd=readRDS("../rv10o/fit.cox.rv10o.1.99.RDS") coef.odd=coef.glmnet(fit.cv.odd,s=fit.cv.odd$lambda.1se) fit.cv.even=readRDS("../rv10c.15/fit.cv.rv10c.15.0.98.RDS") #fit.cv.even=readRDS("../rv10o/fit.cox.rv10o.0.98.RDS") coef.even=coef.glmnet(fit.cv.even,s=fit.cv.odd$lambda.1se) coef=(coef.even+coef.odd)/2 #coef=rbind(0.,coef) # add placeholder for intercept so we can use in same way as glm coefficients. co=array(coef[21:1548],dim=c(191,8)) dimnames(co)[[1]]=gsub(".a","",rownames(coef)[21:211],fixed=T) dimnames(co)[[2]]=c(".12-.37y",".37-.62y",".62-1.12y","1.12-1.62y","1.62-2.62y","2.62-3.62y","3.62-5.62y","5.62-7.62y") de=array(rep(coef[2:20]/8,8),dim=c(19,8)) dimnames(de)[[1]]=rownames(coef)[2:20] dimnames(de)[[2]]=c(".12-.37y",".37-.62y",".62-1.12y","1.12-1.62y","1.62-2.62y","2.62-3.62y","3.62-5.62y","5.62-7.62y") c=rbind(de,co) rs=rowSums(c) load("rv10c.17.1.99.RData") #load("data.rv10n.even.RData" y=yss d=cbind(y,m) load("rv10n.0.98.RData") y=yss d2=cbind(y,m) rm(m) d=rbind(d,d2) rm(d2) rownames(co)[57]="PHQ2.not.administered" rownames(co)[61]="PHQ9.not.administered" rownames(co)[87]="Sleepx" names(rs)[76]="PHQ2.not.administered" names(rs)[80]="PHQ9.not.administered" names(rs)[106]="Sleepx" colnames(d)[c(109,300,491,682,873,1064,1255,1446)]="Sleepx" #fixing mistaken label so doesnt match twice. thing="SuicideIdeation";thing="Anxiolytics_Rx";thing="SAE_Falls" Cout=matrix(data=0.,nrow=dim(co)[1],ncol=9) for (i in seq(1,dim(co)[1])){ thing=rownames(co)[i] cols=grep(paste(c("^",thing),collapse=""),colnames(d)) #cc=apply(d[,..cols],1,max) if (length(grep("dminist",thing))==0) {cc=apply(d[,..cols],1,max)} else {cc=apply(d[,..cols],1,min)} df=data.table(d$OC,cc) t=table(df$V1,df$cc) cx=rs[grep(paste(c("^",thing),collapse=""),names(rs))] pD=100*t[2,2]/(t[2,2]+t[1,2]) pSA=100*t[3,2]/(t[3,2]+t[1,2]) pOD=100*t[4,2]/(t[4,2]+t[1,2]) pC=100*sum(t[2:4,2])/(sum(t[2:4,2])+t[1,2]) pCC=100*sum(t[2:4,2])/(sum(t[2:4,2])+t[1,2]) / (100*sum(t[2:4,1])/(sum(t[2:4,1])+t[1,1]) ) out=c(cx,exp(cx),sum(t[2:4,2]),t[1,2],pC,pCC,pD,pSA,pOD) names(out)=c("Coef","Adl.OR","N.cases","N.cntrls","%cases","%cases/%controls","%S.ndi","%SA","%OD") print(round(out,2)) Cout[i,]=out } colnames(Cout)=names(out) rownames(Cout)=rownames(co) write.csv(Cout,file="States.rv10c.17.long.csv") Cout=matrix(data=0.,nrow=19,ncol=9) for (i in seq(1,19)){ thing=colnames(d)[(i+3)] #if (i == 1) {thing="All";cx=0;t=table(d$OC);s=sum(t);pD=100*t[2]/s;pSA=100*t[3]/s;pOD=100*t[4]/s;pC=pD+pSA+pOD} cx=rs[grep(paste(c("^",thing),collapse=""),names(rs))] j=i+3 df=data.table(d$OC,d[,..j]) ;colnames(df)[2]="V2" t=table(df$V1,df$V2) pD=100*t[2,2]/(t[2,2]+t[1,2]) pSA=100*t[3,2]/(t[3,2]+t[1,2]) pOD=100*t[4,2]/(t[4,2]+t[1,2]) pC=100*sum(t[2:4,2])/(sum(t[2:4,2])+t[1,2]) pCC=100*sum(t[2:4,2])/(sum(t[2:4,2])+t[1,2]) / (100*sum(t[2:4,1])/(sum(t[2:4,1])+t[1,1])) out=c(cx,exp(cx),sum(t[2:4,2]),t[1,2],pC,pCC,pD,pSA,pOD) names(out)=c("Coef","Adl.OR","N.cases","N.cntrls","%cases","%cases/%controls","%S.ndi","%SA","%OD") print(round(out,2)) Cout[i,]=out } colnames(Cout)=names(out) rownames(Cout)=colnames(d)[4:22] write.csv(Cout,file="States.rv10c.17.Xsec.csv") } ###################################### look at coefficients from glm and cox models ############################# if (Run == "coeff"){ #look at glm coefficients for intervals library(data.table) library(glmnet) library(survival) options("width"=420) pop1=read.csv("States.rv10c.17.Xsec.csv") pop2=read.csv("States.rv10c.17.long.csv") pop=rbind(pop1,pop2) wd=as.vector(pop$N.cases+pop$N.cntrls) fit.cv.odd=readRDS("../rv10c.17/fit.cv.rv10c.17.1.99.RDS") #fit.cv.odd=readRDS("../rv10o/fit.cox.rv10o.1.99.RDS") #for rv8m coef.odd=coef.glmnet(fit.cv.odd,s=fit.cv.odd$lambda.1se) fit.cv.even=readRDS("../rv10c.17/fit.cv.rv10c.17.0.98.RDS") #fit.cv.even=readRDS("../rv10o/fit.cox.rv10o.0.98.RDS") #for rv8m coef.even=coef.glmnet(fit.cv.even,s=fit.cv.odd$lambda.1se) coef=(coef.even+coef.odd)/2 #coef=rbind(0.,coef) # add placeholder for intercept so we can use in same way as glm coefficients. co=array(coef[21:1548],dim=c(191,8)) dimnames(co)[[1]]=gsub(".a","",rownames(coef)[21:211],fixed=T) dimnames(co)[[2]]=c(".12-.37y",".37-.62y",".62-1.12y","1.12-1.62y","1.62-2.62y","2.62-3.62y","3.62-5.62y","5.62-7.62y") de=array(rep(coef[2:20]/8,8),dim=c(19,8)) dimnames(de)[[1]]=rownames(coef)[2:20] dimnames(de)[[2]]=c(".12-.37y",".37-.62y",".62-1.12y","1.12-1.62y","1.62-2.62y","2.62-3.62y","3.62-5.62y","5.62-7.62y") c=rbind(de,co) rs=rowSums(c) par(mar=c(7,4,2,2)+4.5) barplot(t(c[1:75,]),las=2,main="RV10c.17; Jan 2017 cohort, glm model, all cohorts, census, vitals, labs") barplot(t(c[76:142,]),las=2,main="RV10c.17; Jan 2017 cohort, all cohorts, surveys, Dx" ,ylim=c(-1,2)) barplot(t(c[143:210,]),las=2,main="RV10c.17; Jan 2017 cohort, glm model, all cohorts, Health, medications, BH services") barplot(t(c[1:75,]),wd[1:75],las=2,main="RV10c.17; Jan 2017 cohort, glm model, all cohorts, census, vitals, labs") barplot(t(c[76:142,]),wd[76:142],las=2,main="RV10c.17; Jan 2017 cohort, glm model, all cohorts, surveys, Dx" ,ylim=c(-1,2)) #w=c(77:79,81:142);barplot(t(c[w,]),wd[w],las=2,main="RV10n; case-control event driven model, all cohorts, surveys, Dx" ,ylim=c(-1,2)) barplot(t(c[143:210,]),wd[143:210],las=2,main="RV10c.17; Jan 2017 cohort, glm model, all cohorts, Health, medications, BH services") #look at glm coefficients for cox model # rv8m fit.cv=readRDS("../rv8m/fit.cv.rv8m.2.100.RDS") cox.cv.odd=readRDS("../rv10m/fit.cox.rv10m.1.99.RDS") #for rv8m coef.odd=coef.glmnet(cox.cv.odd,s=cox.cv.odd$lambda.1se) cox.cv.even=readRDS("../rv10m/fit.cox.rv10m.0.98.RDS") #for rv8m coef.even=coef.glmnet(cox.cv.even,s=cox.cv.even$lambda.1se) coef=(coef.even+coef.odd)/2 co=array(coef[20:1547],dim=c(191,8)) #for cox dimnames(co)[[1]]=gsub(".a","",rownames(coef)[20:210],fixed=T) dimnames(co)[[2]]=c(".12-.37y",".37-.62y",".62-1.12y","1.12-1.62y","1.62-2.62y","2.62-3.62y","3.62-5.62y","5.62-7.62y") de=array(rep(coef[1:19]/8,8),dim=c(19,8)) #for cox dimnames(de)[[1]]=rownames(coef)[1:19] #for cox dimnames(de)[[2]]=c(".12-.37y",".37-.62y",".62-1.12y","1.12-1.62y","1.62-2.62y","2.62-3.62y","3.62-5.62y","5.62-7.62y") c=rbind(de,co) rs=rowSums(c) par(mar=c(7,4,2,2)+4.5) barplot(t(c[1:75,]),las=2,main="RV8v.17; cox model for visit patients, 2017, census, vitals, and labs, all cohorts") barplot(t(c[76:142,]),las=2,main="RV8v.17; cox model for visit patients, 2017, all cohorts, surveys, Dx",ylim=c(-1,2)) barplot(t(c[143:210,]),las=2,main="RV8v.17; cox model for visit patients, 2017, all cohorts, Health, medications, BH services") #gindow=c(0.12,0.37,0.62,1.12,1.62,2.62,3.62,5.67,7.62) #window=c(0.25,0.25,0.5,0.5,1,1,2,2) y=ye$OC; y[y==1]=0;y[y==2]=0; y[y==3]=1 # choose sub-outcome pr=predict(fit.cv,newx=m,s=fit.cv$lambda.1se)#,type="response") nD=length(which(y==1)) nL=length(which(y==0)) histL=hist(pr[which(y==0)],br=seq(-15,31,.1)) histD=hist(pr[which(y==1)],br=seq(-15,31,.1)) plot(histL$mids,(1-cumsum(histL$counts)/nL))#,ylim=c(0,1)) points(histD$mids,1-cumsum(histD$counts)/nD,col='red') plot(1-cumsum(histL$counts)/nL,1-cumsum(histD$counts)/nD) abline(c(0,0),c(1,1)) L=1-cumsum(histL$counts)/nL D=1-cumsum(histD$counts)/nD s=smooth.spline(L,D) t=seq(.005,.995,.01) p=predict(s,t) print(c(" AUC: ",sum(p$y))) } ############ Combine models to maximize NDI deaths predicted in C17 cohort library(data.table) # library(glmnet) # library(survival) options("width"=420) xx=load("rv10c.17/rv10c.17.1.99.RData") ysso=yss; mo=m xx=load("rv10c.17/rv10c.17.0.98.RData") yss=rbind(yss,ysso) ; m=rbind(m,mo) rm(mo) n=m[,c(91:182,282:373,473:564,664:755,855:946,1046:1137,1237:1328,1428:1519)] usage=rowSums(n) nn1=read.csv("rv10c.17/neuralnetwork_RV10c.17_scores.csv") nn2=read.csv("rv10c.17/randomforest_RV10c.17_scores.csv") nn3=read.csv("rv10c.17/tabnet_RV10c.17_scores.csv") Ptw3=match(rownames(m),nn3$PatientICN) score=nn3$pred_score[Ptw3] score[score<0.0001] = 0.0001 score[score> 0.9999] = 0.9999 pr3=log(score/(1-score)) cox.cv.odd=readRDS("./rv10c.15/fit.cv.rv10c.15.2y.odd.RDS") #for rv8m coef.odd=coef.glmnet(cox.cv.odd,s=cox.cv.odd$lambda.1se) cox.cv.even=readRDS("./rv10c.15/fit.cv.rv10c.15.2y.even.RDS") #for rv8m coef.even=coef.glmnet(cox.cv.even,s=cox.cv.even$lambda.1se) coefC15All=(coef.even+coef.odd)/2 c=coefC15All[2:1548] # c=coefMH[1:1547] prAll.c15=m%*%c +coefC15All[1] # ; prMH=m %*% c df=data.table(m[,1:5],pr1,pr2,pr3,prRcc,prMH,prV,prAll.c15,prNDI.c15,prNDI.c17,usage,yss) colnames(df)[6:14]=c("prNNRcc","prRFRcc","prTNRcc","prGLMRcc","prCoxMH","prCoxV","prAll.c15","prNDI.c15","prNDI.c17") w=grep("chronic",colnames(m)) #; w=grep("deation",colnames(m)) ; w=grep("ttempt",colnames(m)) ; w=grep("Overdose",colnames(m)) ; x=apply(m[,w],1,sum) df$Nchronic=x #; df$NIdeation=x ; df$NSAttempt=x ; df$NOD=x df$PHQ2.eq.5.6=0 ; col=grep("PHQ2.eq.5.6",colnames(m)) ; w=which(rowSums(m[,col]) >0) ;df$PHQ2.eq.5.6[w]=1 df$PoC_ED =0 ; col=grep("PoC_ED",colnames(m)) ; w=which(rowSums(m[,col]) >0) ;df$PoC_ED[w]=1 df$SAE_Falls =0 ; col=grep("SAE_Falls",colnames(m)) ; w=which(rowSums(m[,col]) >0) ;df$SAE_Falls[w]=1 df$SuicideIdeation =0 ; col=grep("SuicideIdeation",colnames(m)) ;w=which(rowSums(m[,col]) >0) ;df$SuicideIdeation[w]=1 df$Anxiolytics =0 ; col=grep("Anxiolytics",colnames(m)) ; w=which(rowSums(m[,col]) >0) ;df$Anxiolytics[w]=1 df$OpioidForPain =0 ;col=grep("OpioidForPain",colnames(m)) ; w=which(rowSums(m[,col]) >0) ;df$OpioidForPain[w]=1 df$MDD =0 ;col=grep("MDD",colnames(m)) ;w=which(rowSums(m[,col]) >0) ;df$MDD[w]=1 df$BIPOLAR =0 ;col=grep("BIPOLAR",colnames(m)) ; w=which(rowSums(m[,col]) >0) ;df$BIPOLAR[w]=1 df$Schiz =0 ;col=grep("Schiz",colnames(m)) ;w=which(rowSums(m[,col]) >0) ;df$Schiz[w]=1 df$PHQ2.NA =0 ;col=grep("PHQ2.not.a",colnames(m)) ;w=which(rowSums(m[,col]) >0) ;df$PHQ2.NA[w]=1 df$PHQ2.eq.3.4 =0 ;col=grep("PHQ2.eq.3.4",colnames(m)) ;w=which(rowSums(m[,col]) >0) ;df$PHQ2.eq.3.4[w]=1 df$PHQ2.eq.1.2 =0 ;col=grep("PHQ2.eq.1.2",colnames(m)) ;w=which(rowSums(m[,col]) >0) ;df$PHQ2.eq.1.2[w]=1 df$Q9.NA =0 ;col=grep("PHQ9.not",colnames(m)) ;w=which(rowSums(m[,col]) >0) ;df$Q9.NA[w]=1 df$Q9.1 =0 ;col=grep("Q9.score=1",colnames(m)) ;w=which(rowSums(m[,col]) >0) ;df$Q9.1[w]=1 df$Q9.2 =0 ;col=grep("Q9.score=2",colnames(m)) ;w=which(rowSums(m[,col]) >0) ;df$Q9.2[w]=1 df$Q9.3 =0 ;col=grep("Q9.score=3",colnames(m)) ;w=which(rowSums(m[,col]) >0) ;df$Q9.3[w]=1 df$Q9.4 =0 ;col=grep("Q9.score=4",colnames(m)) ;w=which(rowSums(m[,col]) >0) ;df$Q9.4[w]=1 saveRDS(df,file="/mnt/mvpfs/VA_MVP011_CDW/mcmahonb/rv10c.17/AllScores.RDS") df=readRDS("AllScores.RDS") df[Nchronic>1,"Nchronic"]=1 ; df[NIdeation>1,"NIdeation"]=1 ;df[NSAttempt>1,"NSAttempt"]=1 ;df[NOD>1,"NOD"]=1 mm=as.matrix(data.table(df[,1:6],df[,19:22])) # reopt RccNN mm=as.matrix(data.table(df[,1:5],df[,7],df[,19:22])) # reopt RccRF mm=as.matrix(data.table(df[,1:12],df[,15],df[,19:22])) #opt ensemble for combined outcome mm=as.matrix(data.table(df[,1:11],df[,13],df[,15],df[,19:22])) #opt ensemble for NDI mm=as.matrix(data.table(df[,1:5],df[,13],df[,15],df[,19:22])) #opt C15 for NDI mm=as.matrix(data.table(df[,1:5],df[,13],df[,15],df[,19:22])) #opt C15 for All w1=seq(1,544949); w2=seq(544950,1089899) m1=mm[w1,]; m2=mm[w2,] y1=df[w1,OC]; y2=df[w2,OC] y1[y1>1]=1 ; y2[y2 > 1]=1 ww1=1+y1*1; ww2=1+y2*1 fit1=glm(y1~m1,weights=ww1,family='binomial') ;fit2=glm(y2~m2,weights=ww2,family='binomial') l=length(coef(fit1)); c1=coef(fit1)[2:l]; c2=coef(fit2)[2:l] pr1=m1 %*% c2 +coef(fit2)[1] ; pr2=m2 %*% c1 +coef(fit1)[1] y=c(y1,y2) ; pr=c(pr1,pr2) # df$RccNNndi=pr df$RccRFAll=pr df$EnsAll=pr df$EnsNDI=pr df$c15NDI=pr df$c15All=pr dff=data.table(df,y,pr) setkey(dff,pr) t=tail(dff,1090);table(t$OC) t=tail(dff,10900);table(t$OC) t=tail(dff,109000);table(t$OC) plot(df[OC==0,RccRFAll],df[OC==0,RccNNndi],cex=.3) ; points(df[OC==2,RccRFAll],df[OC==2,RccNNndi],cex=.3,col='green') ; points(df[OC==3,RccRFAll],df[OC==3,RccNNndi],cex=.3,col='blue') ; points(df[OC==1,RccRFAll],df[OC==1,RccNNndi],cex=1,col='red',pch=20) plot(df[OC==0,EnsAll],df[OC==0,EnsNDI],cex=.3) ; points(df[OC==2,EnsAll],df[OC==2,EnsNDI],cex=.3,col='green') ; points(df[OC==3,EnsAll],df[OC==3,EnsNDI],cex=.3,col='blue') ; points(df[OC==1,EnsAll],df[OC==1,EnsNDI],cex=1,col='red',pch=20) plot(df[OC==0,prNDI.c17],df[OC==0,c15NDI],cex=.3) ; points(df[OC==2,prNDI.c17],df[OC==2,c15NDI],cex=.3,col='green') ; points(df[OC==3,prNDI.c17],df[OC==3,c15NDI],cex=.3,col='blue') ; points(df[OC==1,prNDI.c17],df[OC==1,c15NDI],cex=1,col='red',pch=20) ######################################## REOPTIMIZE ########################## if (Run == "reoptimize"){ library(data.table) library(glmnet) library(survival) options("width"=420) # read in tabnet scores and align with order of patients from input data, assumped to be in ps #xx=load("../rv10c.15/rv10c.15.1.99.RData")#;load("rv10o.1.99.RData") load("rv10m.1.99.RData") xx=load("../rv10c.17/rv10c.17.1.99.RData")#;load("rv10o.1.99.RData") load("rv10m.1.99.RData") yss[OCDate>19,"OC"]=0 ; y=yss$OC;y[y>1]= 0 #ICNb=read.csv("MACE6b.even.ICN.csv") # ;ICNb=read.csv("MACE6b.odd.ICN.csv") nn=data.table(read.csv("../rv10c.17/neuralnetwork_RV10c.17_scores.csv",stringsAsFactors=F)) #neuralnetwork_rv9n_scores.csv",stringsAsFactors=F,colClasses=c("character","integer","numeric","integer"))) nn=data.table(read.csv("../rv10c.17/neuralnetwork_RV10c.17_scores.csv",stringsAsFactors=F)) #training for ML models is on even cohorts, so test against odd cohorts. nn=data.table(read.csv("tabnet_RV10m_2_scores.csv",stringsAsFactors=F)) nn=data.table(read.csv("randomforest_RV10m_2_scores.csv",stringsAsFactors=F)) nn=data.table(read.csv("../rv10c.17/neuralnetwork_RV10c.17_2_scores.csv",stringsAsFactors=F)) Ptw=match(rownames(m),nn$PatientICN) #which(is.na(Ptw)==T) fit.cv=readRDS("../rv10n/fit.cv.rv10n.0.98.RDS") lscores=predict(fit.cv,newx=m,s=fit.cv$lambda.1se) softmax=readRDS("rv8n.nn.rv8.17scores.RDS") scores=softmax$score_2-softmax$score_0 #ps=ICNb[,1] #Ptw=match(rownames(m),nn$PatientICN) #Ptw=Ptw[is.na(Ptw)==F]=Ptw[is.na(Ptw)==F] score=nn$pred_score[Ptw] score[score<0.0001] = 0.0001 score[score> 0.9999] = 0.9999 #yy=y$OC[Ptw] pr=log(score/(1-score)) pr=log(scores[Ptw,4]/(1-scores[Ptw,4])) + km[,4]*2.5 + (km[,214]+km[,458]+km[,702]+km[,946])*0.28 + (km[,172]+km[,416]+km[,660]+km[,984])*0.23 + (km[,154]+km[,398]+km[,642]+km[,886])*0.15 # tabnet, LDL-5yr # Pain.NA is m[,39],m[,230],m[,421],m[,612],m[,803],m[,994],m[,1185],m[,1376] # SI is m[,90],m[,281],m[,472],m[,663],m[,854],m[,1045],m[,1236],m[,1427] # SA =s m[,89],m[,280],m[,471],m[,662],m[,853],m[,1044],m[,1235],m[,1426] # OD is m[,97],m[,288],m[,479],m[,670],m[,861],m[,1052],m[,1243],m[,1434] data=data.table(lscores,m[,1:5],m[,39],m[,230],m[,421],m[,612],m[,803],m[,994],m[,1185],m[,1376],m[,90],m[,281],m[,472],m[,663],m[,854],m[,1045],m[,1236],m[,1427],m[,89],m[,280],m[,471],m[,662],m[,853],m[,1044],m[,1235],m[,1426],m[,97],m[,288],m[,479],m[,670],m[,861],m[,1052],m[,1243],m[,1434]) #data=data.table(log(scores[Ptw,4]/(1-scores[Ptw,4])),km[,1:5],km[,181],km[,425],km[,669],km[,913],km[,1157],km[,1401],km[,1645],km[,1889],km[,191],km[,435],km[,679],km[,923],km[,1167],km[,1411],km[,1655],km[,1899], km[,199],km[,443],km[,687],km[,931],km[,1175],km[,1419],km[,1663],km[,1907], km[,214],km[,458],km[,702],km[,946],km[,1190],km[,1434],km[,1678],km[,1922] ,km[,28],km[,272],km[,516],km[,760],km[,1004],km[,1248],km[,1492],km[,1736] ,km[,174],km[,418],km[,662],km[,906],km[,1150],km[,1394],km[,1638],km[,1882] )#, km[,57],km[,301],km[,545],km[,789]) d=as.matrix(data) cv.fit=cv.glmnet(d,y,family="binomial") #fit=glm(y~km,family="binomial") pr=predict(cv.fit,d,s=cv.fit$lambda.min) #pr=log(scores[Ptw,4]/(1-scores[Ptw,4])) + km[,4]*2.5 + (km[,214]+km[,458]+km[,702]+km[,946])*0.28 + (km[,172]+km[,416]+km[,660]+km[,984])*0.23 + (km[,154]+km[,398]+km[,642]+km[,886])*0.15 # tabnet, LDL-5yr coef=coef.glmnet(fit.cv,s=fit.cv$lambda.1se) c=coef[2:1958] pr=km%*%c +coef[1] pr=km%*%c - km[,4]*1. - (km[,214]+km[,458]+km[,702]+km[,946])*0.32 pr=predict(fit.cv,newx=m,s=fit.cv$lambda.1se)#,type="response") nD=length(which(y==1)) nL=length(which(y==0)) histL=hist(pr[which(y==0)],br=seq(-18,15,.1)) histD=hist(pr[which(y==1)],br=seq(-18,15,.1)) plot(histL$mids,(1-cumsum(histL$counts)/nL))#,ylim=c(0,1)) points(histD$mids,1-cumsum(histD$counts)/nD,col='red') plot(1-cumsum(histL$counts)/nL,1-cumsum(histD$counts)/nD) abline(c(0,0),c(1,1)) L=1-cumsum(histL$counts)/nL D=1-cumsum(histD$counts)/nD s=smooth.spline(L,D) t=seq(.005,.995,.01) p=predict(s,t) print(c(" AUC: ",sum(p$y))) } # library(ROCR) pred=prediction(pr,y) pp=performance(pred,"auc") ; auc=pp@y.values ;auc pp=performance(pred,"aucpr") ; aucpr=pp@y.values; aucpr main="c15 NDI-trained, c17 reopt" ; pr=df$c15NDI ; y=df$OC ; br=0.2; y[y>1]=0 #C15GLMndi.hist.lin.png main="Rcc-NNet Comb OC trained, c17 reopt"; pr=df$RccNNndi ; y=df$OC ; br=0.2; y[y>1]=0 #RccNNndi.hist.lin.png main="Ensemble Comb OC trained, c17 reopt"; pr=df$EnsNDI ; y=df$OC ; br=0.2; y[y>1]=0 #EnsNDI.hist.lin.png main="c15 GLM Comb OC-trained, c17 reopt" ; pr=df$c15All ; y=df$OC ; br=0.1; y[y>1]=1 #C15GLMAll.hist.lin.png main="Rcc RF Comb OC-trained, c17 reopt" ; pr=df$RccRFAll ; y=df$OC ; br=0.1; y[y>1]=1 #RccRFAll.hist.lin.png main="Ensemble Comb OC-trained, c17 reopt" ; pr=df$EnsAll ; y=df$OC ; br=0.1; y[y>1]=1 #EnsAll.hist.lin.png # pr=df$RccNNndi pr=df$RccRFAll pr=df$EnsAll pr=df$EnsNDI pr=df$c15NDI ## prNNRcc prRFRcc prTNRcc prGLMRcc prCoxMH prCoxV prAll.c15 pr=df$prNDI.c15 pr=df$prNDI.c17 yss=data.table(df$Tpr,df$OCDate,df$OC) ; colnames(yss)=c("Tpr","OCDate","OC") hNDI=hist(pr[which(yss$OC==1)],br=seq(-15,11,.1),plot=T) #make histogram of suboutcomes hSA=hist(pr[which(yss$OC==2)],br=seq(-15,11,.1),plot=T) hOD=hist(pr[which(yss$OC==3)],br=seq(-15,11,.1),plot=T) hl=hist(pr[which(y!=1)],br=seq(-15,11,.1),plot=T) hd=hist(pr[which(y==1)],br=seq(-15,11,.1),plot=T) plot(hl$mids,hl$density,col='blue',xlim=c(-7,3),log='y',pch=20,ylim=c(0.001,0.8),xlab="log(Pr/(1-Pr))",ylab="normalized counts",main=main,cex=2) points(hd$mids,hd$density,col='red',xlim=c(-7,4),pch=20,cex=2) lines(hNDI$mids,hNDI$density,col='magenta',xlim=c(-7,4),pch=20,lwd=3) lines(hSA$mids,hSA$density,col='cyan',xlim=c(-7,4),pch=20,lwd=3) lines(hOD$mids,hOD$density,col='black',xlim=c(-7,4),pch=20,lwd=3) legend("topright",legend=c("Control","Combined Outcome","NDI suicide","Suicide Attempt","Overdose"),text.col=c("blue","red","magenta","cyan","black"),cex=2) lines(hl$mids,(1-cumsum(hd$counts+hl$counts)/(nL+nD)),col='black',lty=2) abline(v=3.2) # abline(v=2.35) top 1% abline(v=4.6) # abline(v=4.95) top 0.1% lines(c(2,6),2*c(.15,.0003),lwd=3,col='red') lines(c(2,6),c(.15,.0003),lwd=3,col='magenta') lines(c(2,6),c(.15/5,.0003/20),lwd=3,col='blue') #for nn on C-17 lines(c(-1,2),c(.05,.006),lwd=3,col='magenta') lines(c(-1,3.8),c(.12/13,.0003/10),lwd=3,col='blue') lines(c(-2,6),2*c(.12,.0003),lwd=3,col='red') abline(v=-0.8) abline(v=1.9) xx=load("refit.visit.c17.even.RData") > pre=pr;rm(m);ysse=yss;de=d > xx=load("refit.MH.c17.odd.RData") > xx=load("refit.visit.c17.odd.RData") > yss=rbind(yss,ysse);rm(m);d=rbind(d,de);pr=rbind(pr,pre);df=data.table(yss,pr,d);colnames(df)[5]='s2';setkey(df,s1) > table(df[1088809:1089899,"OC"]) 0 1 2 3 371 7 646 67 > table(df[1079000:1089899,"OC"]) 0 1 2 3 7166 96 2827 811 #subgroup analysis MACE6b if (Run == "subgroup"){ library(data.table) br=.1 df=readRDS("scores.RDS") # prr=df$RccNNndi prr=df$RccRFAll prr=df$EnsAll prr=df$EnsNDI prr=df$c15NDI prr=df$c15All ## prNNRcc prRFRcc prTNRcc prGLMRcc prCoxMH prCoxV prr=df$prAll.c15 pr=df$prNDI.c15 prr=df$prNDI.c17 # Nchronic NIdeation NSAttempt NOD PHQ2.eq.5.6 PoC_ED SAE_Falls SuicideIdeation Anxiolytics OpioidForPain MDD BIPOLAR Schiz PHQ2.NA PHQ2.eq.3.4 PHQ2.eq.1.2 Q9.NA Q9.1 Q9.2 Q9.3 Q9.4 s=seq(1,240)/20-8 ylab="Probability in C-17 cohort" main="c15 NDI-trained, c17 reopt" ; prr=df$c15NDI ; y=df$OC ; br=0.2; y[y>1]=0 main="Rcc-NNet Comb OC trained, c17 reopt"; prr=df$RccNNndi ; y=df$OC ; br=0.2; y[y>1]=0 main="Ensemble Comb OC trained, c17 reopt"; prr=df$EnsNDI ; y=df$OC ; br=0.2; y[y>1]=0 main="c15 GLM Comb OC-trained, c17 reopt" ; prr=df$c15All ; y=df$OC ; br=0.1; y[y>1]=1 main="Rcc RF Comb OC-trained, c17 reopt" ; prr=df$RccRFAll ; y=df$OC ; br=0.1; y[y>1]=1 main="Ensemble Comb OC-trained, c17 reopt" ; prr=df$EnsAll ; y=df$OC ; br=0.1; y[y>1]=1 # ; off=log10(table(y)[1]/table(y)[2]) w=seq(1,dim(df)[1]) ;pr=prr ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) plot(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,xlim=c(-6,4.5),ylim=c(0.001,1),log='y',main=main,xlab="score",ylab=ylab,cex.axis=1.5) w=which(df$PHQ2.eq.5.6 >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='red') #PHQ2.eq.5.6 w=which(df$PoC_ED >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='green') #PoC_ED w=which(df$SAE_Falls >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='blue') #SAE_Falls w=which(df$NIdeation >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='cyan') #SuicideIdeation lines(s+.65,exp(s)/(1+exp(s)),col='black',lw=2) legend("topleft",legend=c("All","PHQ2.eq.5.6","PoC_ED","SAE_Falls","SuicideIdeation"),text.col=c("black","red","green","blue","cyan"),cex=1.5) nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,.1),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,.1),plot=F) w=seq(1,dim(df)[1]) ;pr=prr ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) plot(histD$mids,histD$density,col="black",type='b',pch=20,lwd=4,main=main,xlim=c(-6,4.5),ylim=c(0,.6),xlab="score",ylab=ylab,cex.axis=1.5) w=which(df$PHQ2.eq.5.6 >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="red",type='b',pch=20,lwd=4) w=which(df$PoC_ED >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="green",type='b',pch=20,lwd=4) w=which(df$SAE_Falls >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="blue",type='b',pch=20,lwd=4) w=which(df$NIdeation >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="cyan",type='b',pch=20,lwd=4) legend("topright",legend=c("All","PHQ2.eq.5.6","PoC_ED","SAE_Falls","SuicideIdeation"),text.col=c("black","red","green","blue","cyan"),cex=1.5) main="c15 NDI-trained, c17 reopt" ; prr=df$c15NDI ; y=df$OC ; br=0.2; y[y>1]=0 #C15GLMndi.demo.calib.png main="Rcc-NNet Comb OC trained, c17 reopt"; prr=df$RccNNndi ; y=df$OC ; br=0.2; y[y>1]=0 #RccNNndi.demo.calib.png main="Ensemble Comb OC trained, c17 reopt"; prr=df$EnsNDI ; y=df$OC ; br=0.2; y[y>1]=0 #EnsNDI.demo.calib.png main="c15 GLM Comb OC-trained, c17 reopt" ; prr=df$c15All ; y=df$OC ; br=0.1; y[y>1]=1 #C15GLMAll.demo.calib.png main="Rcc RF Comb OC-trained, c17 reopt" ; prr=df$RccRFAll ; y=df$OC ; br=0.1; y[y>1]=1 #RccRFAll.demo.calib.png main="Ensemble Comb OC-trained, c17 reopt" ; prr=df$EnsAll ; y=df$OC ; br=0.1; y[y>1]=1 #EnsAll.demo.calib.png #by demographics w=which(df$Female ==0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) plot(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,xlim=c(-6,4.5),ylim=c(0.001,1),log='y',main=main,xlab="score",ylab=ylab,cex.axis=1.5) w=which(df$Female >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='red') # w=which(df$Rac >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='green') # w=which(df$Eth >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='blue') # w=which(df$Young >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='cyan') # w=which(df$Old >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='magenta') # lines(s+.65,exp(s)/(1+exp(s)),col='black',lw=2) legend("topleft",legend=c("Male","Female","Black","Hispanic","Young","Old"),text.col=c("black","red","green","blue","cyan","magenta"),cex=1.5) nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,.1),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,.1),plot=F) w=which(df$Female ==0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) plot(histD$mids,histD$density,col="black",type='b',pch=20,lwd=4,main=main,xlim=c(-6,4.5),ylim=c(0,.6),xlab="score",ylab=ylab,cex.axis=1.5) w=which(df$Female >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="red",type='b',pch=20,lwd=4) w=which(df$Rac >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="green",type='b',pch=20,lwd=4) w=which(df$Eth >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="blue",type='b',pch=20,lwd=4) w=which(df$Young >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="cyan",type='b',pch=20,lwd=4) w=which(df$Old >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="magenta",type='b',pch=20,lwd=4) legend("topright",legend=c("Male","Female","Black","Hispanic","Young","Old"),text.col=c("black","red","green","blue","cyan","magenta"),cex=1.5) main="c15 NDI-trained, c17 reopt" ; prr=df$c15NDI ; y=df$OC ; br=0.2; y[y>1]=0 #C15GLMndi.MH.xx.png main="Rcc-NNet Comb OC trained, c17 reopt"; prr=df$RccNNndi ; y=df$OC ; br=0.2; y[y>1]=0 #RccNNndi.MH.xx.png main="Ensemble Comb OC trained, c17 reopt"; prr=df$EnsNDI ; y=df$OC ; br=0.2; y[y>1]=0 #EnsNDI.MH.xx.png main="c15 GLM Comb OC-trained, c17 reopt" ; prr=df$c15All ; y=df$OC ; br=0.1; y[y>1]=1 #C15GLMAll.MH.xx.png main="Rcc RF Comb OC-trained, c17 reopt" ; prr=df$RccRFAll ; y=df$OC ; br=0.1; y[y>1]=1 #RccRFAll.MH.xx.png main="Ensemble Comb OC-trained, c17 reopt" ; prr=df$EnsAll ; y=df$OC ; br=0.1; y[y>1]=1 #EnsAll.MH.xx.png ##by mental w=seq(1,dim(df)[1]) ;pr=prr ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) plot(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,xlim=c(-6,4.5),ylim=c(0.001,1),log='y',main=main,xlab="score",ylab=ylab,cex.axis=1.5) w=which(df$Anxiolytics >0 ) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='red') # w=which(df$OpioidForPain >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='green') # w=which(df$MDD >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='blue') # w=which(df$BIPOLAR >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='cyan') # w=which(df$Schiz >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='magenta') # lines(s+.65,exp(s)/(1+exp(s)),col='black',lw=2) legend("topleft",legend=c("All","Anxiolytics","OpioidForPain","MDD","Bipolar","Schizophrenia"),text.col=c("black","red","green","blue","cyan","magenta"),cex=1.5) nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,.1),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,.1),plot=F) w=seq(1,dim(df)[1]) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) plot(histD$mids,histD$density,col="black",type='b',pch=20,lwd=4,main=main,xlim=c(-6,4.5),ylim=c(0,.6),xlab="score",ylab=ylab,cex.axis=1.5) w=which(df$Anxiolytics >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="red",type='b',pch=20,lwd=4) w=which(df$OpioidForPain >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="green",type='b',pch=20,lwd=4) w=which(df$MDD >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="blue",type='b',pch=20,lwd=4) w=which(df$BIPOLAR >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="cyan",type='b',pch=20,lwd=4) w=which(df$Schiz >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="magenta",type='b',pch=20,lwd=4) legend("topright",legend=c("All","Anxiolytics","OpioidForPain","MDD","Bipolar","Schizophrenia"),text.col=c("black","red","green","blue","cyan","magenta"),cex=1.5) main="c15 NDI-trained, c17 reopt" ; prr=df$c15NDI ; y=df$OC ; br=0.2; y[y>1]=0 #C15GLMndi.usage.calib.png main="Rcc-NNet Comb OC trained, c17 reopt"; prr=df$RccNNndi ; y=df$OC ; br=0.2; y[y>1]=0 #RccNNndi.usage.calib.png main="Ensemble Comb OC trained, c17 reopt"; prr=df$EnsNDI ; y=df$OC ; br=0.2; y[y>1]=0 #EnsNDI.usage.calib.png main="c15 GLM Comb OC-trained, c17 reopt" ; prr=df$c15All ; y=df$OC ; br=0.1; y[y>1]=1 #C15GLMAll.usage.calib.png main="Rcc RF Comb OC-trained, c17 reopt" ; prr=df$RccRFAll ; y=df$OC ; br=0.1; y[y>1]=1 #RccRFAll.usage.calib.png main="Ensemble Comb OC-trained, c17 reopt" ; prr=df$EnsAll ; y=df$OC ; br=0.1; y[y>1]=1 #EnsAll.usage.calib.png #by Dx code usage w=which(df$usage <2) ;pr=prr ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) plot(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,xlim=c(-6,4.5),ylim=c(0.001,1),log='y',main=main,xlab="score",ylab=ylab,cex.axis=1.5) w=which(df$usage >=2 & df$usage < 50 ) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='red') # w=which(df$usage >=50 & df$usage < 100 ) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='green') # w=which(df$usage >= 100) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='blue') # lines(s+.65,exp(s)/(1+exp(s)),col='black',lw=2) legend("topleft",legend=c("usage < 2","2 <= usage < 50","50 <= usage < 100","usage >= 100"),text.col=c("black","red","green","blue"),cex=1.5) nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,.1),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,.1),plot=F) w=which(df$usage <2) ;pr=prr ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) plot(histD$mids,histD$density,col="black",type='b',pch=20,lwd=4,main=main,xlim=c(-6,4.5),ylim=c(0,.6),xlab="score",ylab=ylab,cex.axis=1.5) w=which(df$usage >=2 & df$usage < 50 ) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="red",type='b',pch=20,lwd=4) w=which(df$usage >=50 & df$usage < 100 ) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="green",type='b',pch=20,lwd=4) w=which(df$usage >= 100) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="blue",type='b',pch=20,lwd=4) legend("topright",legend=c("usage < 2","2 <= usage < 50","50 <= usage < 100","usage >= 100"),text.col=c("black","red","green","blue"),cex=1.5) main="c15 NDI-trained, c17 reopt" ; prr=df$c15NDI ; y=df$OC ; br=0.2; y[y>1]=0 #C15GLMndi.Q9.calib.png main="Rcc-NNet Comb OC trained, c17 reopt"; prr=df$RccNNndi ; y=df$OC ; br=0.2; y[y>1]=0 #RccNNndi.Q9.calib.png main="Ensemble Comb OC trained, c17 reopt"; prr=df$EnsNDI ; y=df$OC ; br=0.2; y[y>1]=0 #EnsNDI.Q9.calib.png main="c15 GLM Comb OC-trained, c17 reopt" ; prr=df$c15All ; y=df$OC ; br=0.1; y[y>1]=1 #C15GLMAll.Q9.calib.png main="Rcc RF Comb OC-trained, c17 reopt" ; prr=df$RccRFAll ; y=df$OC ; br=0.1; y[y>1]=1 #RccRFAll.Q9.calib.png main="Ensemble Comb OC-trained, c17 reopt" ; prr=df$EnsAll ; y=df$OC ; br=0.1; y[y>1]=1 #EnsAll.Q9.calib.png #by PHQ9 Q9 w=which(df$Q9.NA >0 ) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) plot(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,xlim=c(-6,4.5),ylim=c(0.001,1),log='y',main=main,xlab="score",ylab=ylab,cex.axis=1.5) w=which(df$Q9.1 >0 ) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='blue') # w=which(df$Q9.2 >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='green') # w=which(df$Q9.3 >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='cyan') # w=which(df$Q9.4 >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='red') # lines(s+.65,exp(s)/(1+exp(s)),col='black',lw=2) legend("topleft",legend=c("All","Q9 = 1","Q9 = 2","Q9 = 3","Q9 = 4"),text.col=c("black","blue","green","cyan","red"),cex=1.5) nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,.1),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,.1),plot=F) w=which(df$Q9.NA >0 ) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) plot(histD$mids,histD$density,col="black",type='b',pch=20,lwd=4,main=main,xlim=c(-6,4.5),ylim=c(0,.6),xlab="score",ylab=ylab,cex.axis=1.5) w=which(df$Q9.1 >0 ) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="blue",type='b',pch=20,lwd=4) w=which(df$Q9.2 >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="green",type='b',pch=20,lwd=4) w=which(df$Q9.3 >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="cyan",type='b',pch=20,lwd=4) w=which(df$Q9.4 >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="red",type='b',pch=20,lwd=4) legend("topright",legend=c("All","Q9 = 1","Q9 = 2","Q9 = 3","Q9 = 4"),text.col=c("black","blue","green","cyan","red"),cex=1.5) #by PHQ2 w=which(df$PHQ2.NA >0 ) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) plot(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,xlim=c(-6,4.5),ylim=c(0.001,1),log='y',main=main,xlab="score",ylab=ylab,cex.axis=1.5) w=which(df$PHQ2.NA ==0 & df$PHQ2.eq.1.2 ==0 &df$PHQ2.eq.3.4 ==0 &df$PHQ2.eq.5.6 ==0 &) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='red') # w=which(df$PHQ2.eq.1.2 >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='green') # w=which(df$PHQ2.eq.3.4 >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='blue') # w=which(df$PHQ2.eq.5.6 >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) points(histL$mids,histD$counts/(histL$counts+histD$counts),cex=sqrt(histL$counts+histD$counts)/20,pch=20,col='cyan') # lines(s+.65,exp(s)/(1+exp(s)),col='black',lw=2) legend("topleft",legend=c("PHQ2.NA","PHQ2 = 0","PHQ2 = 1-2","PHQ2 = 3-4","PHA2 = 5-6"),text.col=c("black","red","green","blue","cyan"),cex=1.5) nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,.1),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,.1),plot=F) w=which(df$PHQ2.NA >0 ) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) plot(histD$mids,histD$density,col="black",type='b',pch=20,lwd=4,main=main,xlim=c(-6,4.5),ylim=c(0,.6),xlab="score",ylab=ylab,cex.axis=1.5) w=which(df$PHQ2.NA ==0 & df$PHQ2.eq.1.2 ==0 &df$PHQ2.eq.3.4 ==0 &df$PHQ2.eq.5.6 ==0 &) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="red",type='b',pch=20,lwd=4) w=which(df$PHQ2.eq.1.2 >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="green",type='b',pch=20,lwd=4) w=which(df$PHQ2.eq.3.4 >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="blue",type='b',pch=20,lwd=4) w=which(df$PHQ2.eq.5.6 >0) ; pr=prr[w] ;nD=length(which(y[w]==1)) ;nL=length(which(y[w]==0)) ;histL=hist(pr[which(y[w]==0)],br=seq(-10,11,br),plot=F); histD=hist(pr[which(y[w]==1)],br=seq(-10,11,br),plot=F) lines(histD$mids,histD$density,col="cyan",type='b',pch=20,lwd=4) legend("topright",legend=c("PHQ2.NA","PHQ2 = 0","PHQ2 = 1-2","PHQ2 = 3-4","PHA2 = 5-6"),text.col=c("black","red","green","blue","cyan"),cex=1.5) #subgroup analysis MACE6b # histograms with age xx=load("rv10c.17.1.99.RData") ysso=yss; mo=m xx=load("rv10c.17.0.98.RData") ysse=yss; me=m m=rbind(me,mo); yss=rbind(ysse,ysso) d=readRDS("DoB.all.RDS") d$age=17-d$BirthDate d$PatientICN=as.character(d$PatientICN) thing = "PHQ9.Q9.score=4" w=grep(thing,colnames(m)) x=apply(m[,w],1,max) var=data.table(rownames(m),yss$OC,x) colnames(var)=c("PatientICN","OC","PHQ9.Q9.score_4") things=c("PHQ2.eq.5.6", "SuicideIdeation", "Antidepressant_Rx", "PoC_ED", "BIPOLAR", "Anxiolytics_Rx", "OpioidForPain_Rx", "MH_CatchAll", "MHOC_OtherMH","EH_METCANCR","OUD","SAE_Acet","SAE_Falls","SuicideAttempt","OpioidOverdose") for (i in seq(1,15)){ thing=things[i] print(thing) w=grep(thing,colnames(m)) x=apply(m[,w],1,max) var[,thing] = x } new=merge(d,var,by="PatientICN") new[age>100,"age"]=100 h0=hist(new[OC==0,age],br=seq(10,100)) h1=hist(new[OC==1,age],br=seq(10,100)) h2=hist(new[OC==2,age],br=seq(10,100)) h3=hist(new[OC==3,age],br=seq(10,100)) h4=hist(new[PHQ9.Q9.score_4==1,age],br=seq(10,100)) h5=hist(new[PHQ2.eq.5.6==1,age],br=seq(10,100)) h6=hist(new[SuicideIdeation==1,age],br=seq(10,100)) h7=hist(new[Antidepressant_Rx==1,age],br=seq(10,100)) h8=hist(new[PoC_ED==1,age],br=seq(10,100)) h9=hist(new[BIPOLAR==1,age],br=seq(10,100)) h10=hist(new[Anxiolytics_Rx==1,age],br=seq(10,100)) h11=hist(new[OpioidForPain_Rx==1,age],br=seq(10,100)) h12=hist(new[MH_CatchAll==1,age],br=seq(10,100)) h13=hist(new[MHOC_OtherMH==1,age],br=seq(10,100)) h14=hist(new[EH_METCANCR==1,age],br=seq(10,100)) h15=hist(new[OUD==1,age],br=seq(10,100)) h16=hist(new[SAE_Acet==1,age],br=seq(10,100)) h17=hist(new[SAE_Falls==1,age],br=seq(10,100)) h18=hist(new[SuicideAttempt==1,age],br=seq(10,100)) h19=hist(new[OpioidOverdose==1,age],br=seq(10,100)) out=data.table(h0$mids,(h0$counts+h1$counts+h2$counts+h3$counts),h0$counts,h1$counts,h2$counts,h3$counts,h4$counts,h5$counts,h6$counts,h7$counts,h8$counts,h9$counts,h10$counts,h11$counts,h12$counts,h13$counts, h14$counts,h15$counts,h16$counts,h17$counts,h18$counts,h19$counts) colnames(out)=c("age","total_OC","controls","suicide","SA","OD","PHQ9Q9_4","PHQ2.eq.5.6", "SuicideIdeation", "Antidepressant_Rx", "PoC_ED", "BIPOLAR", "Anxiolytics_Rx", "OpioidForPain_Rx", "MH_CatchAll", "MHOC_OtherMH", "EH_METCANCR","OUD","SAE_Acet","SAE_Falls","SuicideAttempt","OpioidOverdose") saveRDS(out,"hist_subOC_subgroup.RDS") colnames(out)=c("age","total_OC","controls","suicide","SA","OD","PHQ9Q9_4","PHQ2.eq.5.6", "SuicideIdeation", "Antidepressant_Rx", "PoC_ED", "BIPOLAR", "Anxiolytics_Rx", "OpioidForPain_Rx", "MH_CatchAll", "MHOC_OtherMH", "EH_METCANCR")