Matched case centile differences

setwd("c:/Users/mqbpjhr4/Documents")
library(lubridate)
alldata<-read.table("LMStableedited.csv",header=TRUE, sep=",")
length(alldata$ChildID)
alldata<-subset(alldata,!is.na(alldata$BMI),)
length(alldata$ChildID) #44 records removed
alldata$ChildID<-as.factor(alldata$ChildID)
alldata<-alldata[order(alldata$ChildID,alldata$AssessmentDate),]
alldata<-alldata[floor(alldata$ageattest/52.12)<12&floor(alldata$ageattest/52.12)>3,]
levels(alldata$ChildID) #63339 children
alldata$PostcodeDecile<-ceiling(alldata$PostcodeDecile)

alldata$DateOfBirth<-(as.Date(paste(alldata$DateOfBirth,"-01",sep="")))
alldata$AssessmentDate<-as.Date(paste(alldata$AssessmentDate,"-01",sep=""))
alldata$assyr<-year(as.POSIXlt(alldata$AssessmentDate, format="%Y-%m-%d"))
alldata$ageattestyr<-difftime(alldata$AssessmentDate,alldata$DateOfBirth,unit="weeks")/52.12
alldata$ageattestyr<-as.integer(floor(alldata$ageattestyr))
alldata<-alldata[!alldata$AcademicYear=="2011/2012"&!alldata$AcademicYear=="2012/2013",]
alldata$AcademicYear <- droplevels(alldata$AcademicYear)
attach(alldata)
alldata$ChildID<-as.factor(alldata$ChildID)
alldata$SchoolCode<-as.factor(alldata$SchoolCode)
alldata$IsAccountRegistered<-as.factor(alldata$IsAccountRegistered)
alldata<-alldata[order(alldata$ChildID,alldata$AssessmentDate),]
names(alldata)

#REMOVE POTENTIAL ERRONEOUS SAME DAY VALUES
d<-duplicated(alldata[,c(1,4)])
alldata2<-alldata[!d,]
length(alldata$ChildID)-length(alldata2$ChildID) #134 pairs of entries from the same date for the same child
# check duplicate entries from the same day- often there is an erroneous record with a correct record
#d4<-alldata2[alldata2$ChildID %in% alldata2$ChildID[d],]
#write.csv(d4,file="d4.csv")
#Manually review for ID reuse, save subset of reused IDs.
d4<-read.csv("d4.csv")

alldata2<-alldata2[!d,] #remove duplicate daily entries
length(alldata$ChildID)-length(alldata2$ChildID) #134 second entries from the same date removed where two daily entries are the same or differ within 1cm/0.5 kg
alldata2<-alldata2[!alldata2$ChildID %in% d4$ChildID,] #remove full records of children with inconsistent multiple daily records
length(alldata$ChildID)-length(alldata2$ChildID) 
#63 second entries removed from IDs potential reused child Ior associated with erroneous values

d2<-duplicated(alldata2[,c(1,3)])
alldata2<-alldata2[!d2,]
length(alldata$ChildID)-length(alldata2$ChildID) #308 second and later entries removed from the same assessment year
alldata<-alldata2
attach(alldata)

#DRAW IN BMI THRESHOLD TO WORK OUT CENTILES FOR OVERWEIGHT, OBESE, SEVERELY OBESE
alldata$weightcat<-1#normal
alldata$weightcat<-ifelse(alldata$Centile_BMI>=90.879,2,alldata$weightcat)#overweight
alldata$weightcat<-ifelse(alldata$Centile_BMI>=97.725,3,alldata$weightcat)#obese
alldata$weightcat<-ifelse(alldata$Centile_BMI<2.275,4,alldata$weightcat)#underweight
alldata$weightcat<-ifelse(alldata$Centile_BMI>=99.617,5,alldata$weightcat)#severe
alldata$weightcat<-ifelse(alldata$Centile_BMI<0.383,6,alldata$weightcat)#very underweight
table(alldata$weightcat)
alldata$NoRecords<-1
attach(alldata)

#Assign academic year
alldata$Record1<-as.numeric(ifelse(alldata$AcademicYear=="2013/2014",1,0))
alldata$Record2<-as.numeric(ifelse(alldata$AcademicYear=="2014/2015",1,0))
alldata$Record3<-as.numeric(ifelse(alldata$AcademicYear=="2015/2016",1,0))
alldata$Record4<-as.numeric(ifelse(alldata$AcademicYear=="2016/2017",1,0))

