# 1. An example for fitting spatial generalized linear mixed models with random projections to binary observations.

### a. Simulate date using parameter values: $$\nu = 2.5, \sigma^2 = 1, \phi = 0.2, \beta = (1,1).$$

library(fastLaplace)
if(requireNamespace("mgcv")){
sigma2 = 1
phi = 0.2
beta.true = c(1,1)
n = 400
n.pred = 100
coords.all<- matrix(runif((n+n.pred)*2),ncol=2,nrow=n+n.pred) # simulate data locations
X.all <- matrix(runif((n+n.pred)*2),ncol=2,nrow=(n+n.pred))
suppressMessages(require(fields))
dist.all <- fields::rdist(coords.all,coords.all) # compute distance matrix
matern <- function(phi,mat.dist){
K = (1+sqrt(5)/phi*mat.dist+ 5/(3*phi^2)*mat.dist^2)*exp(-sqrt(5)/phi*mat.dist)
return(K)
} # matern 2.5
V.all <- sigma2*matern(phi,dist.all) # compute covariance matrix
set.seed(1)
r.e.all <- mgcv::rmvn(1,rep(0,nrow(coords.all)),V.all) # simulate random effects
pi.all <- X.all%*%beta.true + r.e.all # linear model
p.all <- exp(pi.all)/(1+exp(pi.all)) # compute the probability of z = 1 for binary process
Y.all <- sapply(p.all, function(x) sample(0:1, 1, prob = c(1-x, x))) # simulate binary observations
} else{
stop("Package \"mgcv\" needed to generate a simulated data set")
}
#> Loading required namespace: mgcv
Y <- as.matrix(Y.all[1:n],nrow = n)
X <- X.all[1:n,]
coords <- coords.all[1:n,]
data <- data.frame(cbind(Y,X))
colnames(data) = c("Y","X1","X2")
mod.glm <- glm(Y~-1+X1+X2,family="binomial",data=data)
mod.glm.esp <- predict(mod.glm,data, type="response")
mod.glm.s2 <- var(Y - mod.glm.esp)
mod.glm.phi <- 0.1*max(dist(coords))
startinit <- c(mod.glm$coef,log(mod.glm.s2),log(mod.glm.phi)) names(startinit) <- c("X1","X2","logsigma2","logphi") ### b. Fit model. result.bin <- fsglmm(Y~-1+X1+X2, kappa=2.5, inits = startinit, data = data,coords = coords, family = "binomial", ntrial = 1, offset = NA,method.optim = "CG", method.integrate = "NR", control = list(maxit=1000,ndeps=rep(1e-2,4),reltol=0.01),rank = 50) result.bin$summary
#> Maximum likelihood estimation
#>
#> Call:
#> bbmle::mle2(minuslogl = nlikSGLMM, start = inits, method = method.optim,
#>     data = list(Y = Y, X = X, mat.dist = mat.dist, ntrial = ntrial,
#>         family = family, method = method.integrate, kappa = kappa,
#>         offset = offset, U1 = U1, rank = rank), vecpar = TRUE,
#>     skip.hessian = TRUE, control = control)
#>
#> Coefficients:
#>           Estimate Std. Error z value Pr(z)
#> X1         1.89994    0.44509      NA    NA
#> X2         0.48805    0.41668      NA    NA
#> logsigma2 -0.29190    0.53395      NA    NA
#> logphi    -1.63619    0.37567      NA    NA
#>
#> -2 log L: 401.2075

### c. Compute predicted random effects.

X.pred <- X.all[(n+1):(n+n.pred),]
coords.pred <- coords.all[(n+1):(n+n.pred),]
pred.sglmm(result.bin,data=X.pred,coords=coords.pred)
#>   [1] 0.7173744 0.9011669 0.8174600 0.7143686 0.6964783 0.9054776 0.7366131
#>   [8] 0.7031458 0.9269817 0.9051901 0.8931835 0.4344861 0.8930961 0.7268835
#>  [15] 0.6794283 0.4742220 0.8197480 0.8347397 0.9134523 0.3202977 0.8738509
#>  [22] 0.6427999 0.7457062 0.5924379 0.9058852 0.8466672 0.6078437 0.7995875
#>  [29] 0.5693686 0.6750436 0.9437651 0.7544525 0.7625658 0.9218002 0.6513120
#>  [36] 0.8064205 0.9295038 0.9182620 0.6105966 0.8203830 0.5307209 0.8577690
#>  [43] 0.7125925 0.7156856 0.8385185 0.9215306 0.9043965 0.6113454 0.6571189
#>  [50] 0.5723820 0.8917714 0.7052825 0.6184694 0.7408175 0.8327119 0.7398930
#>  [57] 0.9234421 0.6905265 0.7738696 0.9242138 0.9136343 0.8689911 0.7954075
#>  [64] 0.4641804 0.8459859 0.4840722 0.8223483 0.5877907 0.8536701 0.8729155
#>  [71] 0.8922757 0.7047227 0.7456766 0.9367009 0.8174338 0.6571998 0.9012756
#>  [78] 0.7099313 0.8327465 0.6874149 0.9082839 0.9439709 0.5455138 0.9463924
#>  [85] 0.9189615 0.5550945 0.8069886 0.9413874 0.5549568 0.7051886 0.9503762
#>  [92] 0.9232483 0.8562200 0.6091460 0.7334006 0.9171320 0.6331991 0.5884551
#>  [99] 0.8798014 0.8443174

