####  get min. effect. dose ######
getMED <- function(obj,clinRel=0,gamma=0.05) {
    placebo.resp=obj$probs[1]
    index1=obj$probs>(placebo.resp+clinRel)
    index1=match(1,index1)
    LB=obj$probs - qnorm(1-gamma/2)*obj$se.probs
    index2=LB>placebo.resp
    index2=match(1,index2)
    obj$dose[max(index1,index2)]
}

#### fits cand. models to observed dose-resp data  ######

fitmodel <- function(dose, resp, M, clinRel=0, gamma=0.05, step.size=min(dose[-1]/10)) {
    if (is.null(M$model)) {M$model=cbind(dose); attr(M$model,"power")=1}     #if no model specified, assume linpred=dose
    ### fit pow and expo models ###
    fit=try(glm(resp~., family=M$family, data=as.data.frame(M$model)), silent=TRUE)       # fit power model
    sign=NA
    if (length(fit)>1) {
            d=seq(min(dose),max(dose),by=step.size)
            if (!is.null(attr(M$model,"power"))) {
            newdata=do.call("pow",list(dose=d,p=attr(M$model,"power"), dmax=attr(M$model,"dmax"),
                                       off=attr(M$model,"off"), scale=attr(M$model,"scale"))) }
            if (!is.null(attr(M$model,"expo"))) {
            newdata=do.call("expo",list(dose=d,order=attr(M$model,"expo"),
                                scale=attr(M$model,"scale"), off=attr(M$model,"off"))) }
            if (!is.null(attr(M$model,"loglog"))) {
            newdata=do.call("loglog",list(dose=d,p=attr(M$model,"power"), dmax=attr(M$model,"dmax"),
                                       off=attr(M$model,"off"), scale=attr(M$model,"scale"))) }

                                
            fitprob=predict(fit,type="response",se.fit=TRUE,newdata=as.data.frame(newdata))
            maxi=which.max(abs(fitprob$fit-fitprob$fit[1]))  
            slope=(fitprob$fit[maxi]-fitprob$fit[1])>0      # positive dose effect?
            sign=ifelse(slope>0,1,-1)                       # for signed T
            fit$T <- sign*(fit$null.deviance-fit$deviance) - 2*(fit$df.null-fit$df.residual)
            fit$pred <- list(dose=d,probs=fitprob$fit,se.probs=fitprob$se.fit)
            #fit$newdata<-as.data.frame(newdata)
            fit$MED <- getMED(fit$pred,clinRel,gamma)
            fit$type <- "GLM"
            fit$sign <- sign

    }
    return(fit)
}


##### fits cand. models to permuted data  ###
permfit <- function(obj, resp) {
    T=NA
    if (length(obj)<2) return(T)
        fit=try(do.call("glm",list(formula=resp~., family=obj$family, data=obj$data)), silent=TRUE)       # fit power model
        if (length(fit)>1 && fit$converged==TRUE) {
            fitprob=predict(fit,type="response",se.fit=TRUE,newdata=obj$newdata)
            maxi=which.max(abs(fitprob$fit-fitprob$fit[1]))  
            slope=(fitprob$fit[maxi]-fitprob$fit[1])>0      # positive dose effect?
            sign=ifelse(slope>0,1,-1)                       # for signed T
            T=sign*(fit$null.deviance-fit$deviance) - 2*(fit$df.null-fit$df.residual)  #signed, penalized T statistic
        }
    return(T)
}


#########################################
###    Fit candidate models            ###
### and find permutation distribution  ###
##########################################