#Print BMI at each test
alldata$BMI1<-as.numeric(ifelse(alldata$AcademicYear=="2013/2014",BMI,NA))
alldata$BMI2<-as.numeric(ifelse(alldata$AcademicYear=="2014/2015",BMI,NA))
alldata$BMI3<-as.numeric(ifelse(alldata$AcademicYear=="2015/2016",BMI,NA))
alldata$BMI4<-as.numeric(ifelse(alldata$AcademicYear=="2016/2017",BMI,NA))

#Print weightcat at each age
alldata$weightcat1<-as.numeric(ifelse(alldata$AcademicYear=="2013/2014",weightcat,NA))
alldata$weightcat2<-as.numeric(ifelse(alldata$AcademicYear=="2014/2015",weightcat,NA))
alldata$weightcat3<-as.numeric(ifelse(alldata$AcademicYear=="2015/2016",weightcat,NA))
alldata$weightcat4<-as.numeric(ifelse(alldata$AcademicYear=="2016/2017",weightcat,NA))
attach(alldata)

#Print centile at each age
alldata$cent1<-as.numeric(ifelse(alldata$AcademicYear=="2013/2014",Centile_BMI,NA))
alldata$cent2<-as.numeric(ifelse(alldata$AcademicYear=="2014/2015",Centile_BMI,NA))
alldata$cent3<-as.numeric(ifelse(alldata$AcademicYear=="2015/2016",Centile_BMI,NA))
alldata$cent4<-as.numeric(ifelse(alldata$AcademicYear=="2016/2017",Centile_BMI,NA))
attach(alldata)

#Print date at each age
alldata$d1<-as.numeric(ifelse(alldata$AcademicYear=="2013/2014",AssessmentDate,NA))
alldata$d2<-as.numeric(ifelse(alldata$AcademicYear=="2014/2015",AssessmentDate,NA))
alldata$d3<-as.numeric(ifelse(alldata$AcademicYear=="2015/2016",AssessmentDate,NA))
alldata$d4<-as.numeric(ifelse(alldata$AcademicYear=="2016/2017",AssessmentDate,NA))
attach(alldata)

#Print z at each age
alldata$z1<-as.numeric(ifelse(alldata$AcademicYear=="2013/2014",SDS_BMI,NA))
alldata$z2<-as.numeric(ifelse(AcademicYear=="2014/2015",SDS_BMI,NA))
alldata$z3<-as.numeric(ifelse(alldata$AcademicYear=="2015/2016",SDS_BMI,NA))
alldata$z4<-as.numeric(ifelse(alldata$AcademicYear=="2016/2017",SDS_BMI,NA))

#Determine how many measurements taken per child

alldata$NoRecords<-(Record1+Record2+Record3+Record4)
Recs<-aggregate(alldata$NoRecords, by=list(Category=ChildID), FUN=sum)
names(Recs)[names(Recs)=="x"] <- "NoRecs"
names(Recs)[names(Recs)=="Category"] <- "ChildID"

##############################################################################

#Merge to have BMI at different ages on the same line
#Determine how many measurements taken per child

Recs2<-Recs[Recs$NoRecs>1,]
multiple<-alldata[alldata$ChildID %in% Recs2$ChildID,] #children measured more than once
length(unique(multiple$ChildID)) #34768 children
Ethnicity2<-multiple$Ethnicity
levels(Ethnicity2)[1:3]<-"1" #white
levels(Ethnicity2)[10:12]<-"2" #black
levels(Ethnicity2)[2:3]<-"3" #mixed black and white ethnicities
levels(Ethnicity2)[c(5:8,10)]<-"4" #asian
levels(Ethnicity2)[c(3:4,7)]<-"5" #other
levels(Ethnicity2)[c(6:7)]<-NA
levels(Ethnicity2)
multiple$Eth2<-as.numeric(as.character(Ethnicity2))

multiplea<-multiple[!is.na(multiple$Eth2),c(1, 44)]
attach(multiple)
multiple<-with(multiple, multiple[order(ChildID),])

multiple$ChildID<-as.factor(multiple$ChildID)
multiple$IsAccountRegistered<-as.numeric(multiple$IsAccountRegistered)-1
names(multiple)
multiplebase<-multiple[,c(1,12)]
multiple1<-multiple[!is.na(multiple$cent1),c(1, 32)]
multiple2<-multiple[!is.na(multiple$cent2),c(1, 33)]
multiple3<-multiple[!is.na(multiple$cent3),c(1, 34)]
multiple4<-multiple[!is.na(multiple$cent4),c(1, 35)]
multiple5<-multiple[c(1, 10)]
library(plyr)