# 2. An example for fitting spatial generalized linear mixed models with random projections to negative binomial observations.

### a. Simulate date using parameter values: $$\nu = 2.5, \sigma^2 = 1, \phi = 0.2, \beta = (1,1), \zeta = 2.$$

if(requireNamespace("mgcv")){
sigma2 = 1
phi = 0.2
prec = 2
beta.true = c(1,1)
n = 400
n.pred = 100
coords.all<- matrix(runif((n+n.pred)*2),ncol=2,nrow=n+n.pred) # simulate data locations
X.all <- matrix(runif((n+n.pred)*2),ncol=2,nrow=(n+n.pred))
suppressMessages(require(fields))
dist.all <- fields::rdist(coords.all,coords.all) # compute distance matrix
matern <- function(phi,mat.dist){
K = (1+sqrt(5)/phi*mat.dist+ 5/(3*phi^2)*mat.dist^2)*exp(-sqrt(5)/phi*mat.dist)
return(K)
} # matern 2.5
V.all <- sigma2*matern(phi,dist.all) # compute covariance matrix
set.seed(1)
r.e.all <- mgcv::rmvn(1,rep(0,nrow(coords.all)),V.all) # simulate random effects
mu.all <- exp(X.all%*%beta.true+r.e.all)
Y.all <- rnbinom(length(mu.all), mu=mu.all,size=prec)
} else {
stop("Package \"mgcv\" needed to generate a simulated data set")
}
if(requireNamespace("MASS")){
Y <- as.matrix(Y.all[1:n],nrow = n)
X <- X.all[1:n,]
coords <- coords.all[1:n,]
data <- data.frame(cbind(Y,X))
colnames(data) = c("Y","X1","X2")
mod.glm <- MASS::glm.nb(Y~-1+X1+X2,data=data)
mod.glm.esp <- predict(mod.glm, data)
mod.glm.s2 <- var( log(Y+1) - mod.glm.esp)
mod.glm.phi <- 0.1*max(dist(coords))
mod.glm.prec <- mod.glm$theta startinit <- c(mod.glm$coef,log(mod.glm.s2),log(mod.glm.phi),log(mod.glm.prec))
names(startinit) <- c("X1","X2","logsigma2","logphi","logprec")
} else {
stop("Package \"MASS\" needed to set the initial parameters")
}

### b. Fit model.

result.nb <- fsglmm(Y~-1+X1+X2, kappa=2.5, inits = startinit, data = data,coords = coords, family = "negative.binomial", offset = NA,method.optim = "CG", method.integrate = "NR", control = list(maxit=1000,ndeps=rep(1e-2,5),reltol=0.01),rank = 50)
q = 400
M = eig$vectors[,c(1:q)] Q.s = t(M) %*% Q %*% M tau = 6 Sigma = solve(tau*Q.s) set.seed(1) delta.s = mgcv::rmvn(1, rep(0,q), Sigma) lambda = exp( X%*%beta + M%*%delta.s ) Z = c() for(j in 1:n^2){Z[j] = rpois(1,lambda[j])} Y = as.matrix(Z,ncol=1) data = data.frame("Y"=Y,"X"=X) colnames(data) = c("Y","X1","X2") } else { stop("Packages \"ngspatial\" and \"mgcv\" are needed to generate a simulated data set") } #> Loading required namespace: ngspatial ### b. Fit model linmod <- glm(Y~-1+X1+X2,data=data,family="poisson") # Find starting values linmod$coefficients
#>       X1       X2
#> 1.044368 1.041890
starting <- c(linmod$coefficients,"logtau"=log(1/var(linmod$residuals)) )
result.pois.disc <- fsglmm.discrete(Y~-1+X1+X2, inits = starting, data=data,family="poisson",ntrial=1, method.optim="BFGS", method.integrate="NR", rank=50, A=A)
result.pois.disc\$summary
#> Maximum likelihood estimation
#>
#> Call:
#> bbmle::mle2(minuslogl = nlikSGLMM.discrete, start = inits, method = method.optim,
#>     data = list(Y = Y, X = X, family = family, method = method.integrate,
#>         ntrial, offset = offset, M = M, MQM = MQM, rank = rank),
#>     vecpar = TRUE, skip.hessian = FALSE, control = list(maxit = 1000))
#>
#> Coefficients:
#>        Estimate Std. Error z value Pr(z)
#> X1     0.991924   0.050628      NA    NA
#> X2     1.031817   0.050312      NA    NA
#> logtau 1.724322   0.346589      NA    NA
#>
#> -2 log L: 3480.376