permT <- function(dose, resp, models, perms=0, clinRel=0, alpha=0.025, gamma=0.05, label=names(models), trace=NULL, seed = NULL, steps=min(dose[-1])/10) {

#######################################################################
#### Get data in long form (important for permutations)             ###
#### First column: dose, Second column: resp (binary or continuous) ###
#######################################################################
if (dim(as.matrix(resp))[2]>1)  {
        n=rowSums(resp)
        nc=c(0,cumsum(n))
        obs=matrix(NA,nc[length(dose)+1],2)
        for (i in 1:length(dose)) {
            temp1=rep(dose[i],n[i])
            temp2=c(rep(1,resp[i,1]),rep(0,resp[i,2]))
            obs[(nc[i]+1):nc[i+1],]=cbind(temp1,temp2)
            }
        } else obs=cbind(dose,resp) 

m=length(models)
if (is.null(label)) label=paste("M", 1:m, sep = "")
fits=vector("list",m); names(fits)=label
obsT=aic=MED=df=matrix(NA,m,1)     
conv=pos.effect=matrix(0,m,1)
if (!is.null(seed)) set.seed(seed)
for (i in 1:m) {
    fit=fitmodel(dose,resp,models[[i]], clinRel=clinRel, gamma=gamma, step.size=steps)
    if (length(fit)>1) {
        conv[i]=(fit$converged==TRUE)
        fits[[i]]=fit
        obsT[i]=fit$T #signed T statistic
        aic[i]=fit$aic
        MED[i]=fit$MED
        df[i]=fit$df.null-fit$df.residual
        pos.effect[i]=fit$sign
    }
}   

perm.T=rbind(matrix(NA,perms,m),t(obsT)) # holds maximum penalized deviance for each perm;observed in last row
conv=rbind(matrix(0,perms,m),t(conv))    # convergence of the candidate models
getmin<-function(r) {if (all(is.na(r))) return(NA) else return(min(r,na.rm=TRUE))}
getmax<-function(r) {if (all(is.na(r))) return(NA) else return(max(r,na.rm=TRUE))}
perm=1
while (perm<=perms) {
    if (all(pos.effect==FALSE)) break
    T=matrix(NA,m,1)
    ############ permutations ##############
    obs.perm=sample(obs[,2],replace=FALSE)
    if (dim(as.matrix(resp))[2]>1) resp.perm=t(table(obs.perm,obs[,1]))[,2:1] else resp.perm=obs.perm          # new number of successes
    for (i in 1:m) {
        T[i]=permfit(fits[[i]],resp.perm) #signed T statistic
        conv[perm,i]=is.finite(T[i])
    }   
    perm.T[perm,]=t(T)
    if (!is.null(trace) && any(perm==seq(0,perms,trace))) {
        perm.T.int=perm.T[c(1:perm,perms+1),]
        rawP=((perm+2)-apply(perm.T.int,2,rank,na.last="keep"))/(perm+1) #matrix of raw P-values across all permutations
        perm.P=rawP[perm+1,]
        minP=apply(rawP,1,min,na.rm=TRUE) # minimum P-value across models
        c=quantile(minP,alpha)      # critical value
        #cat("Permutation",perm,": perm.P =",round(perm.P,4),",crit. val. minP =",round(c,4),"\n")
        cat("Permutation",perm,": crit. val. for minP =",round(c,5),"\n")
        flush.console()
        #smart permutation:
        #if (min(perm.P,na.rm=TRUE)>(c+2*sqrt(alpha*(1-alpha)/perm))) {perm.T=perm.T.int;perms=perm;break} 
        #if (min(perm.P,na.rm=TRUE)<(c-2*sqrt(alpha*(1-alpha)/perm))) {perm.T=perm.T.int;perms=perm;break} 
    }
    perm=perm+1
} #end looping over permutations

asympt.P=0.5*(1-pchisq(obsT+2*df,df))+0.5*pchisq(-(obsT+2*df),df)
if (perms>0) {
    rawP=((perms+2)-apply(perm.T,2,rank,na.last="keep"))/(perms+1) #matrix of raw P-values across all permutations
    perm.P=rawP[perms+1,]
    minP=apply(rawP,1,min,na.rm=TRUE) # minimum P-value across models
    c=quantile(minP,alpha)      # critical value
    minP.old=minP
    adj.P=matrix(NA,m,1)
    getmin<-function(r) {if (all(is.na(r))) return(NA) else return(min(r,na.rm=TRUE))}
    ind=order(perm.P,na.last = NA)
    j=1
    for (i in ind) {        #step-down minP procedure
        adj.P[i]=mean(minP<=perm.P[i],na.rm=TRUE) 
        if (j<(m-1)) minP=apply(rawP[, -ind[1:j]],1,getmin) else minP=rawP[, -ind[1:j]]
        if (j>1 && adj.P[ind[j]]<adj.P[ind[j-1]]) adj.P[ind[j]]=adj.P[ind[j-1]]
        j=j+1
        }
    }
index<-is.finite(MED)
weight=exp(0.5*obsT)*index/sum(exp(0.5*obsT[index==1]))
wMED=ifelse(sum(index)==0,NA,sum(MED*weight,na.rm=TRUE))
if (perms==0) results=list(model=label, aic=aic, obsT=obsT, alpha=alpha, perms=perms, df=df, pos.effect=pos.effect,
             asympt.P=asympt.P, MED=MED, wMED=wMED,weight=weight, conv=t(conv), fits=fits, dose=dose, clinRel=clinRel,  gamma=gamma)
    else results=list(model=label, aic=aic, obsT=obsT, minP=minP.old, crit.value=c, alpha=alpha, perms=perms, df=df, pos.effect=pos.effect,
             asympt.P=asympt.P, perm.T=perm.T, perm.P=perm.P, adj.P=adj.P, MED=MED, wMED=wMED,
             weight=weight, conv=apply(conv,2,mean), fits=fits, dose=dose, clinRel=clinRel,  gamma=gamma)
class(results)="bindr"
results
}

################################################
####### Summary functions print and plot  ######
######  for object "bindr"                ######
################################################