multiple6<-multiple[!is.na(multiple$d1),c(1, 36)]
multiple7<-multiple[!is.na(multiple$d2),c(1, 37)]
multiple8<-multiple[!is.na(multiple$d3),c(1, 38)]
multiple9<-multiple[!is.na(multiple$d4),c(1, 39)]
multiple10<-multiple[!is.na(multiple$weightcat1),c(1, 28)]
multiple11<-multiple[!is.na(multiple$weightcat2),c(1, 29)]
multiple12<-multiple[!is.na(multiple$weightcat3),c(1, 30)]
multiple13<-multiple[!is.na(multiple$weightcat4),c(1, 31)]
multiple14<-multiple[!is.na(multiple$z1),c(1, 40)]
multiple15<-multiple[!is.na(multiple$z2),c(1, 41)]
multiple16<-multiple[!is.na(multiple$z3),c(1, 42)]
multiple17<-multiple[!is.na(multiple$z4),c(1, 43)]
attach(multiple)
multiple18<-aggregate(multiple$ageattestyr,by=list(multiple$ChildID),FUN=min)
colnames(multiple18)<- c("ChildID","AgeF")
multiple19<-aggregate(as.numeric(multiple$Gender),by=list(ChildID),FUN=min)
colnames(multiple19)<- c("ChildID","Gender")
attach(multiple)

multi<-merge(multiplebase,multiplea,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple1,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple2,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple3,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple4,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple5,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple6,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple7,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple8,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple9,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple10,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple11,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple12,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple13,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple14,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple15,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple16,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple17,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple18,by="ChildID",all.x=TRUE)
multi<-merge(multi,multiple19,by="ChildID",all.x=TRUE)
attach(multi)

multi2<-aggregate(multi,by=list(multi$ChildID),FUN=mean,na.rm=TRUE)
multi2$ChildID<-multi2$Group.1
multi<-multi2
save(multi,file="multi.rda")
###############################################

#CALCULATE WEIGHT CATEGORY TRANSITIONS BETWEEN FIRST AND LAST MEASUREMENTS

length(multi$ChildID)#34713
attach(multi)

multi$measure1<-as.numeric(multi$weightcat1)
multi$measure1<-ifelse(is.na(multi$measure1)&!is.na(multi$weightcat2),as.numeric(multi$weightcat2),multi$measure1)
multi$measure1<-ifelse(is.na(multi$measure1)&!is.na(multi$weightcat3),as.numeric(multi$weightcat3),multi$measure1)
multi$measure1<-ifelse(is.na(multi$measure1)&!is.na(multi$weightcat4),paste(multi$weightcat4),multi$measure1)
table(multi$measure1)

multi$measure2<-NA
multi$measure2<-ifelse(!is.na(multi$weightcat1)&!is.na(multi$weightcat2),as.numeric(multi$weightcat2),multi$measure2)
multi$measure2<-ifelse(is.na(multi$measure2)&(!is.na(multi$weightcat1)|!is.na(multi$weightcat2))&(!is.na(multi$weightcat3)),as.numeric(multi$weightcat3),multi$measure2)
multi$measure2<-ifelse(is.na(multi$measure2)&(!is.na(multi$weightcat1)|!is.na(multi$weightcat2)|!is.na(multi$weightcat3))&(!is.na(multi$weightcat4)),as.numeric(multi$weightcat4),multi$measure2)
table(multi$measure2)

table(multi$measure1,multi$measure2)

#CALCULATE CENTILE DIFFERENCES BETWEEN FIRST AND LAST MEASUREMENTS

multi$cen1<-as.numeric(multi$cent1)
multi$cen1<-ifelse(is.na(multi$cen1)&!is.na(multi$cent2),as.numeric(multi$cent2),multi$cen1)
multi$cen1<-ifelse(is.na(multi$cen1)&!is.na(multi$cent3),as.numeric(multi$cent3),multi$cen1)
multi$cen1<-ifelse(is.na(multi$cen1)&!is.na(multi$cent4),as.numeric(multi$cent4),multi$cen1)

