Vignette_FLaplace

Sangwan Lee, Jaewoo Park

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)
result.nb$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         0.71409    0.18226      NA    NA
#> X2         0.79465    0.17990      NA    NA
#> logsigma2  0.16537    0.41144      NA    NA
#> logphi    -1.48301    0.20974      NA    NA
#> logprec    0.61702    0.11935      NA    NA
#> 
#> -2 log L: 1973.482

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.nb,data=X.pred,coords=coords.pred)
#>   [1]  1.0796784  2.5511352  6.2148177  6.4954555  1.6994834  1.1658579
#>   [7]  3.4941327  2.7992210  1.9913884 13.2394374 11.4312459  9.5249584
#>  [13]  5.6843326  3.2474477  1.3924241 11.1788155  6.3858692 15.6510471
#>  [19] 30.5209369 32.5874261 33.8961389  2.3121162  1.0648111  1.4463879
#>  [25]  1.6027887 15.2776189  1.6313418  4.2806683  7.8745649  0.4147299
#>  [31]  3.3787403 26.5047551  2.2723055  3.3217126  2.9620488  3.7096393
#>  [37]  2.1696730  1.0133735 11.9323421  6.8693548  2.8557717  1.1899519
#>  [43]  3.1238762  6.8792812  5.6720593  7.1592524  2.3571801  1.1341541
#>  [49] 23.1115252  0.7456393 19.5059490  1.2288857  9.0283915  4.2649047
#>  [55]  1.4327406  1.2139899  1.9578175  6.5373381  0.8973056  8.8158868
#>  [61]  2.1944846  4.5888146  2.2190761  0.7039620  2.2549576  1.4654028
#>  [67]  3.0896944  0.8932928  3.5688911  0.9714466  3.9925341  3.0443346
#>  [73] 12.2336647  4.8261169  0.9185217  4.7886871  0.4901198  3.0649941
#>  [79] 14.9865994  2.4464208  9.9471592  1.5416763  1.2350628  1.6442479
#>  [85]  0.5339071  2.8120669  9.0368075  1.5598224  3.3950542  4.6703657
#>  [91]  1.1676898  3.2353459  2.0608186  1.5531116 11.2147363 18.4602356
#>  [97]  0.8459434  3.4027514  3.8329304 12.2404823

3. An example for fitting spatial generalized linear mixed models with random projections to poisson observations (discrete spatial domain).

a. Simulate date using parameter values: $ = (1,1), = 6.$

if(requireNamespace("ngspatial")&
   requireNamespace("mgcv")){
  n = 30
  A = ngspatial::adjacency.matrix(n)
  Q = diag(rowSums(A),n^2) - A
  x = rep(0:(n - 1) / (n - 1), times = n) 
  y = rep(0:(n - 1) / (n - 1), each = n) 
  X = cbind(x, y)                                 # Use the vertex locations as spatial covariates.
  beta = c(1, 1)                                  # These are the regression coefficients.
  P.perp = diag(1,n^2) - X%*%solve(t(X)%*%X)%*%t(X)
  eig = eigen(P.perp %*% A %*% P.perp)
  eigenvalues = eig$values
  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