Conditional profile likelihood for balanced accuracy

proflike.R

llprofci<-function(tpos,fpos,tneg,fneg){

    ndis<-tpos+fneg
    nnot<-tneg+fpos

    lljoint<-function(a,d){
	k<-length(a)
	se<-a+d
	sp<-a-d
	rval <-dbinom(tpos,ndis,se,log=TRUE)+dbinom(tneg,nnot,sp,log=TRUE)	
	rval[is.nan(rval)]<- -Inf
	rval
    }
    
    llprof<-function(a){
        epsilon<-0.0001
	limits<-c(-a+epsilon, 1-a-epsilon)
        if (diff(limits)<=0) return(Inf)
	optimise(function(d) -2*lljoint(a,d),limits)$objective
    }
   
    xrange<-seq(0.5,.99,length=100)
    llprofcurve<-sapply(xrange, llprof)
    devmin<-min(llprofcurve,na.rm=TRUE)
    range(xrange[llprofcurve<=devmin+qchisq(.95,1)])
}


llcprofci<-function(tpos,fpos,tneg,fneg){

    ndis<-tpos+fneg
    nnot<-tneg+fpos

    lljoint<-function(a,d){
	k<-length(a)
	se<-a+d
	sp<-a-d
	rval <-dbinom(tpos,ndis,se,log=TRUE)+dbinom(tneg,nnot,sp,log=TRUE)	
	rval[is.nan(rval)]<- -Inf
	rval
    }

    lldinfo<-function(a,d){
	k<-length(a)
	se<-a+d
	sp<-a-d
        
        
        vardiff<-(se*(1-se)/ndis+sp*(1-sp)/nnot)/4
        info<-1/vardiff
        log(info)
    }

    
    llcprof<-function(a){
        epsilon<-0.0001
	limits<-c(-a+epsilon, 1-a-epsilon)
        if (diff(limits)<=0) return(Inf)
	llopt<-optimise(function(d) -2*lljoint(a,d),limits)
        jdd<-lldinfo(a,llopt$minimum)
        llopt$objective+jdd
    }
   
    xrange<-seq(0.5,.99,length=500)
    llcprofcurve<-sapply(xrange, llcprof)
    devmin<-min(llcprofcurve,na.rm=TRUE)
    range(xrange[llcprofcurve<=devmin+qchisq(.95,1)])
}


logit<-function(p)  log(p/(1-p))

expit<-function(eta) exp(eta)/(1+exp(eta))
                          

logitci<-function(tpos,fpos,tneg,fneg){

    ndis<-tpos+fneg
    nnot<-tneg+fpos

    sens<-tpos/ndis
    spec<-tneg/nnot

    bla<-(sens+spec)/2

    sdbla<-sqrt(sens*(1-sens)/ndis+spec*(1-spec)/nnot)/2

    sdlogit<-sdbla/(bla*(1-bla))

    logitbla<-logit(bla)

    expit(c(logitbla-1.96*sdlogit,logitbla+1.96*sdlogit))

}


normci<-function(tpos,fpos,tneg,fneg){

    ndis<-tpos+fneg
    nnot<-tneg+fpos

    sens<-tpos/ndis
    spec<-tneg/nnot

    bla<-(sens+spec)/2

    sdbla<-sqrt(sens*(1-sens)/ndis+spec*(1-spec)/nnot)/2


    c(bla-1.96*sdbla,bla+1.96*sdbla)

}