multi$cen2<-NA
multi$cen2<-ifelse(!is.na(multi$cent1)&!is.na(multi$cent2),as.numeric(multi$cent2),as.numeric(multi$cen2))
multi$cen2<-ifelse(is.na(multi$cen2)&(!is.na(multi$cent1)|!is.na(multi$cent2))&!is.na(multi$cent3),as.numeric(multi$cent3),multi$cen2)
multi$cen2<-ifelse(is.na(multi$cen2)&(!is.na(multi$cent1)|!is.na(multi$cent2)|!is.na(multi$cent3))&(!is.na(multi$cent4)),as.numeric(multi$cent4),multi$cen2)

summary(multi$cen2)

#CALCULATE Z SCORE DIFFERENCES BETWEEN MEASUREMENTS

multi$zs1<-as.numeric(multi$z1)
multi$zs1<-ifelse(is.na(multi$z1)&!is.na(multi$z2),as.numeric(multi$z2),multi$zs1)
multi$zs1<-ifelse(is.na(multi$z1)&!is.na(multi$z3),as.numeric(multi$z3),multi$zs1)
multi$zs1<-ifelse(is.na(multi$z1)&!is.na(multi$z4),as.numeric(multi$z4),multi$zs1)

multi$zs2<-NA
multi$zs2<-ifelse(!is.na(multi$z1)&!is.na(multi$z2),as.numeric(multi$z2),as.numeric(multi$zs2))
multi$zs2<-ifelse(is.na(multi$zs2)&(!is.na(multi$z1)|!is.na(multi$z2))&!is.na(multi$z3),as.numeric(multi$z3),multi$zs2)
multi$zs2<-ifelse(is.na(multi$zs2)&(!is.na(multi$z1)|!is.na(multi$z2)|!is.na(multi$z3))&(!is.na(multi$z4)),as.numeric(multi$z4),multi$zs2)

multi$ds1<-as.numeric(multi$d1)
multi$ds1<-ifelse(is.na(multi$ds1)&!is.na(multi$d2),as.numeric(multi$d2),as.numeric(multi$ds1))
multi$ds1<-ifelse(is.na(multi$ds1)&!is.na(multi$d3),as.numeric(multi$d3),as.numeric(multi$ds1))
multi$ds1<-ifelse(is.na(multi$ds1)&!is.na(multi$d4),as.numeric(multi$d4),as.numeric(multi$ds1))

multi$ds2<-NA
multi$ds2<-ifelse(!is.na(multi$d1)&!is.na(multi$d2),as.numeric(multi$d2),as.numeric(multi$ds2))
multi$ds2<-ifelse(is.na(multi$ds2)&(!is.na(multi$d1)|!is.na(multi$d2))&!is.na(multi$d3),as.numeric(multi$d3),as.numeric(multi$ds2))
multi$ds2<-ifelse(is.na(multi$ds2)&(!is.na(multi$d1)|!is.na(multi$d2)|!is.na(multi$d3))&(!is.na(multi$d4)),as.numeric(multi$d4),as.numeric(multi$ds2))

multi$cendiff<-multi$cen2-multi$cen1
days<-as.numeric(multi$ds2)-as.numeric(multi$ds1)
multi$cendiffs<-multi$cendiff/days*365.25

multi$zdiff<-multi$zs2-multi$zs1
multi$zdiffs<-multi$zdiff/days*365.25

multi$diffnorm<-abs(as.numeric(multi$cen1)-50)-abs(as.numeric(multi$cen2)-50) #if positive, marks a move towards the normal
multi$diffnorms<-multi$diffnorm/days*365.25
multi$diffznorm<-(abs(0-multi$zs1)-abs(0-multi$zs2)) #if positive, marks a move towards the normal
attach(multi)

multi$dir<-ifelse(multi$cendiff>0,1,0)
multi$dirz<-ifelse(multi$zdiff>0,1,0)

table(multi$IsAccountRegistered)

multi<-multi[multi$cen1>=2.275,]
multi$bmidecile<-findInterval(multi$cen1,c(10,20,30,40,50,60,70,80,90))
multi$cat1<-findInterval(multi$cen1,c(10,20,30,40,50,60,70,80,90))
multi$cat2<-findInterval(multi$cen1,c(10,20,30,40,50,60,70,80,90))
multi$imdquintile<-findInterval(multi$PostcodeDecile,c(0,2,4,6,8,10))
multi$AgeF<-floor(multi$AgeF)
attach(multi)

