Nothing
##############################################################################################
Q.d.poisson <- function(d,Xbeta,Y,det.Sigma,inv.Sigma,M){
##############################################################################################
eta = as.numeric(Xbeta + M%*%d)
dens = sum(stats::dpois(Y,exp(eta),log=TRUE)) +
my.gauss.d(d,det.Sigma = det.Sigma, inv.Sigma = inv.Sigma)
return(dens)
}
##############################################################################################
Q1.d.poisson <- function(d,Xbeta,Y,det.Sigma, inv.Sigma,M){
##############################################################################################
eta = as.numeric(Xbeta + M%*%d)
grad <- t(t(M)%*%Y) - t(exp(eta))%*%M - t(d)%*%inv.Sigma
return(grad)
}
##############################################################################################
Q2.d.poisson <- function(d, Xbeta,Y,det.Sigma,inv.Sigma,M){
##############################################################################################
eta = as.numeric(Xbeta + M%*%d)
nHess = t(M)%*%diag(exp(eta))%*%M + inv.Sigma
Hess <- -nHess
return(Hess)
}
##############################################################################################
my.gauss.d <- function(d, det.Sigma, inv.Sigma){
##############################################################################################
n <- length(d)
dens <- (-n/2) * log(2*pi) + (1/2)*det.Sigma - (1/2) * t(d) %*% inv.Sigma %*% d
return(dens)
}
##############################################################################################
newton_raphson.discrete <- function(initial, score, hessian, tolerance = 0.001, max.iter=100, m.dim,...){
##############################################################################################
solution <- matrix(NA, max.iter, m.dim)
solution[1,] <- initial
for(i in 2:max.iter){
HSS <- hessian(initial, ...)
SCR <- t(score(initial, ...))
solution[i,] <- initial - solve(HSS, SCR)
initial <- solution[i,]
convergence <- abs(solution[i,] - solution[i-1,])
if(all(convergence < tolerance)==TRUE) break
}
output <- list()
output[1][[1]] <- HSS
output[2][[1]] <- initial
return(output)
}
##############################################################################################
fLaplace.discrete <- function(Q.b, gr, hess, method, optimizer, m.dim, ...){
##############################################################################################
log.integral <- -sqrt(.Machine$double.xmax)
Init <- rep(0,m.dim)
if (method=="NR"){
temp <- try(newton_raphson.discrete(initial = Init, score=gr, hessian = hess, m.dim = m.dim, max.iter = 100,...), silent = TRUE)}
if(class(temp)!='try-error' & method=="NR"){
value <- Q.b(d = temp[2][[1]] , ...)
log.integral <- value + ((m.dim/2)*log(2*pi) - 0.5*determinant(-temp[1][[1]])$modulus)}
return(log.integral)
}
##############################################################################################
nlikSGLMM.discrete <- function(par, Y, X, nugget, family, method, ntrial=1, offset=NA,
M,Q,MQM,rank){
##############################################################################################
I = -sqrt(.Machine$double.xmax)
n <- length(Y)
n.beta <- dim(X)[2]
beta <- as.numeric(par[1:n.beta])
Xbeta <- X%*%beta
if(is.na(offset)[1] != TRUE){Xbeta <- cbind(X,log(offset))%*%c(beta,1)}
tau <- as.numeric(exp(par[(n.beta+1)]))
det.Sigma <- rank*log(tau)
inv.Sigma <- tau * MQM
if(family == "poisson"){
I <- fLaplace.discrete(Q.d.poisson, gr = Q1.d.poisson, hess = Q2.d.poisson, method= method, optimizer="BFGS",
m.dim = rank, M = M,
Xbeta = Xbeta, Y = Y, det.Sigma = det.Sigma, inv.Sigma = inv.Sigma)}
return(-I)
}
#' Fitting Projection Based Laplace Approximation for Spatial Generalized Linear Mixed Model
#' @description \code{fsglmm.discrete}
#' @param formula an object of class "formula."
#' @param inits starting values for the parameters.
#' @param data a data frame containing variables in the model.
#' @param family a character string of the error distribution and link function to be used in the model.
#' @param ntrial a numeric vector for binomial model.
#' @param offset this is used to specify an a priori a known component to be included in the linear predictor during fitting.
#' @param method.optim the method to be used for outer optimization. "CG" for Conjugate Gradient Method.
#' @param method.integrate the method to be used for inner optimization. "NR" for Newton Raphson Method.
#' @param rank an integer of 'rank' to be used for projections. Default is 5 percent of observations.
#' @param A an adjacency matrix
#'
#'
#'
#' @return a list containing the following components:
#' @return \code{summary} a summary of the fitted model
#' @return \code{mle2} an object of class "mle2"
#' @return \code{Delta} a matrix containing the estimated random effects of the reduced dimensional model.
#' @return \code{M} the projection matrix used.
#'
#'
#' @import bbmle
#' @export
#'
#' @examples
#' \donttest{
#' 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)
#' }
#'
#' }
#'
fsglmm.discrete <- function(formula, inits, data, family, ntrial=1,
method.optim, method.integrate, rank=NULL, A, offset = NA){
formula <- stats::as.formula(formula)
mf <- stats::model.frame(formula, data=data)
Y <- stats::model.response(mf)
X <- stats::model.matrix(formula, data=data)
names <- c(colnames(X), "logtau")
n.beta <- dim(X)[2]
n.obs <- dim(X)[1]
if(is.null(rank)) rank <- floor(n.obs*0.05)
P.perp = diag(1,n.obs) - X%*%solve(t(X)%*%X)%*%t(X)
eig = RSpectra::eigs_sym(P.perp %*% A %*% P.perp,k=rank,which="LA")
eigenvalues = eig$values
M = eig$vectors[,c(1:rank)]
Q = diag(rowSums(A),n.obs) - A
MQM = t(M) %*% Q %*% M
names(inits) <- parnames(nlikSGLMM.discrete) <- names
estimate <- bbmle::mle2(nlikSGLMM.discrete, start=inits,
vecpar=TRUE,
method=method.optim,
control = list(maxit=1000),
skip.hessian = FALSE,
data = list(Y=Y, X=X, family= family, method = method.integrate, ntrial,
offset = offset, M = M, MQM = MQM, rank = rank))
n.pars <- length(coef(estimate))
summary.estimate <- summary(estimate)
summary.estimate@coef[,1] <- summary.estimate@coef[,1]
output.se <- suppressWarnings(sqrt(diag(estimate@vcov)))
if(sum(is.na(output.se)) != 0 ) {
summary.estimate@coef[,2] <- NA
}
summary.estimate@coef[,3] <- NA
summary.estimate@coef[,4] <- NA
output <- list()
output[1][[1]] <- summary.estimate
output[2][[1]] <- estimate
coef_final <- summary.estimate@coef[,1]
beta <- as.numeric(coef_final[1:n.beta])
tau <- as.numeric(exp(coef_final[n.beta+1]))
det.Sigma <- rank*log(tau)
inv.Sigma <- tau * MQM
if(family == "poisson"){
Xbeta = X%*%beta[1:n.beta]
if(is.na(offset)[1] != TRUE) { Xbeta <- cbind(X,log(offset))%*%c(beta,1)}
Init <- rep(0,rank)
gr = Q1.d.poisson
hess = Q2.d.poisson
DELTA.HAT <- newton_raphson.discrete(initial = Init, score = gr, hessian = hess,
tolerance = 0.001, max.iter = 100,
m.dim=rank,
Y=Y, Xbeta = Xbeta, det.Sigma=det.Sigma,inv.Sigma=inv.Sigma,M=M)[2][[1]]
}
output[3][[1]] <- as.matrix(DELTA.HAT)
output[4][[1]] <- M
names(output) <- c("summary","mle2","Delta","M")
return(output)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.