require(geoR)
require(bmstdr)
#upload entomological data (ento) and malaria data. Attach malaria data to ento table.
#Preparation for modelling

ento$LogGf1=log(ento$AnGambiae+1)
ento$LogFf1=log(ento$AnFunestus+1)
ento$LogNf1=log(ento$AnNili+1)
ento$LogMR=log(ento$Malariancidence)

#priors for range and sigma
gang=as.geodata(ento$AnGambiae, coords.col = 2:3, data.col = 5)
gangexv=variog(gang)
gangthv=eyefit(gangexv)

ganf=as.geodata(ento$AnFunestus, coords.col = 2:3, data.col = 7)
ganfexv=variog(ganf)
ganfthv=eyefit(ganfexv)

gann=as.geodata(ento$AnNili, coords.col = 2:3, data.col = 9)
gannexv=variog(gann)
gannthv=eyefit(gannexv)

gmal=as.geodata(ento$Malariancidence, coords.col = 2:3, data.col = 32)
gmalexv=variog(gmal)
gmalthv=eyefit(gmalexv)

ranges=list(gangthv[[1]]$practicalRange,ganfthv[[1]]$practicalRange,gannthv[[1]]$practicalRange,gmalthv[[1]]$practicalRange)
sigmas=list(gangthv[[1]]$cov.pars[1],ganfthv[[1]]$cov.pars[1],gannthv[[1]]$cov.pars[1],gmalthv[[1]]$cov.pars[1])

formula4=vector("list",4)
formula4[[1]]=as.formula(paste("LogGf1~offset(Week)+",paste(colnames(ento)[c(42,33,52)],collapse="+")))
formula4[[2]]=as.formula(paste("LogFf1~offset(Week)+",paste(colnames(ento)[c(23,21)],collapse="+")))
formula4[[3]]=as.formula(paste("LogNf1~offset(Week)+",paste(colnames(ento)[c(25)],collapse="+")))
formula4[[4]]=as.formula(paste("LogMR~offset(Week)+",paste(colnames(ento)[c(41,38,52)],collapse="+")))


