PERMIT AKI algorithm- CPRD ready

ssh mqbpjhr4@incline256.itservices.manchester.ac.uk
module load apps/gcc/R/3.3.3
R

load("crea.rep.rda")

library (taRifx)
library (dplyr)
library(tcltk)
library(lubridate)
require(dplyr)
library(survival)

##################################################################################
#BASED ON THE NHS UK AKI ALGORITHM, OMITTING INDICATORS FROM URINE ALBUMIN
###################################################################################

#SELECT ELIGIBLE PATIENTS
attach(crea.rep)
crea.rep$KDmark<-ifelse(!is.na(crea.rep$MDRDeGFR)&crea.rep$MDRDeGFR<60,1,0)
KD<-crea.rep[crea.rep$KDmark>0,]
a<-crea.rep[crea.rep$PatientID %in% KD$PatientID,]
#The dataset is cleaned again at the end to make sure that CKD patients are not classed as having an AKI
na.omit(a)

a$CKDate<-ifelse(a$CKDStage>0,a$event.date,NA) 
smalltab<-a[!is.na(a$CKDate),c("PatientID","CKDate")]
first<-smalltab %>%
   group_by(PatientID) %>%  
   arrange(CKDate) %>% 
   slice(which.min(as.numeric(CKDate))) %>%
as.data.frame    
head(first)
  
a<-merge(a,as.data.frame(first),all.x=TRUE)
a$event.date<-as.Date(as.character(a$event.date),format="%Y-%m-%d")
a$CKDate<-as.Date(as.character(a$CKDate),format="%Y-%m-%d")

#CKDate is the date at which CKD is first diagnosed

a<-a[!is.na(a$PatientID) & (is.na(a$CKDate)|a$CKDate>a$event.date),]
attach(a)
summary(a$CKDate,na.rm=TRUE)
length(a$PatientID)
#We only have 1559 tests to assess

#index<-which(event.date>=as.Date(ymd("2007-01-01")))
#The index start date should be 1 year prior to the first recorded creatinine test of the study window
#a<-a[index,]
length(a$PatientID) #458030 potentially qualifying tests
levels(as.factor(a$Ethnicity)) 

#CALCULATING MDRD eGFR
o<-ifelse(a$Ethnicity=="4",1.212,1)
o<-ifelse((is.na(a$Ethnicity)|a$Ethnicity==""),1,o)
p<-ifelse(a$Gender=="F",0.742,1)
a$MDRDeGFR<-(175*((a$Creatinine/88.42)^-1.154))*(a$Age^-0.203)*o*p

a<-a[,c("PatientID","event.date","Creatinine","Age", "Gender","MDRDeGFR")]
attach(a)
names(a)
########################################################################### CHECKED

#Order creatinine values and find the max per day per patient
destring(c(PatientID,Creatinine,Age,Gender,event.date,MDRDeGFR))
a<-unique(a)
length(a$PatientID) #Remaining entries
#1526
########################################################################## CHECKED

require(dplyr)
sorted <- a %>% 
          arrange(PatientID,event.date) %>%
          group_by(PatientID) %>%
          mutate(protagonist=row_number())
a<-data.frame(sorted)
attach(a)

#GIVE EACH TEST AN ORDER NUMBER PER PATIENT SO WE CAN LATER LOOP THROUGH THEM
####################################################################################### CHECKED
# PART 2 - LOOP FOR AKI FLAGS 
#LOOP FORMATION FOR COMPARING EACH TEST WITH PREVIOUS TESTS
#CREATE A MARKER VARIABLE FOR EACH TIME AKI EALERT CRITERIA SATISFIED

newAKI<-0
AKIyear<-0
AKIweek<-0
max(protagonist) #number of loop iterations needed.

#LOOP for x values, where x = the highest number of tests a patient had in the index year (max(protagonist)) 
mightydate<-NA
mightyvalue<-NA
protagonistdos<-NA
yearmarker<-0
weekmarker<-0
Dmarker<-0
yearmedian<-0
weekmin<-0
Dmin<-0
yearvalue<-0
weekvalue<-0
Dvalue<-0
D<-0
RV1<-0
RV2<-0
protagonistdos<-0
protagonistvalue<-0

a$ref.low<-ifelse(a$Gender=="M",62,44)
a$ref.up<-ifelse(a$Gender=="M",115,97)
a <-a[order(PatientID, event.date),] 
attach(a)

