## read data
glaucoma <- read.csv("netmeta-example.csv")
## discard empty columns at right
glaucoma<-glaucoma[,1:9]

library(nlme)
library(rmeta)
to.trial<-function(dat,suffix=""){
        if (NROW(dat)==3) {
                return(rbind(to.trial(dat[c(1,2),],".a"),
                             to.trial(dat[c(1,3),],".b"),
                             to.trial(dat[c(2,3),],".c")))
        }
        rewt<-if (dat$arms[1]==3) 1/sqrt(2) else 1
        rval<-data.frame( Trial=paste(dat$Trial.Name[1],suffix,sep=""),
                    Trt1=dat$Treatment[1],Trt2=dat$Treatment[2],
                    PercChange1=dat$PercChange[1], PercChange2=dat$PercChange[2],
                    SE1=dat$SDPercChange[1]*rewt, SE2=dat$SDPercChange[2]*rewt
                   )
        return(rval)
}


mangle<-function(data,mu=1,sigma=2,treat=3,control=4,referent="none"){
    dd<-data.frame(mu=data[,mu],sigma=data[,sigma])
    trts<-sort(unique(c(as.character(data[,treat]),
            as.character(data[,control]))))
    for(i in seq(along=trts)){
        dd<-cbind(dd,(data[,treat]==trts[i])-(data[,control]==trts[i]))
    }
    names(dd)<-c("mu","sigma",trts)


    trts<-make.names(names(dd)[-(1:2)])
    #trts<-c("0",trts)
    trts<-c("0",trts[!(trts %in% referent)])
    ff<-paste("mu",paste(trts,collapse="+"),sep="~")
    ff<-as.formula(ff,env=parent.frame())
	
    tp<-function(x) paste(names(dd)[-(1:2)][ x[-(1:2)]!=0],collapse=":")
    trtpair<-apply(dd,1,tp)
    dd$trtpair<-factor(trtpair)


    names(dd)<-make.names(names(dd))
    row.names(dd)<-row.names(data)
    
     attr(dd,"formula")<-ff
    dd
}

nice.summary<-function(lmeobj){
  cat("\n",deparse(substitute(lmeobj)),"\n")
  
  x<-summary(lmeobj)
  
  print(summary(x$modelStruct), sigma = x$sigma, reEstimates = x$coef$random, 
        verbose = FALSE)
  
  cfs<-fixed.effects(lmeobj)
  vv<-lmeobj$varFix
  
  diff<- cfs
  SE<-sqrt(diag(vv))
  
  df<-data.frame(diff=round(diff,2),
                 low=round(diff-1.96*SE,2),
                 up=round(diff+1.96*SE,2))
  df$pvalue<-round(2*pnorm(-abs(diff/SE)),3)
  rownames(df)<-names(cfs)
  print(df)
  
}


glaucoma<-glaucoma[order(glaucoma$Trial.Number),]

## peak and trough are outcomes, not treatments

peak<-subset(glaucoma,regexpr("peak",as.character(Treatment))>0)
trough<-subset(glaucoma,regexpr("trough",as.character(Treatment))>0)
nrow(peak)+nrow(trough)==nrow(glaucoma)

peak$Treatment<-sub(" peak","",as.character(peak$Treatment))
trough$Treatment<-sub(" trough","",as.character(trough$Treatment))

## numbers of arms:
table(table(peak$Trial.Number))
table(table(trough$Trial.Number))

peak$arms<-table(peak$Trial.Number)[match(as.character(peak$Trial.Number),
                                          names(table(peak$Trial.Number)))]

trough$arms<-table(trough$Trial.Number)[match(as.character(trough$Trial.Number),
                                              names(table(trough$Trial.Number)))]


peaktrials<-do.call("rbind",by(peak,peak$Trial.Number,to.trial))
peaktrials$DiffPercChange<-with(peaktrials, PercChange1-PercChange2)
peaktrials$SEdiff<-with(peaktrials, sqrt(SE1^2+SE2^2))
peakdata<-mangle(peaktrials,mu="DiffPercChange",sigma="SEdiff",
                 treat="Trt1",control="Trt2",referent="timolol")

peak.lme<-lme(attr(peakdata,"formula"),~1|trtpair,data=peakdata,
              weights=varConstPower(form=~sigma, fixed=list(power=1)))
peak.lm<-lm(attr(peakdata,"formula"),data=peakdata,weight=1/sigma^2)

summary(peak.lm)
nice.summary(peak.lme)
metaplot(summary(peak.lme)$tTable[,1], summary(peak.lme)$tTable[,2],
          lab=names(fixed.effects(peak.lme)), 
          col=meta.colors(lines="thistle"),
          ylab="Intervention",xlab="Intraocular pressure difference (%)")
title(main="Peak: compared to timolol")
nr<-2+length(fixed.effects(peak.lme))
text(1,-nr, "Timolol better",adj=0)
text(-1,-nr, "Timolol worse",adj=1)
            

troughtrials<-do.call("rbind",by(trough,trough$Trial.Number,to.trial))
troughtrials$DiffPercChange<-with(troughtrials, PercChange1-PercChange2)
troughtrials$SEdiff<-with(troughtrials, sqrt(SE1^2+SE2^2))
troughdata<-mangle(troughtrials,mu="DiffPercChange",sigma="SEdiff",
                 treat="Trt1",control="Trt2",referent="timolol")

trough.lme<-lme(attr(troughdata,"formula"),~1|trtpair,data=troughdata,
              weights=varConstPower(form=~sigma, fixed=list(power=1)))
trough.lm<-lm(attr(troughdata,"formula"),data=troughdata,weight=1/sigma^2)

summary(trough.lm)
nice.summary(trough.lme)
metaplot(summary(trough.lme)$tTable[,1], summary(trough.lme)$tTable[,2],
          lab=names(fixed.effects(trough.lme)), 
          col=meta.colors(lines="thistle"),
          ylab="Intervention",xlab="Intraocular pressure difference (%)")
title(main="Trough: compared to timolol")
nr<-2+length(fixed.effects(trough.lme))
text(1,-nr, "Timolol better",adj=0)
text(-1,-nr, "Timolol worse",adj=1)