fastLaplace (version 0.0.2)

fsglmm.discrete: Fitting Projection Based Laplace Approximation for Spatial Generalized Linear Mixed Model

Description

fsglmm.discrete

Usage

fsglmm.discrete(
  formula,
  inits,
  data,
  family,
  ntrial = 1,
  method.optim,
  method.integrate,
  rank = NULL,
  A,
  offset = NA
)

Arguments

formula

an object of class "formula."

inits

starting values for the parameters.

data

a data frame containing variables in the model.

family

a character string of the error distribution and link function to be used in the model.

ntrial

a numeric vector for binomial model.

method.optim

the method to be used for outer optimization. "CG" for Conjugate Gradient Method.

method.integrate

the method to be used for inner optimization. "NR" for Newton Raphson Method.

rank

an integer of 'rank' to be used for projections. Default is 5 percent of observations.

A

an adjacency matrix

offset

this is used to specify an a priori a known component to be included in the linear predictor during fitting.

Value

a list containing the following components:

summary a summary of the fitted model

mle2 an object of class "mle2"

Delta a matrix containing the estimated random effects of the reduced dimensional model.

M the projection matrix used.

Examples

# NOT RUN {
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)
beta = c(1, 1)
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")
linmod <- glm(Y~-1+X1+X2,data=data,family="poisson") # Find starting values
linmod$coefficients
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)
}

# }
# NOT RUN {
# }