The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.