############################################################################################################# CHECKED
pb<-tkProgressBar(title="Identifying AKI events",min=0,max=max(protagonist),width=300) #indicates progress through the loop
for (i in 1:max(protagonist)) {
setTkProgressBar(pb,i,label=paste(round(i/max(protagonist)*100, 0),"% done"))
mightydate<-as.Date(ifelse(protagonist==protagonist[i] & PatientID==PatientID[i],event.date[i],mightydate),origin="1970-01-01")
protagonistdos<-ave(mightydate,PatientID,FUN=max)
#Apply each test date to all other tests from the same patient to find the time difference

yearmarker[i]<-ifelse(as.Date(protagonistdos)-event.date[i]<=365 & (protagonistdos-event.date[i]>90),1,yearmarker)
weekmarker[i]<-ifelse(as.Date(protagonistdos)-event.date[i]<=7 & (protagonistdos-event.date)[i]>=0,1,weekmarker) 
Dmarker[i]<-ifelse(as.Date(protagonistdos)-event.date[i]<=2 & (protagonistdos-event.date)[i]>=0,1,Dmarker) 

yearvalue[i]<-ifelse(yearmarker[i]==1,Creatinine[i],NA)
weekvalue[i]<-ifelse(weekmarker[i]>0,Creatinine[i],NA)
Dvalue[i]<-ifelse(Dmarker[i]>0,Creatinine[i],NA)

#apply the optimal reference creatinine (which will later be chosen from yearmedian, weekmin or Dmin) across all of a patient's samples
yearmedian[i]<-ave(as.numeric(yearvalue),PatientID==PatientID,FUN = function(x) median(x, na.rm = TRUE))
weekmin[i]<-ave(as.numeric(weekvalue),PatientID==PatientID,FUN = function(x) min(x, na.rm = TRUE))
Dmin[i]<-ave(as.numeric(weekvalue),PatientID==PatientID,FUN = function(x) min(x, na.rm = TRUE))

#finds the index (investigated) creatinine of key interest, applies this value across all patient samples
mightyvalue<-NA
mightyvalue[i]<-Creatinine[i]
protagonistvalue<-ave(mightyvalue,PatientID,FUN = function(x) max(x, na.rm = TRUE))

#Finds the 48 hour difference in 2 samples taken within a 48 hour window
D<-NA
D[i]<-ifelse(as.numeric(Dmarker[i])>0,abs(Creatinine[i]-Dvalue[i]),D[i])
RV1[i]<-mightyvalue[i]/weekmin[PatientID==PatientID]
RV2[i]<-mightyvalue[i]/yearmedian[PatientID==PatientID]
RV_ratio<-ifelse(is.na(RV1),RV2,RV1)
##########################################################################################################CHECKED

AKIyear[i]<-ifelse((mightyvalue[i]/yearmedian[i])>=1.5 & !is.na(yearmedian[i]) & (mightyvalue[i]>(3*RV_ratio[i])|RV_ratio[i]>=1.5),1,0)
AKIweek[i]<-ifelse((mightyvalue[i]/weekmin[i])>=1.5 & !is.na(weekmin[i]) & (mightyvalue[i]>(3*RV_ratio[i])|RV_ratio[i]>=1.5),1,0)
AKIyear[i]<-ifelse(is.na(AKIyear[i]) & !is.na(yearmedian[i]) & RV_ratio[i]<1.5 & D[i]>26,1,AKIyear[i])
AKIweek[i]<-ifelse(is.na(AKIweek[i]) & weekmin[i]>0 & RV_ratio[i]<1.5 & D[i]>26,1,AKIweek[i])}

#END OF LOOP
close(pb)
##########################################################################################################CHECKED

#Split into AKI episodes
#DIVIDE AKI INTO 90 DAY SESSIONS

#DIVIDE AKI INTO 90 DAY SESSIONS
#Assign an akicounter value of 1 to all AKI blood tests
a<-a[AKIweek==1|AKIyear==1,]
a<-a[order(PatientID,event.date),]
a$AKIcounter<-1
a$PatientID<-as.factor(a$PatientID)
attach(a)
library(data.table)
a<-as.data.table(a)
a[, lag.date := c(NA, event.date[-.N]), by = PatientID]
a$add<-ifelse(as.numeric(event.date)-a$lag.date>90,1,0)
attach(a)

#Add +1 to aki counter if the date of the sample was >90 days after the first sample with an akicounter value of 1
a$add<-ifelse(is.na(a$add),0,a$add)
a$AKIcounter==1

for(x in 2:length(a$AKIcounter)){

a$AKIcounter[x]<-ifelse(PatientID[x]==PatientID[x-1],
a$AKIcounter[x-1]+a$add[x],1)}
summary(a$AKIcounter)
names(a)
a<-a[,c("PatientID","event.date","AKIcounter")]
save(a,file="AKIrecords.rda")
############################################################################################CHECKED
###############################################ATTACH AKI COUNTER TO CREA.REP TABLE
load("crea.rep.rda")
load("AKIrecords.rda")
indx1 <-neardate(crea.rep$PatientID, a$PatientID, crea.rep$event.date, a$event.date, 
                   best="prior")
crea.rep$NoAKIepisodes<-a[indx1, "AKIcounter"]
crea.rep$NoAKIepisodes<-unlist(crea.rep$NoAKIepisodes)
crea.rep$NoAKIepisodes<-ifelse(is.na(crea.rep$NoAKIepisodes),0,crea.rep$NoAKIepisodes)
save(crea.rep,file="crea.rep.rda")