fjgs=function(fx,dx,cx,nn=1000,bn=200,NN,ranges,sigmas){
	#fx             formula in list format, each elemnt is one process
	#dx             data
	#cx             coordinates
	#nn             MCMC sample size x iteration
	#bn             initial iterations to discard x iteration
	#NN             total number of iteration
	#ranges 	list of ranges of the same length as fx
	#sigmas 	marginal standard deviation of the field.


adpb=function (formula, data, 
    validrows = NULL,coords = 4:5, prior.beta0 = 0, prior.M = 1e-04, prior.sigma2 = c(2, 
        2), prior.tau2 = c(2, 0.1), prior.phi.param = NULL, cov.model = "exponential", 
    n.report = 500, verbose = FALSE, mchoice = TRUE, N = 5000, 
    burn.in = 1000, rseed = 44, plotit = TRUE) 
{
    set.seed(rseed)
    if (length(coords) == 2)coords <- as.matrix(unique(data[, coords]))
    if (length(validrows) > 0) {
        fdat <- data[(-validrows), ]
        vdat <- data[validrows, ]
        r <- nrow(vdat)
        n <- nrow(fdat)
        u <- getXy(formula = formula, data = vdat)
        xpreds <- u$X
        vdaty <- u$y
        vcoords <- coords[validrows, ]
        fcoords <- coords[-validrows, ]
    }
    else {
        fdat <- data
        n <- nrow(fdat)
        fcoords <- coords
    }
    alldistmat <- dist_mat(coords, coordtype)
    max.d <- max(alldistmat)
    distmat <- dist_mat(fcoords, coordtype)
    u <- getXy(formula = formula, data = fdat)
    X <- u$X
    y <- u$y
    p <- ncol(X)
    xnames <- colnames(X)[-1]
    if (length(prior.phi.param) < 2)prior.phi.param <- c(0.25, 1) * 3/max.d
    priors <- list(beta.Norm = list(prior.beta0, diag(prior.M^(-1), 
        p)), phi.Unif = prior.phi.param, sigma.sq.IG = prior.sigma2, 
        tau.sq.IG = prior.tau2)
    phistart <- mean(prior.phi.param)
    starting <- list(phi = phistart, sigma.sq = 1, tau.sq = 1)
    tuning <- list(phi = 0.3, sigma.sq = 0.3, tau.sq = 0.3)
    newy <- y
    fdat$newy <- newy
    newformula <- update(formula, newy ~ .)
    modfit <- spBayes::spLM(formula = newformula, data = fdat, 
        coords = fcoords, starting = starting, tuning = tuning, 
        priors = priors, cov.model = cov.model, n.samples = N, 
        verbose = verbose, n.report = N/n.report)
    modrecover <- spBayes::spRecover(modfit, start = burn.in + 
        1, verbose = FALSE)
    samps <- data.frame(modrecover$p.beta.recover.samples, modrecover$p.theta.recover.samples)
    params <- get_parameter_estimates(samps)
    u <- getXy(formula = formula, data = data)
    allX <- u$X
    fits <- as.vector(allX %*% params[1:p, 1])
    allres <- list(params = params, fit = modfit, max.d = max.d, 
        fitteds = fits)
    if (mchoice) {
        sn <- n
        betas <- samps[, 1:p]
        tau_sq <- samps$tau.sq
        sigma_sq <- samps$sigma.sq
        itmax <- length(tau_sq)
        yrep <- matrix(NA, nrow = itmax, ncol = n)
        xbeta <- t(X %*% t(betas))
        loglik <- matrix(NA, nrow = itmax, ncol = n)
        log_full_like_vec <- numeric()
        for (it in 1:itmax) {
            sigma2 <- sigma_sq[it]
            tau2 <- tau_sq[it]
            phi <- samps$phi[it]
            Sigma <- diag(tau2, nrow = sn, ncol = sn) + sigma2 * 
                exp(-phi * distmat)
            Qmat <- solve(Sigma)
            meanvec <- xbeta[it, ]
            meanmult <- diag(1/diag(Qmat), nrow = sn, ncol = sn) %*% 
                Qmat
            condmean <- newy - meanmult %*% (newy - meanvec)
            condvar <- 1/diag(Qmat)
            logden <- mnormt::dmnorm(newy, mean = meanvec, varcov = Sigma, 
                log = TRUE)
            log_full_like_vec[it] <- logden
            loglik[it, ] <- dnorm(newy, mean = condmean, sd = sqrt(condvar), 
                log = TRUE)
            yrep[it, ] <- condmean + rnorm(n) * sqrt(condvar)
        }
        sigma2 <- mean(sigma_sq)
        tau2 <- mean(tau_sq)
        Sigma <- diag(tau2, nrow = sn, ncol = sn) + sigma2 * 
            exp(-phi * distmat)
        meanxbeta <- apply(xbeta, 2, mean)
        logden <- mnormt::dmnorm(newy, mean = meanxbeta, varcov = Sigma, 
            log = TRUE)
        log_full_like_at_thetahat <- logden
        yrepmeans <- as.vector(apply(yrep, 2, mean))
        yrepvars <- as.vector(apply(yrep, 2, var))
        gof <- sum((newy - yrepmeans)^2)
        penalty <- sum(yrepvars)
        pmcc_results <- list(gof = gof, penalty = penalty, pmcc = gof + 
            penalty)
        waic_results <- calculate_waic(loglik)
        dic_results <- calculate_dic(log_full_like_at_thetahat, 
            log_full_like_vec)
        sbBayesmod <- c(unlist(dic_results), unlist(waic_results), 
            unlist(pmcc_results))
        allres$mchoice <- sbBayesmod
        allres$logliks <- list(log_full_like_at_thetahat = log_full_like_at_thetahat, 
            log_full_like_vec = log_full_like_vec, loglik = loglik)
        allres$spDiag.mchoice <- spBayes::spDiag(modrecover)
        if (verbose) {
            print(allres$mchoice)
            print(allres$spDiag.mchoice)
        }
    }
    if (length(validrows) > 0) {
        modpred <- spBayes::spPredict(modfit, pred.covars = xpreds, 
            pred.coords = vcoords, start = burn.in + 1, verbose = FALSE)
        ypreds <- modpred$p.y.predictive.samples
        predsums <- get_validation_summaries(t(ypreds))
        yvalidrows <- data.frame(vdat, predsums)
        b <- calculate_validation_statistics(vdaty, ypreds)
        if (plotit) 
            obs_v_pred_plot(vdaty, predsums)
        allres$stats <- b$stats
        allres$yobs_preds <- yvalidrows
        allres$valpreds <- t(ypreds)
    }
    allres$prior.phi.param <- prior.phi.param
    
    u <- getXy(formula = formula, data = data)
    y <- u$y
    allres$residuals <- y - allres$fitteds
    allres$sn <- nrow(data)
    allres$tn <- 1
    allres$formula <- formula
    allres$scale.transform <- scale.transform
    allres$call <- match.call()
    allres(results$params) <- c("mean", "sd", "2.5%", "97.5%")
    class(allres) <- "bmstdr"
    allres
}


	pn=length(fx)
	x=xx=logl=fitbeta=vector("list",pn)
	y=fs=numeric(pn)
	r=sample(1:1000000,NN)
	accepted=accepted2=prior.sigma2=prior.tau2=prior.phi=prior.sigma22=prior.tau22=prior.phi22=r*0
	prior.phi.i=as.numeric(quantile(dist(cx[,1:2]),prob=c(0.025,0.975)))
	priorphi=prior.phi.i

	for(i in 1:pn){ fitbeta[[i]]=coefficients(glm(fx[[i]],data=dx))
                        fs[i]=length(fitbeta[[i]])}

	for(i in 1:NN){
        id=0

        if(i==1){
                for(ii in 1:pn){
                print(c(i,ii))
                #print(fitbeta)
                        if(ii==1){
                                x[[ii]]=adpB(fx[[ii]], data=dx,coordtype="lonlat", coords=as.matrix(cx),prior.beta0=fitbeta[[ii]],prior.M=0.1, 
                                        mchoice=T, prior.phi.param=prior.phi.i,N=nn,burn.in=bn,prior.range=ranges[[ii]],prior.sigma=sigmas[[ii]],
                                        package="spBayes",rseed=r[i],n.report=0,verbose=F)
                                y[ii]=nrow(x[[ii]][[1]])
                                prior.sigma2.i=c(x[[ii]][[1]][c(y[ii]-2),1],x[[ii]][[1]][c(y[ii]-2),2])
                                prior.tau2.i=c(x[[ii]][[1]][c(y[ii]-1),1],x[[ii]][[1]][c(y[ii]-1),2])
                                prior.phi.i=c(x[[ii]][[1]][y[ii],3],x[[ii]][[1]][y[ii],4])
                                fitbeta[[ii]]=x[[ii]][[1]][1:fs[ii],1]

                        }
                        else{
                                x[[ii]]=adpB(fx[[ii]], data=dx,coordtype="lonlat", coords=as.matrix(cx),prior.beta0=fitbeta[[ii]],prior.M=0.1, 
                                        prior.sigma2=prior.sigma2.i,prior.tau2=prior.tau2.i,prior.phi.param=prior.phi.i, 
                                        mchoice=T,N=nn,burn.in=bn,prior.range=ranges[[ii]],prior.sigma=sigmas[[ii]],
                                        package="spBayes",rseed=r[i],n.report=0,verbose=F)
                                y[ii]=nrow(x[[ii]][[1]])
                                prior.sigma2.i=c(x[[ii]][[1]][c(y[ii]-2),1],x[[ii]][[1]][c(y[ii]-2),2])
                                prior.tau2.i=c(x[[ii]][[1]][c(y[ii]-1),1],x[[ii]][[1]][c(y[ii]-1),2])
                                prior.phi.i=c(x[[ii]][[1]][y[ii],3],x[[ii]][[1]][y[ii],4])
                                fitbeta[[ii]]=x[[ii]][[1]][1:fs[ii],1]
                        }
                        logl[[ii]]=rep(NA,NN)
                        logl[[ii]][i]=x[[ii]][[6]][[1]]
                }
        }
        else{
                for(ii in 1:pn){
 
                        xx[[ii]]=try(adpB(fx[[ii]], data=dx,coordtype="lonlat", coords=as.matrix(cx),prior.beta0=fitbeta[[ii]],prior.M=0.1, 
                                prior.sigma2=prior.sigma2.i,prior.tau2=prior.tau2.i,prior.phi.param=prior.phi.i, 
                                mchoice=T,N=nn,burn.in=bn,prior.range=ranges[[ii]],prior.sigma=sigmas[[ii]],
                                package="spBayes",rseed=r[i],n.report=0,verbose=F),silent=T)

                        if(class(xx[[ii]])=="bmstdr"){
                                prior.sigma2.i=c(xx[[ii]][[1]][c(y[ii]-2),1],xx[[ii]][[1]][c(y[ii]-2),2])
                                prior.tau2.i=c(xx[[ii]][[1]][c(y[ii]-1),1],xx[[ii]][[1]][c(y[ii]-1),2])
                                prior.phi.i=c(xx[[ii]][[1]][y[ii],3],xx[[ii]][[1]][y[ii],4])

                                logl[[ii]][i]=xx[[ii]][[6]][[1]]
                                if(logl[[ii]][i]>logl[[ii]][c(i-1)] | is.na(logl[[ii]][c(i-1)]))id=id+1
                                if(!is.na(logl[[ii]][i]) & !is.na(logl[[ii]][c(i-1)])){
                                        if(logl[[ii]][i]>(logl[[ii]][c(i-1)]-(logl[[ii]][c(i-1)]/10)))id=id+1
                                       if(logl[[ii]][i]>max(logl[[ii]][1:c(i-1)],na.rm=T))id=id+1

                                if(prior.sigma2.i[2]==0 | prior.tau2.i[2]==0 | prior.phi.i[1]==prior.phi.i[2]){
                                        Rx=sample(accepted2[accepted2!=0],1)
                                        prior.sigma2.i=c(prior.sigma2[Rx],prior.sigma22[Rx])
                                        prior.tau2.i=c(prior.tau2[Rx],prior.tau22[Rx])
                                        prior.phi.i=c(prior.phi[Rx],prior.phi22[Rx])
                                }

                                }
                        }
                        else{
                                print(xx[[ii]])
                                id=NA   
                        }
                
                 if(ii==pn){
                   if(!is.na(id)){
                        if(id>2){
                                 x=xx
                                 accepted[i]=max(accepted[c(1:i)])+1
                                 accepted2[i]=i
                                 save(x,file=paste("Best_",i,".Rdata",sep=""))
                                 if(id>3)save(x,file=paste("BestAll_",i,".Rdata",sep=""))
                                 for(iii in 1:pn)fitbeta[[iii]]=xx[[iii]][[1]][1:fs[iii],1]

                        }
                        else{
                                idi=pn
                                prior.sigma2.i=c(x[[idi]][[1]][c(y[idi]-2),1],x[[idi]][[1]][c(y[idi]-2),2])
                                prior.tau2.i=c(x[[idi]][[1]][c(y[idi]-1),1],x[[idi]][[1]][c(y[idi]-1),2])
                                prior.phi.i=c(x[[idi]][[1]][y[idi],3],x[[idi]][[1]][y[idi],4])
                                accepted[i]=max(accepted[c(1:i)])
                        }
                  }
                  else{
                                idi=pn
                                prior.sigma2.i=c(x[[idi]][[1]][c(y[idi]-2),1],x[[idi]][[1]][c(y[idi]-2),2])
                                prior.tau2.i=c(x[[idi]][[1]][c(y[idi]-1),1],x[[idi]][[1]][c(y[idi]-1),2])
                                prior.phi.i=c(x[[idi]][[1]][y[idi],3],x[[idi]][[1]][y[idi],4])
                                accepted[i]=max(accepted[c(1:i)])
                 }

                   prior.sigma2[i]=prior.sigma2.i[1]
                   prior.tau2[i]=prior.tau2.i[1]
                   prior.phi[i]=prior.phi.i[1]
                   prior.sigma22[i]=prior.sigma2.i[2]
                   prior.tau22[i]=prior.tau2.i[2]
                   prior.phi22[i]=prior.phi.i[2]
                }
            }
        }
	
	}#i


results=list(x,prior.sigma2,prior.tau2,prior.phi,logl,fitbeta,accepted,r)
names(results)=c("Final.models","prior.sigma2chain","prior.tau2chain","prior.phiChain","loglChain","fittedbetas","acceptedChain","randomSeeds")

results

}

#Joint modelling
joint=fjgs(fx=formula, dx=ento, cx=ento[,2:3], nn = 1000, bn = 200, NN=100, ranges=ranges,sigmas=sigmas)