bekmr:

### This is the wrapper function that fits a logistic regression model to the binomial outcome
### when the multivariate exposure is handled via KMR

### Notation:
###  y: nx1 vector of 0-1 for the n schools
###  X: nxp design matrix of covariates with first column being the intercept
###  beta: px1 vector of regression coefficients
###  z: nx(n-1) matrix with each row the vector of co-clustering probabilities for a school
###  h.z: nx1 vector with the values of the Gaussian process at "coordinates" z_1, z_2,..z_n
###  sigma.2: marginal variance of the Gaussian process prior on h
###  phi: range parameter of the Gaussian covariance function used to model the covariance
###       of the Gaussian process
###  m.beta: px1 prior mean of the multivariate normal distribution for the vector of
###          regression coefficients beta
###  V: pxp prior covariance matrix of the multivariate normal distribution for the vector of
###          regression coefficients beta
###  iV: inverse of the prior covariance matrix V for the beta vector



### As hyperparameters, I'd use:
###
### - m.beta the px1 vector of all 0's
### - V a diagonal pxp matrix with very large entries in the diagonal
###
### - shape.sigma.2  I would use a number very close to 2, so to have a large prior
###   variance
### - scale.sigma.2  I would use an empirical Bayes approach and maybe compute the average of the variance of all the z's
###
### - lower.phi and upper.phi   I would compute the distance matrix between the z's (see the Dist.mat definition in the other file), and compute the minimum and the maximum of all the non-diagonal elements of Dist.mat. I'd put lower.phi equal to the minimum and upper.phi equal to the maximum of these distances


### Initial values for the MCMC:
### - beta: I would run a logistic regression without the z's, and set the initial values for the beta equal to their estimates under a logistic regression
### - h.z   I would set all of them equal to 0
### - sigma.2   I would set it equal to the average of the variance of all the z's
### - phi    I would set it equal to a number drawn at random from a Uniform distribution with support the interval (lower.phi,upper.phi)


#
#init.beta <-
#init.h.z <-
#init.sigma.2 <-
#init.phi <-
#
#N.iter <-
#
### standard deviation of the proposal distribution used for the h function and teh
### beta function
#delta.h <-
#delta.beta <-
#
#### counter that keeps track of the number of acceptances for the h's and the beta
#accept.h <- 0
#accept.beta <- 0
#
### Number of MCMC iterations
#N.iter <-
#
#source("sampling_fcts_KMR.R")
#
### Number of MCMC iterations used for burnin
#burnin <-
#
#
#for(k in 1:N.iter){
#    if(k==1){
#        init.iR.phi <- comp.gau.icorrfct(Dist.mat.sqrd,init.phi))
#        new.h.z <- sample.h(y,init.h.z,X,init.beta,n,delta.h,init.iR.phi)
#        new.beta <- function(y,new.h.z,X,init.beta,p,delta.beta,m.beta,iV)
#        new.sigma.2 <- sample.sigma.2(new.h.z,n,init.iR.phi,shape.sigma2,scale.sigma2)
#        new.phi <- sample.phi(new.h.z,init.phi,init.iR.phi,Dist.mat.sqrd,lower.phi,upper.phi,n)
#    }
#    
#    if(k!=1){
#        iR.phi <- comp.gau.icorrfct(Dist.mat.sqrd,phi)
#        new.h.z <- sample.h(y,h.z,X,beta,n,delta.h,iR.phi)
#        new.beta <- function(y,new.h.z,X,beta,p,delta.beta,m.beta,iV)
#        new.sigma.2 <- sample.sigma.2(new.h.z,n,iR.phi,shape.sigma2,scale.sigma2)
#        new.phi <- sample.phi(new.h.z,phi,iR.phi,Dist.mat.sqrd,lower.phi,upper.phi,n)
#    }
#    
#    ifelse(sum(new.h.z==h.z)!=n,accept.h <- accept.h+1,accept.h <- accept.h)
#    ifelse(sum(new.beta==beta)!=p,accept.beta <- accept.beta+1,accept.beta <- accept.beta)
#    
#    h.z <- new.h.z
#    beta <- new.beta
#    sigma.2 <- new.sigma.2
#    phi <- new.phi
#
#    ### the treshold for acceptance can be lowered to match what suggested
#    ### in various papers; also the factors by which the sd of the proposals
#    ### are rescaled can be changed
#    if((k <= burnin) & ((k%%100)==0)){
#        if((accept.h/100) >= 0.4){
#            delta.h <- delta.h*1.5
#        }
#        if((accept.h/100) < 0.1){
#            delta.h <- delta.h/1.5
#        }
#        if((accept.beta/100) >= 0.4){
#            delta.beta <- delta.beta*1.5
#        }
#        if((accept.beta/100) < 0.1){
#            delta.beta <- delta.beta/1.5
#        }
#        accept.h <- 0
#        accept.beta <- 0
#    }
#
#    print(paste("Iteration ",k,sep=""))
#}
#
#