#######################################################################
# R function: This function fits a logisitic model to individual 
#             data according to the Bayesian methods discussed in McMahan et al.
# Input:
# Z = A matrix of testing responses whose jth row is of the form (col1=Z_j, col2=cj, col3=assay used, col4:col(4+cj-1)=indices of the individuals assigned to the jth pool) 
# X = Covariate matrix in which the ith row contains the covariate information pertaining to the ith individual
# Y = matrix whose ith row is of the form (col1=0, col2=number of pools ith individual involved in, col3:col(3+col2-1)=pools ith individual was assigned to)
# GI = The number of Gibbs iterates 
# b0 = Initialization value for the regression parameters
# na = Number of different assays being used
# ase(bse) = prior parameters for the testing sensitivities, if these values are to be estimated, note these should be of length na (Assuming a beta prior)
# asp(bsp) = prior parameters for the testing specificities, if these values are to be estimated, note these should be of length na (Assuming a beta prior)
# Se = vector of sensitivity values, if known
# Sp = vector of specificity values, if known
# est.error = logical if FALSE function uses the known Se and Sp to fit model, if TRUE Se and Sp are estimated along with the other parameters
#
#
# Required Packages: mnormt, matrixcalc
#
#
#

Bayes.reg.ind<-function(Y,X,GI,b0, ase=1, bse=1, asp=1, bsp=1, Se=NULL, Sp=NULL, est.error=FALSE, visual=FALSE, step=50){


N<-dim(X)[1]

beta<-matrix(-99,nrow=GI,ncol=length(b0))

if(est.error==TRUE){
Se.mat<-rep(-99,GI)
Sp.mat<-rep(-99,GI)
}


###############################
# Initializing latent variables
Y.t<-rbinom(N,1,exp(X%*%b0)/(1+exp(X%*%b0)))


##########################
#GIBBS SAMPLER STARTS HERE

for(g in 1:GI){


######################################################################################################## 
#Sampling the beta parameters (This assumes a diffuse prior easily changeable see Gamerman (1997)
beta[g,]<-WLS.Metropolis.Hastings.u(b0=b0,X=X,y=Y.t)
b0<-beta[g,]


if(est.error==TRUE){
###################
#Sampling Se and Sp

Se.1<-sum(Y*Y.t)
Se.2<-sum((1-Y)*Y.t)
Sp.1<-sum((1-Y)*(1-Y.t))
Sp.2<-sum(Y*(1-Y.t))

Se.mat[g]<-rbeta(1,ase+Se.1,bse+Se.2)
Sp.mat[g]<-rbeta(1,asp+Sp.1,bsp+Sp.2)

Se<-Se.mat[g]
Sp<-Sp.mat[g]
}



######################
# Sample true status Y

p<-as.vector(exp(X%*%b0)/(1+exp(X%*%b0)))
p1<- p*Se^(Y)*(1-Se)^(1-Y)
p2<- (1-p)*Sp^(1-Y)*(1-Sp)^(Y)
p.star<-p1/(p1+p2)
Y.t<-rbinom(N,1,p.star)




#print(g)

if(visual==TRUE){
if(sum(seq(100,GI,step)==g)==1){
par(mfrow=c(length(b0),1))
for(h in 1:length(b0)){
plot(1:g,beta[1:g,h],type="l")
lines(1:g,cumsum(beta[1:g,h])/(1:g),col="red")
}}
}


} # End Gibbs sampling 


if(est.error==TRUE){
return(list("beta"=beta, "Se"=Se.mat, "Sp"=Sp.mat))
}
if(est.error==FALSE){
return(list("beta"=beta))
}

} # End of function



##########################################################################
# Note to run code need package mnormt
library(mnormt)
library(matrixcalc)


#############################################################################################
# Carries out Metropolis hasting sampling as discussed in Gamerman (1997)
WLS.Metropolis.Hastings.u<-function(b0,X,y){

res<-mt_ct.u(b=b0, y=y, X=X, sig=0)
Mb0<-res$M
Cb0<-res$C
Lb0<-res$L
sig<-res$s


b.star<-as.vector(rmnorm(1,Mb0,Cb0))
qb.star<-dmnorm(b.star,as.vector(Mb0),Cb0)

res<-mt_ct.u(b=b.star, y=y, X=X, sig=sig)
Mb.star<-res$M
Cb.star<-res$C
Lb.star<-res$L

qb0<-dmnorm(b0,as.vector(Mb.star),Cb.star)

ratio.pi<-exp(Lb.star-Lb0)

alpha<-min(1,max(0,((ratio.pi*qb0)/(qb.star))))
accept<-rbinom(1,1,alpha)
if(accept==1){b1<-b.star}
if(accept==0){b1<-b0}

return("b1"=b1)
}

############################################
# Finds the Mean and Covariance of proposal 
# density according to Gammerman (1997)
mt_ct.u<-function(b,y,X,sig){
P<-length(b)
xb<-X%*%b
p<-exp(xb)/(1+exp(xb))

W<-p*(1-p)
Yb<-xb*(p*(1-p))+(y-p)
XWX<-t(X*matrix(W,ncol=P, nrow=length(W),byrow=FALSE))%*%X
if(is.non.singular.matrix(XWX)==FALSE || sig==1){
XWX<-XWX+0.1*diag(rep(1,P))
sig<-1
}

C<-solve(XWX)
C<-(C+t(C))/2
M<-C%*%(t(X)%*%(Yb))
L<- sum(y*(xb)-log(1+exp(xb)))  

return(list("M"=M, "C"=C, "L"=L, "s"=sig))
}