Nothing
##############################################################################################
inv.logit <- function(x) {
##############################################################################################
res = 1/(1+exp(-x))
return(res)
}
##############################################################################################
logit <- function(x) {
##############################################################################################
res = log(x) - log(1-x)
return(res)
}
##############################################################################################
matern.cov <- function(phi,kappa,mat.dist) {
##############################################################################################
if(phi < 1e-8 ) phi = 1e-8
if(kappa == 1/2) {
K = exp(-1/phi*mat.dist)
}
if(kappa == 1.5) {
K = (1+sqrt(3)/phi*mat.dist)*exp(-sqrt(3)/phi*mat.dist)
}
if(kappa == 2.5){
K = (1+sqrt(5)/phi*mat.dist+ 5/(3*phi^2)*mat.dist^2)*exp(-sqrt(5)/phi*mat.dist)
}
if(kappa == 10){
K = exp(-1/(2*phi^2)*mat.dist^2)
}
return(K)
}
##############################################################################################
rp_fl <- function(n=1000,k=100,C,alpha=1,m=50){
##############################################################################################
PHI = C^alpha %*% matrix(stats::rnorm(n*k,0,sd=k^-0.25),ncol=k)
K11 = t(PHI)%*%C%*%PHI
SVDK11 = svd(K11)
V11 = SVDK11$u
L11 = SVDK11$d
C = C%*%PHI%*%V11%*%diag(L11^-0.5)
C.svd = svd(C)
u_m = C.svd$u[,1:m]
d_m = C.svd$d[1:m]
return(list("u" = u_m, "d" = as.matrix(d_m,ncol=1)))
}
##############################################################################################
sigma.rp <- function(phi,kappa,mat.dist,U1,rank){
##############################################################################################
covmat_original <- matern.cov(phi=phi,kappa=kappa,mat.dist=mat.dist)
covmat_rp <- rp_fl(n=dim(mat.dist)[1],k=2*rank,alpha=1,C=covmat_original,m=rank)
m <- rank
um <- covmat_rp$u[,1:m]
signdiag = sign(diag(t(um)%*%U1))
signdiag = as.logical(1-signdiag)
um[,signdiag] = -um[,signdiag]
dm <- diag(covmat_rp$d[1:m,]^2)
return(list(um=um ,dm = dm))
}
##############################################################################################
sigma.spectra <- function(phi,kappa,mat.dist,U1,rank){
##############################################################################################
covmat_original <- matern.cov(phi=phi,kappa=kappa,mat.dist=mat.dist)
covmat_svd <- RSpectra::eigs_sym(covmat_original,k=rank,which="LA")
m <- rank
um <- covmat_svd$vectors[,1:m]
signdiag = sign(diag(t(um)%*%U1))
signdiag = as.logical(1-signdiag)
um[,signdiag] = -um[,signdiag]
dm <- diag(covmat_svd$values[1:m])
return(list(um=um ,dm = dm))
}
##############################################################################################
Q.b.poisson <- function(b,Xbeta,Y, det.Sigma, inv.Sigma,um,dm){
##############################################################################################
eta = Xbeta + um%*%dm^0.5%*%b
dens = sum(stats::dpois(Y,lambda= exp(eta), log=TRUE)) +
my.gauss(b,det.Sigma = det.Sigma, inv.Sigma = inv.Sigma)
return(dens)
}
##############################################################################################
Q1.b.poisson <- function(b,Xbeta,Y,det.Sigma, inv.Sigma,um,dm){
##############################################################################################
grad <- t(dm^0.5%*%t(um)%*%Y) - t(exp(Xbeta + um%*%dm^0.5%*%b))%*%um%*%dm^0.5 - t(b)*inv.Sigma
return(grad)
}
##############################################################################################
Q2.b.poisson <- function(b, Xbeta,Y,det.Sigma,inv.Sigma,um,dm){
##############################################################################################
m = dim(dm)[1]
eta = exp(Xbeta + um%*%dm^0.5%*%b)
nHess = dm^0.5%*%t(um)%*%diag(as.numeric(eta))%*%um%*%dm^0.5 + diag(inv.Sigma,m)
Hess <- -nHess
return(Hess)
}
##############################################################################################
Q.b.binomial <- function(b,Xbeta,Y,det.Sigma,inv.Sigma,ntrial=1,um,dm){
##############################################################################################
eta = Xbeta + um%*%dm^0.5%*%b
p <- inv.logit(eta)
dens = sum(stats::dbinom(Y,size = ntrial, prob = p,log=TRUE)) + my.gauss(b,det.Sigma = det.Sigma, inv.Sigma = inv.Sigma)
return(dens)
}
##############################################################################################
Q1.b.binomial <- function(b,Xbeta,Y,det.Sigma, inv.Sigma,ntrial=1,um,dm){
##############################################################################################
eta = Xbeta + um%*%dm^0.5%*%b
b1 <- ntrial*inv.logit(eta)
grad <- t(dm^0.5%*%t(um)%*%Y) - t(b1)%*%um%*%dm^0.5 - t(b)*inv.Sigma
return(grad)
}
##############################################################################################
Q2.b.binomial <- function(b, Xbeta,Y,det.Sigma,inv.Sigma,ntrial=1,um,dm){
##############################################################################################
m = dim(dm)[1]
eta = as.numeric(Xbeta + um%*%dm^0.5%*%b)
nHess = ntrial*dm^0.5%*%t(um)%*%diag( exp(eta)/(1+exp(eta))^2 )%*%um%*%dm^0.5 + diag(inv.Sigma,m)
Hess <- -nHess
return(Hess)
}
##############################################################################################
Q.b.negbin <- function(b,Xbeta,Y,prec,det.Sigma,inv.Sigma,um,dm){
##############################################################################################
eta = Xbeta + um%*%dm^0.5%*%b
dens = sum(stats::dnbinom(Y,size = prec, mu = exp(eta),log=TRUE)) + my.gauss(b,det.Sigma = det.Sigma, inv.Sigma = inv.Sigma)
return(dens)
}
##############################################################################################
Q1.b.negbin <- function(b,Xbeta,Y,prec,det.Sigma,inv.Sigma,um,dm){
##############################################################################################
et1 = exp(Xbeta)
et2 = exp(um%*%dm^0.5%*%b)
et12 = et1*et2
grad <- t(dm^0.5%*%t(um)%*%Y) - t((prec+Y)*(et12*(et12+prec)^-1))%*%um%*%dm^0.5 - t(b)*inv.Sigma
return(grad)
}
##############################################################################################
Q2.b.negbin <- function(b, Xbeta,Y,prec,det.Sigma,inv.Sigma,um,dm){
##############################################################################################
m = dim(dm)[1]
et1 = exp(Xbeta)
et2 = exp(um%*%dm^0.5%*%b)
et12 = et1*et2
p1 <- et12*((et12+prec)^-1)
p2 <- p1^2
D <- p1-p2
nHess = dm^0.5%*%t(um)%*%diag(as.numeric((prec+Y)*D))%*%um%*%dm^0.5 + diag(inv.Sigma,m)
Hess <- -nHess
return(Hess)
}
##############################################################################################
my.gauss <- function(b, det.Sigma, inv.Sigma){
##############################################################################################
n <- length(b)
dens <- (-n/2) * log(2*pi) -(1/2)*det.Sigma - (1/2)*inv.Sigma * t(b) %*% b
return(dens)
}
##############################################################################################
newton_raphson <- function(initial, score, hessian, tolerance = 0.001, max.iter, 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 <- function(Q.b, gr, hess, method,m.dim=50, ...){
##############################################################################################
log.integral <- -sqrt(.Machine$double.xmax)
Init <- rep(0,m.dim)
temp <- try(newton_raphson(initial = Init, score=gr, hessian = hess,m.dim = m.dim, max.iter = 100,...), silent = TRUE)
if(class(temp)!='try-error'){
value <- Q.b(b = temp[2][[1]] , ...)
log.integral <- value + ((m.dim/2)*log(2*pi) - 0.5*determinant(-temp[1][[1]])$modulus)}
return(log.integral)
}
##############################################################################################
nlikSGLMM <- function(par, Y, X, kappa, mat.dist, family, method,
ntrial=1, offset=NA,U1,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)}
sigma <- exp(as.numeric(par[c(n.beta+1)]))
if(sigma < 1e-8) sigma <- 1e-8
phi <- exp(as.numeric(par[c(n.beta+2)]))
if(family == "negative.binomial") { prec <- exp(as.numeric(par[c(n.beta+3)]))}
Sigma <- sigma.rp(phi=phi,kappa=kappa,mat.dist=mat.dist,U1=U1,rank=rank)
um <- Sigma$um
dm <- Sigma$dm
m <- rank
det.Sigma <- m * log(sigma)
inv.Sigma <- 1 / sigma
if(family == "poisson"){
I <- fLaplace(Q.b.poisson, gr = Q1.b.poisson, hess = Q2.b.poisson, method = "NR",
m.dim = m,um = um,dm = dm,
Xbeta = Xbeta, Y = Y, det.Sigma = det.Sigma, inv.Sigma = inv.Sigma)}
if(family == "binomial"){
I <- fLaplace(Q.b.binomial, gr = Q1.b.binomial, hess = Q2.b.binomial,method = "NR",
m.dim = m, um = um, dm = dm,
Xbeta = Xbeta, Y = Y, det.Sigma = det.Sigma, inv.Sigma = inv.Sigma,
ntrial = ntrial)}
if(family == "negative.binomial"){
I <- fLaplace(Q.b.negbin, gr = Q1.b.negbin, hess = Q2.b.negbin, method = "NR",
m.dim = m, um = um, dm = dm,
Xbeta = Xbeta, prec = prec,Y = Y, det.Sigma = det.Sigma, inv.Sigma = inv.Sigma)}
return(-I)
}
##############################################################################################
nlikSGLMM_pca <- function(par,rank,family,method,kappa,X,Y,mat.dist,ntrial,U1,offset){
##############################################################################################
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)}
sigma <- as.numeric(exp(par[c(n.beta+1)]))
phi <- as.numeric(exp(par[c(n.beta+2)]))
if(family == "negative.binomial") { prec <- exp(as.numeric(par[c(n.beta+3)]))}
Sigma <- sigma.spectra(phi = phi,kappa = kappa,mat.dist = mat.dist,U1=U1,rank=rank)
um <- Sigma$um
dm <- Sigma$dm
m <- rank
det.Sigma <- m * log(sigma)
inv.Sigma <- 1 / sigma
if(family == "poisson"){
I <- fLaplace(Q.b.poisson, gr = Q1.b.poisson, hess = Q2.b.poisson, method = method,
m.dim = m,um = um,dm = dm,
Xbeta = Xbeta, Y = Y, det.Sigma = det.Sigma, inv.Sigma = inv.Sigma)}
if(family == "binomial"){
I <- fLaplace(Q.b.binomial, gr = Q1.b.binomial, hess = Q2.b.binomial,method = method,
m.dim = m, um = um, dm = dm,
Xbeta = Xbeta, Y = Y, det.Sigma = det.Sigma, inv.Sigma = inv.Sigma,
ntrial= ntrial)}
if(family == "negative.binomial"){
I <- fLaplace(Q.b.negbin, gr = Q1.b.negbin, hess = Q2.b.negbin, method = method,
m.dim = m, um = um, dm = dm,
Xbeta = Xbeta, prec = prec,Y = Y, det.Sigma = det.Sigma, inv.Sigma = inv.Sigma)}
return(-I)
}
#' @title Fitting Projection Based Laplace Approximation for Spatial Generalized Linear Mixed Model
#'
#' @description \code{fsglmm} is used to fit reduced-dimensional spatial generalized linear mixed models for continuous spatial domain.
#'
#' @param formula an object of class "formula."
#' @param kappa the smoothness parameter for the matern process (either 0.5, 1.5, 2.5 or 10).
#' @param inits starting values for the parameters.
#' @param data a data frame containing variables in the model.
#' @param coords a matrix of dimension \eqn{N} x 2 representing the longitude and latitude of each observation.
#' @param family a character string of the error distribution and link function to be used in the model.
#' @param ntrial a numeric vector for a 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 control a list of control parameters.
#' @param rank an integer of 'rank' to be used for projections. Default is 5 percent of observations.
#'
#'
#' @return a list containing the following components:
#' @return \code{summary} a summary of the fitted model
#' @return \code{vcov} a symmetric matrix giving an estimate of the Hessian at the solution found.
#' @return \code{mle2} an object of class "mle2"
#' @return \code{family} the family used.
#' @return \code{kappa} the matern smoothness parameter used.
#' @return \code{Delta} a matrix containing the estimated random effects of the reduced dimensional model.
#' @return \code{U} a matrix whose columns contain the estimated eigenvectors of the reduced dimensional model.
#' @return \code{D} a matrix whose diagonal components contain the estimated eigenvalues of the reduced dimensional model.
#' @return \code{coords} the matrix of coordinates used.
#'
#' @import bbmle
#' @importFrom fields rdist
#' @export
#'
#' @references Jaewoo Park and Sangwan Lee - "A Projection-based Laplace Approximation for Spatial Latent Variable Models"
#'
#' @examples
#' \donttest{
#' 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)
#' X.all <- matrix(runif((n+n.pred)*2),ncol=2,nrow=(n+n.pred))
#' dist.all <- fields::rdist(coords.all,coords.all)
#' V.all <- sigma2*(1+sqrt(5)/phi*dist.all+5/(3*phi^2)*dist.all^2)*exp(-sqrt(5)/phi*dist.all)
#' set.seed(1)
#' r.e.all <- mgcv::rmvn(1,rep(0,nrow(coords.all)),V.all)
#' pi.all <- X.all%*%beta.true + r.e.all
#' p.all <- exp(pi.all)/(1+exp(pi.all))
#' Y.all <- sapply(p.all, function(x) sample(0:1, 1, prob = c(1-x, x)))
#' 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")
#' 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",rank = 50)
#' }
#' }
fsglmm <- function(formula, kappa ,inits, data ,coords, family, ntrial=1,offset=NA,
method.optim, method.integrate,rank=NULL, control = list()){
formula <- stats::as.formula(formula)
model.formula <- stats::model.frame(formula, data=data)
Y <- stats::model.response(model.formula)
X <- stats::model.matrix(formula, data=data)
mat.dist <- as.matrix(fields::rdist(coords,coords))
names.par <- c(colnames(X), "logsigma2", "logphi")
n.beta <- dim(X)[2]
n <- dim(X)[1]
if(is.null(rank)) rank <- floor(n*0.05)
sign.cov <- as.matrix(matern.cov(phi=exp(inits['logphi']),kappa=kappa,mat.dist=mat.dist))
sign.temp <- rp_fl(n = n,k=2*rank,alpha=1,C = sign.cov,m=rank)
U1 <- sign.temp$u[,1:rank]
if(family == "negative.binomial"){names.par <- c(names.par,"logprec")}
if(length(control)==0) control = list(maxit=1000,ndeps=rep(1e-3,length(names.par)),reltol=1e-3)
names(inits) <- parnames(nlikSGLMM) <- names.par
estimate <- bbmle::mle2(nlikSGLMM, start=inits,
vecpar=TRUE,
method=method.optim,
control = control,
skip.hessian = TRUE,
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))
coef_final <- estimate@coef
output.hessian.pca <- stats::optimHess(coef_final,fn=nlikSGLMM_pca,rank=rank,kappa=kappa,family=family,
method=method.integrate,X=X,Y=Y,mat.dist=mat.dist,ntrial=ntrial,U1=U1,
offset = offset)
vcov.mat <- try(solve(output.hessian.pca),silent = T)
if(class(vcov.mat)[1] != "try-error"){
output.se <- suppressWarnings(sqrt(diag(solve(output.hessian.pca))))
if(sum(is.na(output.se)) != 0 ) {
output.se <- NA
}
} else{
vcov.mat <- matrix(NA,ncol=length(names.par),nrow=length(names.par))
}
n.pars <- length(coef(estimate))
summary.estimate <- summary(estimate)
summary.estimate@coef[,1] <- summary.estimate@coef[,1]
summary.estimate@coef[,2] <- output.se
summary.estimate@coef[,3] <- NA
summary.estimate@coef[,4] <- NA
output <- list()
output[1][[1]] <- summary.estimate
output[2][[1]] <- vcov.mat
output[3][[1]] <- estimate
sigma <- as.numeric(exp(coef_final["logsigma2"]))
phi <- as.numeric(exp(coef_final["logphi"]))
Sigma <- sigma.rp(phi=phi,kappa=kappa,mat.dist=mat.dist,U1=U1,rank=rank)
UM <- Sigma$um
DM <- Sigma$dm
m <- rank
det.Sigma <- m * log(sigma)
inv.Sigma <- 1 / sigma
if(family == "poisson"){
Xbeta <- X%*%coef_final[1:n.beta]
gr <- Q1.b.poisson
hess <- Q2.b.poisson
DELTA.HAT <- newton_raphson(rep(0,m),score = gr,hessian = hess,m.dim=m,max.iter = 100,
Xbeta=Xbeta, Y=Y, det.Sigma=det.Sigma, inv.Sigma=inv.Sigma,
um=UM, dm=DM)[2][[1]]}
if(family == "binomial"){
Xbeta = X%*%coef_final[1:n.beta]
gr <- Q1.b.binomial
hess <- Q2.b.binomial
DELTA.HAT <- newton_raphson(rep(0,m),score = gr,hessian = hess,m.dim=m,max.iter = 100,
Xbeta=Xbeta, Y=Y, det.Sigma=det.Sigma, inv.Sigma=inv.Sigma,
ntrial=ntrial,um=UM, dm=DM)[2][[1]]}
if(family == "negative.binomial"){
prec <- exp(coef_final["logprec"])
Xbeta <- X%*%coef_final[1:n.beta]
gr <- Q1.b.negbin
hess <- Q2.b.negbin
DELTA.HAT <- newton_raphson(rep(0,m),score = gr,hessian = hess,m.dim=m,max.iter = 100,
Xbeta=Xbeta, Y=Y, prec = prec,det.Sigma=det.Sigma, inv.Sigma=inv.Sigma,
um=UM, dm=DM)[2][[1]]}
output[4][[1]] <- family
output[5][[1]] <- kappa
output[6][[1]] <- as.matrix(DELTA.HAT)
output[7][[1]] <- UM
output[8][[1]] <- DM
output[9][[1]] <- coords
names(output) <- c("summary","vcov","mle2","family","kappa","Delta","U","D","coords")
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.