print.bindr = function(obj) {
    cat("Maximum T:", obj$model[which.max(obj$obsT)],"\n")
    cat("MED estimate for model with maximum T:", obj$MED[which.max(obj$obsT)],"\n")
    if (obj$perms>0) {
        cat("Minimum P-value:", obj$model[which.min(obj$perm.P)],"\n")
        cat("Proof of Concept P-value:",round(min(obj$adj.P,na.rm=TRUE),4),"\n")
        cat("MED estimate(s) for model(s) with minimum P-value:", obj$MED[obj$perm.P==min(obj$perm.P)],"\n")
    }
    cat("Weighted MED estimate:",obj$wMED ,"\n")
    }

summary.bindr <- function(obj) {
    oldClass(obj) <- "summary.bindr"
    obj
}

print.summary.bindr <- function(obj) {
    cat("Permutation analysis of minP:\n\n")
    cat("Results for individual models: \n")
    if (obj$perms==0)  print(data.frame(model=obj$model, aic=obj$aic, obsT=obj$obsT, asympt.P=obj$asympt.P, MED=obj$MED, weight=obj$weight, conv=obj$conv),digits=4)
    else print(data.frame(model=obj$model, aic=obj$aic, obsT=obj$obsT, asympt.P=obj$asympt.P, perm.P=obj$perm.P, adj.P=obj$adj.P, MED=obj$MED, weight=obj$weight, conv=obj$conv),digits=4)
    cat("\nMaximum T:",round(max(obj$obsT,na.rm=TRUE),2)," for model:", obj$model[which.max(obj$obsT)],"\n")
    cat("MED estimate for model with largest T:", obj$MED[which.max(obj$obsT)],"\n")
    if (obj$perms>0) {
        index=which(obj$perm.P == min(obj$perm.P,na.rm=TRUE))
        cat("Minimum P-value(s):",round(obj$perm.P[index],4)," for model(s):", obj$model[index],"\n")
        cat("Critical value:",obj$crit.value,"\n")
        cat("Proof of Concept P-value (",obj$perms,"Permutations ):",round(min(obj$adj.P,na.rm=TRUE),5),"\n")
        cat("MED estimate for model(s) with smallest P-value:", obj$MED[index],"\n")
        if (length(obj$perm.P==min(obj$perm.P,na.rm=TRUE))>1) cat("MED estimate for model with smallest P-value that has largest T:", obj$MED[order(1-obj$adj.P,obj$obsT,decreasing = TRUE)[1]],"\n")
    }
    cat("Weighted MED estimate:",obj$wMED ,"\n")
}

hist.bindr <-function(obj, xlab=quote(paste(min[s],p[s])), main="",...) {
    hist(obj$minP, main=main, xlab=xlab,...)
    points(obj$perm.P,y=rep(0,length(obj$model)),lwd=2,pch=19)
    abline(v=obj$crit.value,lty="dashed")
}




plot.bindr <- function(obj, which.model=order(1-obj$adj.P,obj$obsT,decreasing = TRUE)[1], se.fit=TRUE, plot.wMED=TRUE, legend=TRUE, plot.MED=FALSE, gamma=obj$gamma, ...) {
    if (is.character(which.model)) which.model=match(which.model,obj$model)
    plot(obj$fits[[which.model[1]]]$y~obj$dose,type="p", xlab="dose", ylab="efficacy", col=grey(0.2), pch=19,cex=1.1,...)
    #points((y/n)~dose,col=gray(0.2),pch=19,cex=1.1)
    for (i in 1:length(which.model)) {
        doses=obj$fits[[which.model[i]]]$pred$dose
        probs=obj$fits[[which.model[i]]]$pred$probs
        MED=obj$fits[[which.model[i]]]$MED
        model.label=obj$model[[which.model[i]]]
        lines(probs~doses, type="l", lty=i, lwd=2)
        #legend(x="topright",legend=model.label, bty="n")
        if (se.fit) {
            se.probs=obj$fits[[which.model[i]]]$pred$se.probs
            LB=probs-qnorm(1-gamma/2)*se.probs
            UB=probs+qnorm(1-gamma/2)*se.probs
            lines(LB~doses, type="l", lty=3)
            lines(UB~doses, type="l", lty=3)
        }
        if (length(which.model)==1 && obj$clinRel!=0) {
            abline(h=probs[1]+obj$clinRel, lty=2)
            legend(x=min(obj$dose),y=probs[1]+obj$clinRel, legend="clinical treshold",bty="n")
        }
        if (plot.MED) points(MED,min(obj$fits[[which.model[1]]]$y),pch="I",cex=1.3)
    }   
    if (plot.wMED) points(obj$wMED,min(obj$fits[[which.model[1]]]$y),pch="X",cex=1.3)
    if (legend) legend(x="bottomright",legend=obj$model[which.model],lty=1:length(which.model),lwd=2)
}