multi<-multi[,c(2,4,9,22:44)]
attach(multi)

#BREAKPOINT
##################################
##############
#Find difference between 1st and last measurements


#CREATE SUBSETS BY BMI CATEGORY USING CLINICAL THRESHOLDS
multi$IsAccountRegistered<-as.factor((as.numeric(multi$IsAccountRegistered)))
multi<-multi[!is.na(multi$Eth2)&!is.na(multi$imdquintile),]

attach(multi)
length(multi$ChildID)

save(multi,file="multiformatching.rda")
summary(multi$ds2-multi$ds1)
over<-multi[multi$cen1>=90.879&multi$cen1<97.725,]
normal<-multi[multi$cen1>=2.275&multi$cen1<90.879,]
obese<-multi[multi$cen1>=97.725,]

#MATCH CASES AND CONTROLS
##############################################
library(MatchIt)#missing data cause failure even if they are in fields not included in the matching terms
set.seed(1234)

match.ov<- matchit(IsAccountRegistered~AgeF+Gender+Eth2+imdquintile,data=over, method="exact", ratio=2)
match.ob<- matchit(IsAccountRegistered~AgeF+Gender+Eth2+imdquintile,data=obese, method="exact", ratio=2)
match.no<- matchit(IsAccountRegistered~AgeF+Gender+Eth2+bmidecile+imdquintile,data=normal, method="exact", ratio=2)

matchd.no<- match.data(match.no,group="all", distance="distance",weights="weights", subclass="subclass")
matchd.ob<- match.data(match.ob,group="all", distance="distance",weights="weights", subclass="subclass")
matchd.ov<- match.data(match.ov,group="all", distance="distance",weights="weights", subclass="subclass")

normno<-matchd.no[matchd.no$IsAccountRegistered==0,]
normyes<-matchd.no[matchd.no$IsAccountRegistered==1,]

obyes<-matchd.ob[matchd.ob$IsAccountRegistered==1,]
obno<-matchd.ob[matchd.ob$IsAccountRegistered==0,]

ovyes<-matchd.ov[matchd.ov$IsAccountRegistered==1,]
ovno<-matchd.ov[matchd.ov$IsAccountRegistered==0,]

#CALCULATE CENTILE CHANGE IN MATCHED CHILDREN

mean(ovyes$cendiffs)
mean(ovno$cendiffs)
mean(obno$cendiffs)
mean(obyes$cendiffs)
mean(normyes$cendiffs)
mean(normno$cendiffs)

length(obno$ChildID[obno$diffnorms>0])
length(obno$ChildID[obno$diffnorms==0])
length(obno$ChildID[obno$diffnorms<0])

length(ovno$ChildID[ovno$diffnorms>0])
length(ovno$ChildID[ovno$diffnorms==0])
length(ovno$ChildID[ovno$diffnorms<0])

length(normno$ChildID[normno$diffnorms>0])
length(normno$ChildID[normno$diffnorms==0])
length(normno$ChildID[normno$diffnorms<0])

length(obyes$ChildID[obyes$diffnorms>0])
length(obyes$ChildID[obyes$diffnorms==0])
length(obyes$ChildID[obyes$diffnorms<0])

length(ovyes$ChildID[ovyes$diffnorms>0])
length(ovyes$ChildID[ovyes$diffnorms==0])
length(ovyes$ChildID[ovyes$diffnorms<0])

length(normyes$ChildID[normyes$diffnorms>0])
length(normyes$ChildID[normyes$diffnorms==0])
length(normyes$ChildID[normyes$diffnorms<0])

median(ovyes$cendiffs)
median(ovno$cendiffs)
median(obno$cendiffs)
median(obyes$cendiffs)
median(normyes$cendiffs)
median(normno$cendiffs)

mean(ovyes$diffnorms)
mean(ovno$diffnorms)
mean(obno$diffnorms)
mean(obyes$diffnorms)
mean(normyes$diffnorms)
mean(normno$diffnorms)

median(ovyes$diffnorms)
median(ovno$diffnorms)
median(obno$diffnorms)
median(obyes$diffnorms)
median(normyes$diffnorms)
median(normno$diffnorms)

median(ovyes$zdiffs)
median(ovno$zdiffs)
median(obno$zdiffs)
median(obyes$zdiffs)
median(normyes$zdiffs)
median(normno$zdiffs)