spMvGLM(spBayes)
spMvGLM()所属R语言包：spBayes
Function for fitting multivariate Bayesian generalized linear spatial regression models
函数拟合多元贝叶斯广义线性空间回归模型
译者：生物统计家园网 机器人LoveR
描述----------Description----------
The function spMvGLM fits multivariate Bayesian generalized linear spatial regression models. Given a set of knots, spMvGLM will also fit a predictive process model (see references below).
的功能spMvGLM适合多元贝叶斯广义线性空间回归模型。给定一组的结，spMvGLM也符合预测过程模型（请参阅下面的参考资料）。
用法----------Usage----------
spMvGLM(formula, family="binomial", weights, data = parent.frame(), coords, knots,
amcmc, starting, tuning, priors, cov.model, n.samples, sub.samples,
verbose=TRUE, n.report=100, ...)
参数----------Arguments----------
参数：formula
for a multivariate model with q response variables, this is a list of univariate formulas. See example below.
q响应变量的多变量模型，这是一个单变量的公式列表。见下面的例子。
参数：family
currently only supports binomial and poisson data using the logit and log link functions, respectively.
目前只支持binomial和poisson数据的Logit和log链接功能，分别使用。
参数：weights
an optional n x q matrix of weights to be used in the fitting process when family is binomial. Here the order of the columns correspond to the univariate models in the formula list. Unless otherwise specified a matrix of ones is assumed i.e., one trial at each location.
一个可选的n x q的权重矩阵中要使用的装配过程中，当family是binomial。这里列的顺序对应的公式列表中的单变量模型。除另有规定外矩阵的假设，一项试验，即在每个位置。
参数：data
an optional data frame containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which spMvGLM is called.
一个可选的数据框包含在模型中的变量。如果没有找到，数据，变量environment(formula)，通常是spMvGLM被称为环境。
参数：coords
an n x 2 matrix of the observation coordinates in R^2 (e.g., easting and northing).
n x 2R^2（例如，东向和北向）的观察坐标矩阵。
参数：knots
either a m x 2 matrix of the predictive process knot coordinates in R^2 (e.g., easting and northing) or a vector of length two or three with the first and second elements recording the number of columns and rows in the desired knot grid. The third, optional, element sets the offset of the outermost knots from the extent of the coords extent.
m x 2矩阵的预测过程结坐标在R^2（例如，东向和北向）或长度的向量的记录的数目的列和行，在该第一和第二元件中的两个或三个与所需的结栅格。可选的第三个，元素设置从coords程度的程度的最外节的偏移量。
参数：amcmc
a list with tags n.batch, batch.length, and accept.rate.
标签的列表n.batch，batch.length和accept.rate。
参数：starting
a list with each tag corresponding to a parameter name. Valid list tags are beta, A, phi, and nu. The value portion of each tag is a vector of parameter's starting value. For A the vector is of length (q^2-q)/2+q and phi and nu are of length q. Here, A holds the the lower-triangle elements in column major ordering of the Cholesky square root of the spatial cross-covariance matrix. If the predictive process is used then w must be of length qm; otherwise, it must be of length qn. Alternatively, w can be set as a scalar, in which case the value is repeated.
与每个标记对应的参数名称的列表。有效的列表标签是beta，A，phi，nu。每个标签的价值部分是一个向量参数的初始值。对于A向量的长度为(q^2-q)/2+q和phi和nu是长度q的。在这里，A保持较低三角形的元素在列主要顺序的空间的互协方差矩阵的Cholesky平方根。如果预测过程中被使用，则w必须的长度qm，否则，它必须是长度qn。还有一种方法，w可以被设置为一个标量，在这种情况下，该值是重复的。
参数：tuning
a list with each tag corresponding to a parameter name. Valid list tags are beta, A, phi, nu, and w. The value portion of each tag defines the variance of the Metropolis normal proposal distribution. The tuning value for beta can be a vector of length p or the lower-triangle of the pxp Cholesky square-root of the desired proposal variance matrix. For A, the vector is of length (q^2-q)/2+q and phi and nu are of length q. If the predictive process is used then w must be of length m; otherwise, it must be of length n. Alternatively, w can be set as a scalar, in which the value is repeated.
与每个标记对应的参数名称的列表。有效的列表标签是beta，A，phi，nu，w。每个标签都会正常的提议分布的方差定义的值部分。 beta的调谐值可以是一个矢量长度p或pxp的Cholesky平方根所需的建议方差矩阵的下部三角。对于A，矢量长度(q^2-q)/2+q和phi和nu是长度为q的。如果预测过程中被使用，则w必须的长度m，否则，它必须是长度n。或者，w可以被设置为一个标量，在该值是重复的。
参数：priors
a list with each tag corresponding to a parameter name. Valid list tags are beta.flat, beta.normal, K.IW, Psi.IW, phi.unif, and nu.unif. If beta.normal then covariate specific mean and variance hyperparameters are passed as the first and second list elements, respectively. K is assumed to follow an inverse-Wishart distribution, whereas the spatial decay phi and smoothness nu parameters are assumed to follow Uniform distributions. The hyperparameters of the inverse-Wishart are passed as a list of length two, with the first and second elements corresponding to the df and qxq scale matrix, respectively. The hyperparameters of the Uniform are also passed as a vector of length 2xq with consecutive elements representing the first and second elements corresponding to the lower and upper support in the order of the univariate models given in formula.
与每个标记对应的参数名称的列表。有效的列表标签是beta.flat，beta.normal，K.IW，Psi.IW，phi.unif，nu.unif。如果beta.normal协变量的均值和方差的超传递的第一和第二个列表元素，分别。 K假设遵循逆Wishart分布，而在空间衰减phi和平滑nu参数的假设为均匀分布。逆威沙特的超长度2作为一个列表被传递，与df和qxq规模的矩阵对应的第一和第二元件，分别。作为一个向量的长度也被传递均匀的超2xq代表在formula给出的单变量模型的顺序中的下限和上限的支持对应的第一和第二元件的连续的元素。
参数：cov.model
a quoted key word that specifies the covariance function used to model the spatial dependence structure among the observations. Supported covariance model key words are: "exponential", "matern", "spherical", and "gaussian". See below for details.
一个带引号的关键字，指定的协方差函数用于模拟空间之间的依赖结构的观察。支持的协方差模型的关键词是："exponential"，"matern"，"spherical"和"gaussian"。有关详细信息，请参见下文。
参数：n.samples
the number of MCMC iterations.
MCMC中的迭代数。
参数：sub.samples
a vector of length 3 that specifies start, end, and thin, respectively, of the MCMC samples. The default is c(1, n.samples, 1) (i.e., all samples).
一个向量的长度为3，指定起点，终点，又瘦，分别MCMC样本。默认值是c(1, n.samples, 1)（即，所有样品）。
参数：verbose
if TRUE, model specification and progress of the sampler is printed to the screen. Otherwise, nothing is printed to the screen.
如果TRUE，型号规格和进步的取样打印到屏幕上。否则，一切都被打印到屏幕上。
参数：n.report
the interval to report Metropolis acceptance and MCMC progress.
报告的时间间隔都会接受和MCMC进步。
参数：...
currently no additional arguments.
目前没有任何额外的参数。
Details
详细信息----------Details----------
If a binomial model is specified the response vector is the number of successful trials at each location. Note, this is different from the glm argument which is passed as a rate.
如果binomial模型指定的响应向量的数量在每个位置的成功试验。请注意，这是从glm参数传递的速度不同。
值----------Value----------
An object of class spMvGLM, which is a list with the following tags:
类的一个对象spMvGLM，这是一个与下面的标签列表：
参数：coords
the n x 2 matrix specified by coords.
n x 2的矩阵指定coords。
参数：knot.coords
the m x 2 matrix as specified by knots.
m x 2矩阵指定的knots。
参数：p.samples
a coda object of posterior samples for the defined parameters.
一个coda对象定义的参数后样品。
参数：acceptance
the Metropolis sampling acceptance rate. If amcmc is used then this will be a matrix of each parameter's acceptance rate at the end of each batch. Otherwise, the sampler is a Metropolis with a joint proposal of all parameters, and the acceptance rate is the average over all proposals.
Metropolis抽样验收合格率。如果amcmc被使用，那么这将是一个矩阵的每个参数的录取率在每个批处理结束。否则，采样器是一个国际大都会，一个联合提案的所有参数，平均合格率为提出的所有提案。
参数：acceptance.w
If this is a non-predictive process model and amcmc is used then this will be a matrix of each random spatial effects acceptance rate at the end of each batch.
如果这是一个非预测的过程模型和amcmc被使用，那么这将是一个矩阵的每个随机空间效果接受率在每个批处理结束。
参数：acceptance.w.str
If this is a predictive process model and amcmc is used then this will be a matrix of each random spatial effects acceptance rate at the end of each batch.
如果这是一个预测的过程模型和amcmc使用，那么这将是一个矩阵的每个随机空间效果接受率在每个批处理结束。
参数：sp.effects
a matrix that holds samples from the posterior distribution of the spatial random effects. The rows of this matrix correspond to the n point observations and the columns are the posterior samples.
一个矩阵保存样本的空间随机效应的后验分布。这个矩阵的行对应的n点的观测和列后的样品。
The return object might include additional data used for subsequent prediction and/or model fit evaluation.
返回的对象可能包括额外的数据，用于以后的预测和/或模型拟合的评价。
（作者）----------Author(s)----------
Andrew O. Finley <a href="mailto:finleya@msu.edu">finleya@msu.edu</a>, <br>
Sudipto Banerjee <a href="mailto:baner009@umn.edu">baner009@umn.edu</a>
参考文献----------References----------
approach to quantifying uncertainty in multi-source forest area estimates. Environmental and Ecological Statistics, 15:241–258.
Predictive Process Models for Large Spatial Datasets. Journal of the Royal Statistical Society Series B, 70:825–848.
Hierarchical spatial modeling of additive and dominance genetic variance for large spatial trial datasets. Biometrics. DOI: 10.1111/j.1541-0420.2008.01115.x
large datasets. Computational Statistics and Data Analysis, DOI: 10.1016/j.csda.2008.09.008
参见----------See Also----------
spGGT,spMvLM, spMvGLM
spGGT，spMvLM，spMvGLM
实例----------Examples----------
## Not run: [＃不运行：]
library(MASS)
library(MBA)
library(fields)
################################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃]
##Spatial multivariate binomial[＃空间的多元二项式]
################################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃]
##Generate some data[＃生成一些数据。]
n <- 50
q <- 3
nltr <- q*(q-1)/2+q
coords <- cbind(runif(n,1,100),runif(n,1,100))
theta <- rep(3/50,q)
A <- matrix(0,q,q)
A[lower.tri(A,TRUE)] <- rnorm(nltr,1,1)
K <- A%*%t(A)
Psi <- diag(0,q)
c1 <- mvCovInvLogDet(coords=coords, knots=coords,
cov.model="exponential",
V=K, Psi=Psi, theta=theta,
modified.pp=TRUE, SWM=FALSE)
w <- mvrnorm(1,rep(0,nrow(c1$C)),c1$C)
X <- mkMvX(list(matrix(1,n,1), matrix(1,n,1), matrix(1,n,1)))
beta <- c(-1,0,1)
weight <- 10 ##i.e., trials [＃即试验]
p <- 1/(1+exp(-(X%*%beta+w)))
y <- rbinom(n*q, rep(weight,n*q), prob=p)
y.1 <- y[seq(1,length(y),q)]
y.2 <- y[seq(2,length(y),q)]
y.3 <- y[seq(3,length(y),q)]
##Specify starting values and collect samples[＃指定的初始值，采集样品]
A.starting <- diag(1,q)[lower.tri(diag(1,q), TRUE)]
fit <- glm((y/weight)~X-1, weights=rep(weight, n*q), family="binomial")
beta.starting <- coefficients(fit)
beta.tuning <- t(chol(vcov(fit)))
n.batch <- 200
batch.length <- 50
n.samples <- n.batch*batch.length
m.1 <- spMvGLM(list(y.1~1,y.2~1,y.3~1), weights=matrix(weight,n,q),
family="binomial", coords=coords,
starting=
list("beta"=beta.starting, "phi"=rep(0.06,q),"A"=A.starting, "w"=0),
tuning=
list("beta"=beta.tuning, "phi"=rep(0.05,q),"A"=rep(0.01,nltr), "w"=0.001),
priors=
list("beta.Flat", "phi.Unif"=rep(c(0.03, 0.3),q),"K.IW"=list(q+1, diag(0.1,q))),
amcmc=list("n.batch"=n.batch,"batch.length"=batch.length,"accept.rate"=0.43),
cov.model="exponential", sub.sample=c(5000,n.samples,10),
verbose=TRUE, n.report=100)
print(summary(mcmc(m.1$p.samples)))
beta.hat <- colMeans(m.1$p.samples[,1:q])
w.hat <- rowMeans(m.1$sp.effects)
p.hat <- 1/(1+exp(-(X%*%beta.hat+w.hat)))
p.hat.1 <- p.hat[seq(1,length(p.hat),q)]
p.hat.2 <- p.hat[seq(2,length(p.hat),q)]
p.hat.3 <- p.hat[seq(3,length(p.hat),q)]
##Take a look[＃看一看]
par(mfrow=c(3,2))
surf <- mba.surf(cbind(coords,y.1/weight),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image.plot(surf, main="Observed rates")
surf <- mba.surf(cbind(coords,p.hat.1),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image.plot(surf, main="Fitted rates")
surf <- mba.surf(cbind(coords,y.2/weight),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image.plot(surf)
surf <- mba.surf(cbind(coords,p.hat.2),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image.plot(surf)
surf <- mba.surf(cbind(coords,y.3/weight),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image.plot(surf)
surf <- mba.surf(cbind(coords,p.hat.3),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image.plot(surf)
##Prediction[＃预测]
m <- 200
pred.coords <- cbind(runif(m,1,100),runif(m,1,100))
pred.covars <- mkMvX(list(matrix(1,m,1), matrix(1,m,1), matrix(1,m,1)))
pred <- spPredict(m.1, pred.coords, pred.covars)
pred.p <- rowMeans(pred$y.pred)
pred.p.1 <- pred.p[seq(1,length(pred.p),q)]
pred.p.2 <- pred.p[seq(2,length(pred.p),q)]
pred.p.3 <- pred.p[seq(3,length(pred.p),q)]
##Take a look[＃看一看]
par(mfrow=c(3,2))
surf <- mba.surf(cbind(coords,y.1/weight),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image.plot(surf, main="Observed rates")
surf <- mba.surf(cbind(pred.coords,pred.p.1),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image.plot(surf, main="Fitted rates")
surf <- mba.surf(cbind(coords,y.2/weight),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image.plot(surf)
surf <- mba.surf(cbind(pred.coords,pred.p.2),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image.plot(surf)
surf <- mba.surf(cbind(coords,y.3/weight),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image.plot(surf)
surf <- mba.surf(cbind(pred.coords,pred.p.3),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image.plot(surf)
################################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃]
##Spatial multivariate poisson[＃空间的多元Poisson]
################################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃]
##Generate some data[＃生成一些数据。]
n <- 100
q <- 3
nltr <- q*(q-1)/2+q
coords <- cbind(runif(n,1,100),runif(n,1,100))
theta <- rep(3/50,q)
A <- matrix(0,q,q)
A[lower.tri(A,TRUE)] <- rnorm(nltr,1,1)
K <- A%*%t(A)
Psi <- diag(0,q)
c1 <- mvCovInvLogDet(coords=coords, knots=coords,
cov.model="exponential",
V=K, Psi=Psi, theta=theta,
modified.pp=TRUE, SWM=FALSE)
w <- mvrnorm(1,rep(0,nrow(c1$C)),c1$C)
X <- mkMvX(list(matrix(1,n,1), matrix(1,n,1), matrix(1,n,1)))
beta <- c(-1,0,1)
y <- rpois(n*q, exp(X%*%beta+w))
y.1 <- y[seq(1,length(y),q)]
y.2 <- y[seq(2,length(y),q)]
y.3 <- y[seq(3,length(y),q)]
##Specify starting values and collect samples. For[＃指定的初始值，并收集样本。为]
##a true analysis, several longer chains should be[＃一个真实的分析，应该是几个长链]
##run.[＃运行。]
A.starting <- diag(1,q)[lower.tri(diag(1,q), TRUE)]
beta.starting <- coefficients(glm(y~X-1, family="poisson"))
beta.tuning <- t(chol(vcov(glm(y~X-1, family="poisson"))))
n.samples <- 5000
m.1 <- spMvGLM(list(y.1~1,y.2~1,y.3~1), family="poisson", coords=coords,
starting=
list("beta"=beta.starting, "phi"=rep(0.06,q),
"A"=A.starting, "w"=0),
tuning=
list("beta"=beta.tuning, "phi"=rep(0.01,q),
"A"=rep(0.005,nltr), "w"=0.001),
priors=
list("beta.Flat", "phi.Unif"=rep(c(0.03, 0.3),q),
"K.IW"=list(q+1, diag(0.1,q))),
cov.model="exponential",
n.samples=n.samples, sub.sample=c(5000,n.samples,10),
verbose=TRUE, n.report=500)
m.1$p.samples[,paste("phi_",1:q,sep="")] <-
3/m.1$p.samples[,paste("phi_",1:q,sep="")]
print(summary(mcmc(m.1$p.samples)))
beta.hat <- colMeans(m.1$p.samples[,1:q])
w.hat <- rowMeans(m.1$sp.effects)
y.hat <- exp(X%*%beta.hat+w.hat)
y.hat.1 <- y.hat[seq(1,length(y.hat),q)]
y.hat.2 <- y.hat[seq(2,length(y.hat),q)]
y.hat.3 <- y.hat[seq(3,length(y.hat),q)]
##Take a look[＃看一看]
par(mfrow=c(3,2))
surf <- mba.surf(cbind(coords,y.1),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image(surf, main="Observed counts")
contour(surf, drawlabels=FALSE, add=TRUE)
text(coords, labels=y.1, cex=1, col="blue")
surf <- mba.surf(cbind(coords,y.hat.1),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image(surf, main="Fitted counts")
contour(surf, add=TRUE)
contour(surf, drawlabels=FALSE, add=TRUE)
text(coords, labels=round(y.hat.1,0), cex=1, col="blue")
surf <- mba.surf(cbind(coords,y.2),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image(surf)
contour(surf, drawlabels=FALSE, add=TRUE)
text(coords, labels=y.2, cex=1, col="blue")
surf <- mba.surf(cbind(coords,y.hat.2),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image(surf)
contour(surf, drawlabels=FALSE, add=TRUE)
text(coords, labels=round(y.hat.2,0), cex=1, col="blue")
surf <- mba.surf(cbind(coords,y.3),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image(surf)
contour(surf, drawlabels=FALSE, add=TRUE)
text(coords, labels=y.3, cex=1, col="blue")
surf <- mba.surf(cbind(coords,y.hat.3),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image(surf)
contour(surf, drawlabels=FALSE, add=TRUE)
text(coords, labels=round(y.hat.3,0), cex=1, col="blue")
## End(Not run)[＃（不执行）]
转载请注明:出自 生物统计家园网(http://www.biostatistic.net)。
注：
注1：为了方便大家学习，本文档为生物统计家园网机器人LoveR翻译而成，仅供个人R语言学习参考使用，生物统计家园保留版权。
注2：由于是机器人自动翻译，难免有不准确之处，使用时仔细对照中、英文内容进行反复理解，可以帮助R语言的学习。
注3：如遇到不准确之处，请在本贴的后面进行回帖，我们会逐渐进行修订。
|
发表于 2012-9-30 14:24:30
