R/uniNetLeroux.R

In netcmc: Spatio-Network Generalised Linear Mixed Models for Areal Unit and Network Data

uniNetLeroux = function(formula, 
                        data, 
                        trials = NULL,
                        family = "gaussian",
                        squareSpatialNeighbourhoodMatrix,
                        spatialAssignment,
                        W, 
                        numberOfSamples = 10, 
                        burnin = 0, 
                        thin = 1,
                        seed = 1, 
                        trueBeta = NULL,
                        trueSpatialRandomEffects = NULL,
                        trueURandomEffects = NULL, 
                        trueSpatialTauSquared = NULL, 
                        trueSpatialRho = NULL, 
                        trueSigmaSquaredU = NULL,
                        trueSigmaSquaredE = NULL,
                        covarianceBetaPrior = 10^5,
                        a1 = 0.001, 
                        b1 = 0.001, 
                        a2 = 0.001, 
                        b2 = 0.001, 
                        a3 = 0.001, 
                        b3 = 0.001, 
                        centerSpatialRandomEffects = TRUE, 
                        centerURandomEffects = TRUE) {
  
  UniNetLerouxCall = match.call() 
  
  if(!family %in% c("gaussian", "binomial", "poisson")){
    stop("The family specified is not gaussian, binomial or poisson!")
  } 
  
  if(family == "gaussian"){
    
    output = univariateGaussianNetworkLerouxMH(formula = formula, 
                                             data = data, 
                                             squareSpatialNeighbourhoodMatrix = squareSpatialNeighbourhoodMatrix,
                                             spatialAssignment = spatialAssignment,
                                             W = W, 
                                             numberOfSamples = numberOfSamples, 
                                             burnin = burnin, 
                                             thin = thin,
                                             seed = seed, 
                                             trueBeta = trueBeta,
                                             trueSpatialRandomEffects = trueSpatialRandomEffects,
                                             trueURandomEffects = trueURandomEffects, 
                                             trueSpatialTauSquared = trueSpatialTauSquared, 
                                             trueSpatialRho = trueSpatialRho, 
                                             trueSigmaSquaredU = trueSigmaSquaredU,
                                             trueSigmaSquaredE = trueSigmaSquaredE,
                                             covarianceBetaPrior = covarianceBetaPrior,
                                             a1 = a1, 
                                             b1 = b1, 
                                             a2 = a2, 
                                             b2 = b2, 
                                             a3 = a3, 
                                             b3 = a3, 
                                             centerSpatialRandomEffects = centerSpatialRandomEffects, 
                                             centerURandomEffects = centerURandomEffects)
    
    results = list(call = UniNetLerouxCall,
                   y = output$y,
                   X = output$X,
                   standardizedX = output$standardizedX,
                   squareSpatialNeighbourhoodMatrix = output$squareSpatialNeighbourhoodMatrix ,
                   spatialAssignment = output$spatialAssignment ,
                   W = output$W,
                   samples = cbind(output$betaSamples, "tauSquared" = output$spatialTauSquaredSamples, 
                                   "rho" = output$spatialRhoSamples, "sigmaSquaredU" = output$sigmaSquaredUSamples,
                                   "sigmaSquaredE" = output$sigmaSquaredESamples),
                   betaSamples = output$betaSamples,
                   spatialTauSquaredSamples = output$spatialTauSquaredSamples,
                   spatialRhoSamples = output$spatialRhoSamples,
                   sigmaSquaredUSamples = output$sigmaSquaredUSamples,
                   sigmaSquaredESamples = output$sigmaSquaredESamples,
                   spatialRandomEffectsSamples = output$spatialRandomEffectsSamples,
                   uRandomEffectsSamples = output$uRandomEffectsSamples,
                   acceptanceRates = c(output$betaAcceptanceRate, output$spatialTauSquaredAcceptanceRate,
                                       output$spatialRhoAcceptanceRate, output$sigmaSquaredUAcceptanceRate,
                                       output$sigmaSquaredEAcceptanceRate),
                   spatialRandomEffectsAcceptanceRate = output$spatialRandomEffectsAcceptanceRate,
                   uRandomEffectsAcceptanceRate = output$uRandomEffectsAcceptanceRate,
                   timeTaken = output$timeTaken,
                   burnin = output$burnin,
                   thin = output$thin,
                   DBar = output$DBar,
                   posteriorDeviance = output$posteriorDeviance,
                   posteriorLogLikelihood = output$posteriorLogLikelihood,
                   pd = output$pd,
                   DIC = output$DIC)
    
  } else if(family == "poisson") {
    
    output = univariatePoissonNetworkLeroux(formula = formula, 
                                             data = data, 
                                             squareSpatialNeighbourhoodMatrix = squareSpatialNeighbourhoodMatrix,
                                             spatialAssignment = spatialAssignment,
                                             W = W, 
                                             numberOfSamples = numberOfSamples, 
                                             burnin = burnin, 
                                             thin = thin,
                                             seed = seed, 
                                             trueBeta = trueBeta,
                                             trueSpatialRandomEffects = trueSpatialRandomEffects,
                                             trueURandomEffects = trueURandomEffects, 
                                             trueSpatialTauSquared = trueSpatialTauSquared, 
                                             trueSpatialRho = trueSpatialRho, 
                                             trueSigmaSquaredU = trueSigmaSquaredU,
                                             covarianceBetaPrior = covarianceBetaPrior,
                                             a1 = a1, 
                                             b1 = b1, 
                                             a2 = a2, 
                                             b2 = b2, 
                                             centerSpatialRandomEffects = centerSpatialRandomEffects, 
                                             centerURandomEffects = centerURandomEffects)
    
    results = list(call = UniNetLerouxCall,
                   y = output$y,
                   X = output$X,
                   standardizedX = output$standardizedX,
                   squareSpatialNeighbourhoodMatrix = output$squareSpatialNeighbourhoodMatrix ,
                   spatialAssignment = output$spatialAssignment ,
                   W = output$W,
                   samples = cbind(output$betaSamples, "tauSquared" = output$spatialTauSquaredSamples, 
                                   "rho" = output$spatialRhoSamples, "sigmaSquaredU" = output$sigmaSquaredUSamples),
                   betaSamples = output$betaSamples,
                   spatialTauSquaredSamples = output$spatialTauSquaredSamples,
                   spatialRhoSamples = output$spatialRhoSamples,
                   sigmaSquaredUSamples = output$sigmaSquaredUSamples,
                   spatialRandomEffectsSamples = output$spatialRandomEffectsSamples,
                   uRandomEffectsSamples = output$uRandomEffectsSamples,
                   acceptanceRates = c(output$betaAcceptanceRate, output$spatialTauSquaredAcceptanceRate,
                                       output$spatialRhoAcceptanceRate, output$sigmaSquaredUAcceptanceRate),
                   spatialRandomEffectsAcceptanceRate = output$spatialRandomEffectsAcceptanceRate,
                   uRandomEffectsAcceptanceRate = output$uRandomEffectsAcceptanceRate,
                   timeTaken = output$timeTaken,
                   burnin = output$burnin,
                   thin = output$thin,
                   DBar = output$DBar,
                   posteriorDeviance = output$posteriorDeviance,
                   posteriorLogLikelihood = output$posteriorLogLikelihood,
                   pd = output$pd,
                   DIC = output$DIC)
    
  } else {
    
    output = univariateBinomialNetworkLeroux(formula = formula, 
                                               data = data, 
                                               trials = trials,
                                               squareSpatialNeighbourhoodMatrix = squareSpatialNeighbourhoodMatrix,
                                               spatialAssignment = spatialAssignment,
                                               W = W, 
                                               numberOfSamples = numberOfSamples, 
                                               burnin = burnin, 
                                               thin = thin,
                                               seed = seed, 
                                               trueBeta = trueBeta,
                                               trueSpatialRandomEffects = trueSpatialRandomEffects,
                                               trueURandomEffects = trueURandomEffects, 
                                               trueSpatialTauSquared = trueSpatialTauSquared, 
                                               trueSpatialRho = trueSpatialRho, 
                                               trueSigmaSquaredU = trueSigmaSquaredU,
                                               covarianceBetaPrior = covarianceBetaPrior,
                                               a1 = a1, 
                                               b1 = b1, 
                                               a2 = a2, 
                                               b2 = b2, 
                                               centerSpatialRandomEffects = centerSpatialRandomEffects, 
                                               centerURandomEffects = centerURandomEffects)
    
    results = list(call = UniNetLerouxCall,
                   y = output$y,
                   X = output$X,
                   standardizedX = output$standardizedX,
                   squareSpatialNeighbourhoodMatrix = output$squareSpatialNeighbourhoodMatrix ,
                   spatialAssignment = output$spatialAssignment ,
                   W = output$W,
                   samples = cbind(output$betaSamples, "tauSquared" = output$spatialTauSquaredSamples, 
                                   "rho" = output$spatialRhoSamples, "sigmaSquaredU" = output$sigmaSquaredUSamples),
                   betaSamples = output$betaSamples,
                   spatialTauSquaredSamples = output$spatialTauSquaredSamples,
                   spatialRhoSamples = output$spatialRhoSamples,
                   sigmaSquaredUSamples = output$sigmaSquaredUSamples,
                   spatialRandomEffectsSamples = output$spatialRandomEffectsSamples,
                   uRandomEffectsSamples = output$uRandomEffectsSamples,
                   acceptanceRates = c(output$betaAcceptanceRate, output$spatialTauSquaredAcceptanceRate,
                                       output$spatialRhoAcceptanceRate, output$sigmaSquaredUAcceptanceRate),
                   spatialRandomEffectsAcceptanceRate = output$spatialRandomEffectsAcceptanceRate,
                   uRandomEffectsAcceptanceRate = output$uRandomEffectsAcceptanceRate,
                   timeTaken = output$timeTaken,
                   burnin = output$burnin,
                   thin = output$thin,
                   DBar = output$DBar,
                   posteriorDeviance = output$posteriorDeviance,
                   posteriorLogLikelihood = output$posteriorLogLikelihood,
                   pd = output$pd,
                   DIC = output$DIC)
    
  }
  
  class(results) = "netcmc"
  
  return(results)
  
}