Repeatability analysis R script

A.repeatability.7.19.16.R

rm(list=ls())
library(ICC)
dat<-read.table('data20160211.csv',header=TRUE,sep=';')
library(outliers); library(lme4); library(lmerTest); library(boot)
library(multcomp)
types<-c("ADC1400", "SUB", "T2AX")
structures<-c("PeripheralZone", "TumorROI_PZ_1",  "NormalROI_PZ_1",
              "WholeGland")
get.dat<-function(quantity, structures, types, dat){
    X<-dat[dat$Quantity==quantity & !is.na(dat$Value) &
                (dat$Structure %in% structures) &
                (dat$ImageType %in% types),]
    X$ImageType<-as.factor(as.character(X$ImageType))
    X$Structure<-as.factor(as.character(X$Structure))
    X$PatientID<-as.factor(as.character(X$PatientID))
    return(X)}
get.par<-function(fit){
    mean<-as.numeric(fixef(fit)[1])
    sb<-as.numeric(attr(VarCorr(fit)$PatientID,'stddev'))
    se<-attr(VarCorr(fit),'sc')
    df<-summary(fit)$coefficients[,'df']
    b<-sqrt(df/(df-1/2))
    var.mean<-summary(fit)$vcov
    var.sd<-se^2 *(1-1/b^2)
    sd.rc.pct<-as.numeric(
                 sqrt(2)*qnorm(0.975) * 100*se/mean*
                 sqrt(var.sd/se^2+b^2*var.mean/mean^2))
    baseline<-seq(1,30,2); followup<-seq(2,30,2)
    Value<-slot(fit,"frame")$Value
    PatientID<-slot(fit,"frame")$PatientID
    icc<-sb^2/(sb^2+se^2)
    icc_ci<-ICCest(x=PatientID,y=Value)
    i<-!is.na(Value[baseline]) & !is.na(Value[followup])
    psDiff<-abs(Value[baseline][i]-Value[followup][i])
    psMean<-(Value[baseline][i]+Value[followup][i])/2
    meanDiff<-mean(psDiff)/mean(psMean)*100
    sdDiff<-sd(psDiff)/mean(psMean)*100
    rc<-qnorm(0.975)*sqrt(2)*se
    rc.pct<-rc/mean *100
    rc.pct.L<-rc.pct-qnorm(0.975)*sd.rc.pct
    rc.pct.U<-rc.pct+qnorm(0.975)*sd.rc.pct
    return(list(sb=sb, se=se,
                icc=icc, icc.L=icc_ci$LowerCI, icc.U=icc_ci$UpperCI,
                rc=rc, rc.pct=rc.pct, sd.rc.pct=sd.rc.pct,
                rc.pct.L=rc.pct.L, rc.pct.U=rc.pct.U,
                meanDiff=meanDiff, sdDiff=sdDiff))}
#
#    -- Volume for three structures (not normal), and three image types
dat.vol<-get.dat("Volume", structures[-3], types, dat)
#
#     -- Mean for ADC1400 for four structures
dat.mean<-get.dat("Mean", structures, "ADC1400", dat)
#
tbl<-NULL; rows<-NULL
#
#    -- One-way AOVA for volume, for each image.type, for each structure
quantity<-"Volume"; tp<-NULL
for(structure in structures[-3]){
    for(type in types){
      tmp<-get.dat("Volume", structure, type, dat.vol)
      n<-dim(tmp)[1]; n.sub<-length(levels(tmp$PatientID))
      fit<-lmer(Value ~(1|PatientID),data=tmp)
      par<-get.par(fit)
      tbl<-rbind(tbl,c(n,n.sub,par$sb,par$se,par$rc,par$rc.pct,par$sd.rc.pct,
                       par$rc.pct.L,par$rc.pct.U,par$icc,par$icc.L,
                       par$icc.U,par$meanDiff,par$sdDiff))
      rows<-c(rows,paste(quantity,structure,type,sep='.'))
    }}
quantity<-"Mean"
for(structure in structures){
  for(type in types[1]){
    tmp<-get.dat("Mean",structure,type,dat.mean)
    n<-dim(tmp)[1]; n.sub<-length(levels(tmp$PatientID))
    fit<-lmer(Value ~(1|PatientID),data=tmp)
    par<-get.par(fit)
    tbl<-rbind(tbl,c(n,n.sub,par$sb,par$se,par$rc,par$rc.pct,par$sd.rc.pct,
                 par$rc.pct.L,par$rc.pct.U,par$icc,par$icc.L,
                 par$icc.U,par$meanDiff,par$sdDiff))
    rows<-c(rows,paste(quantity,structure,type,sep='.'))
  }}
rownames(tbl)<-rows; tbl<-round(tbl,2)
colnames(tbl)<-c('N.Values','N.Patients','Between.Pt.SD','Within.Pt.SD','RC','RC.Pct','SD.RC.Pct','RC.Pct.L','RC.Pct.U',
                 'ICC','ICC lCI','ICC uCI','Mean % diff','SD % diff')
#
#    -- Method Comparison Contrasts
C<-matrix(c(0,1,0,
            0,0,1,
            0,1,-1),3,3,TRUE)
rownames(C)<-c("ADC-SUB=0",
               "ADC-T2AX=0",
               "SUB-T2AX=0")
#

tbl2<-NULL; rows<-NULL
#
#       -- Two-way AOVA for volume x image.type, for each structure
quantity<-"Volume"
for(structure in structures[-3]){
  tmp<-dat.vol[dat.vol$Structure==structure,]
  tmp$PatientID<-as.factor(as.character(tmp$PatientID))
  n<-dim(tmp)[1]; n.sub<-length(levels(tmp$PatientID))
  fit<-lmer(Value ~(1|PatientID)+ImageType,data=tmp)
  mean.mean<-fixef(fit)[1]
  sb<-as.numeric(attr(VarCorr(fit)$PatientID,'stddev'))
  se<-attr(VarCorr(fit),'sc')
  icc<-sb^2/(sb^2+se^2)
  rc<-qnorm(0.975)*sqrt(2)*se
  rc.pct<-rc/mean.mean *100
  tbl2<-rbind(tbl2,c(n,n.sub,sb,se,rc,rc.pct,icc))
  rows<-c(rows,paste(quantity,structure,sep='.'))
#
#  Compare methods, adjust for three comparisons
  colnames(C)<-names(fixef(fit))
  mc<-glht(fit,linfct=C)
  print(summary(mc))
}
rownames(tbl2)<-rows
colnames(tbl2)<-c('N.Values','N.Patients','Between.Pt.SD','Within.Pt.SD','RC','RC.Pct','ICC')
tbl2<-round(tbl2,2)