Install Packages & import library

#install.packages("geostatsp")
#install.packages("gambiaUTM")
#if (!requireNamespace("geostatsp", quietly = TRUE)) {
#    install.packages("geostatsp")
#}
#library(geostatsp)
data(package = "geostatsp")
data(gambiaUTM, package = "geostatsp")
str(gambiaUTM)
## Loading required package: terra
## Warning: package 'terra' was built under R version 4.3.3
## terra 1.7.71
## Formal class 'PackedSpatVector' [package "terra"] with 5 slots
##   ..@ type       : chr "points"
##   ..@ crs        : chr "BOUNDCRS[\n    SOURCECRS[\n        PROJCRS[\"unknown\",\n            BASEGEOGCRS[\"unknown\",\n                "| __truncated__
##   ..@ coordinates: num [1:2035, 1:2] 349631 349631 349631 349631 349631 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "x" "y"
##   ..@ index      : num [1:2035, 1:4] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:4] "geom" "part" "hole" "start"
##   ..@ attributes :'data.frame':  2035 obs. of  6 variables:
##   .. ..$ pos    : int [1:2035] 1 0 0 1 0 1 1 0 0 1 ...
##   .. ..$ age    : int [1:2035] 1783 404 452 566 598 590 589 609 791 829 ...
##   .. ..$ netuse : int [1:2035] 0 1 1 1 1 1 1 1 1 0 ...
##   .. ..$ treated: int [1:2035] 0 0 0 0 0 0 0 0 0 0 ...
##   .. ..$ green  : num [1:2035] 40.9 40.9 40.9 40.9 40.9 ...
##   .. ..$ phc    : int [1:2035] 1 1 1 1 1 1 1 1 1 1 ...
library(geostatsp)
## Warning: package 'geostatsp' was built under R version 4.3.3
## Loading required package: Matrix
Y <- gambiaUTM@attributes[[1]]
X <- scale(cbind(gambiaUTM@attributes[[2]],gambiaUTM@attributes[[3]],gambiaUTM@attributes[[4]],
                 gambiaUTM@attributes[[5]],gambiaUTM@attributes[[6]]))
n <- length(Y)

library(rjags)
## Warning: package 'rjags' was built under R version 4.3.3
## Loading required package: coda
## Warning: package 'coda' was built under R version 4.3.3
## 
## Attaching package: 'coda'
## The following objects are masked from 'package:terra':
## 
##     varnames, varnames<-
## Linked to JAGS 4.3.1
## Loaded modules: basemod,bugs
burn <- 1000
iters <- 5000
chains <- 1
mod <- textConnection("model{
  for(i in 1:n){
  Y[i] ~ dbern(pi[i])
  logit(pi[i]) <- beta[1] + X[i,1]*beta[2] + X[i,2]*beta[3] +
                  X[i,3]*beta[4] + X[i,4]*beta[5] + X[i,5]*beta[6]
  }
  for(j in 1:6){beta[j] ~ dnorm(0,0.01)}
}")
data <- list(Y=Y,X=X,n=n)
model <- jags.model(mod,data=data, n.chains = chains,quiet=TRUE)
update(model,burn, progress.bar="none")
samps <- coda.samples(model,variable.names=c("beta"),n.iter = iters,progress.bar="none")
summary(samps)
## 
## Iterations = 2001:7000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##             Mean      SD  Naive SE Time-series SE
## beta[1] -0.63195 0.04825 0.0006823      0.0008405
## beta[2]  0.27352 0.04742 0.0006706      0.0008801
## beta[3] -0.25345 0.05237 0.0007406      0.0012173
## beta[4] -0.12911 0.05982 0.0008460      0.0015516
## beta[5]  0.29539 0.05002 0.0007075      0.0008889
## beta[6] -0.09908 0.05276 0.0007461      0.0011747
## 
## 2. Quantiles for each variable:
## 
##            2.5%     25%      50%      75%     97.5%
## beta[1] -0.7273 -0.6641 -0.63267 -0.59919 -0.536197
## beta[2]  0.1826  0.2415  0.27334  0.30572  0.366771
## beta[3] -0.3543 -0.2885 -0.25239 -0.21905 -0.149144
## beta[4] -0.2441 -0.1693 -0.13000 -0.08906 -0.010283
## beta[5]  0.1981  0.2613  0.29610  0.32972  0.390024
## beta[6] -0.1988 -0.1353 -0.09988 -0.06379  0.005845
plot(samps)

set.seed(0820)
Y <- c(64,72,55,27,75,24,28,66,40,13)
N <- c(75,95,63,39,83,26,41,82,54,16)
q <- c(0.845,0.847,0.880,0.674,0.909,0.899,0.770,0.801,0.802,0.875)

X <- log(q)-log(1-q) #X=logit(q)
## Fit first model
library(rjags)
data <- list(Y=Y,N=N,X=X)
params <- c("beta")

model_string <- textConnection("model{
  # Likelihood
  for(i in 1:10){
  Y[i] ~ dbinom(p[i],N[i])
  logit(p[i]) <- beta[1] + beta[2]*X[i]
  }
  
  
  # Priors
  beta[1] ~ dnorm(0,0.01)
  beta[2] ~ dnorm(0,0.01)
}")
model <- jags.model(model_string,data=data,n.chains = 2,quiet = TRUE)
update(model,10000,progress.bar = "none")
samples1 <- coda.samples(model,variable.names = params,thin = 5, n.iter = 20000)
summary(samples1)
## 
## Iterations = 11005:31000
## Thinning interval = 5 
## Number of chains = 2 
## Sample size per chain = 4000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##            Mean     SD Naive SE Time-series SE
## beta[1] -0.1312 0.4164 0.004655       0.013896
## beta[2]  0.9711 0.2557 0.002858       0.008603
## 
## 2. Quantiles for each variable:
## 
##            2.5%     25%     50%    75%  97.5%
## beta[1] -0.9392 -0.4081 -0.1366 0.1452 0.6954
## beta[2]  0.4728  0.7992  0.9739 1.1431 1.4763
plot(samples1)

b1 <- c(samples1[[1]][,1],samples1[[2]][,1])
b2 <- c(samples1[[1]][,2],samples1[[2]][,2])
### Logistic Regression for Gambia data - Bernoulli likelihood ###
#library(geoR)
#Y <- gambia[,3]
#X <- scale(gambia[,4:8])
#n <- length(Y)