spGGT(spBayes)
spGGT()所属R语言包：spBayes
Function for fitting univariate and multivariate
功能配件单因素和多
译者：生物统计家园网 机器人LoveR
描述----------Description----------
The function spGGT fits Gaussian univariate and multivariate Bayesian spatial regression models. In most cases, we encourage users to use spLM or spLM over spGGT. These functions offer a simplified interface and additional functionality to handle large data sets.
的功能spGGT符合高斯单因素和多因素贝叶斯空间回归模型。在大多数情况下，我们鼓励用户使用spLM或spLM经spGGT。这些功能提供了一个简化的界面和附加功能处理大型数据集。
用法----------Usage----------
spGGT(formula, data = parent.frame(), coords, run.control,
var.update.control, beta.update.control, cov.model,
verbose=TRUE, ...)
参数----------Arguments----------
参数：formula
for a univariate model, this is a symbolic description of the regression model to be fit. For a multivariate model, this is a list of univariate formulas. See example below.
单变量模型中，这是一个象征性的回归模型是适合的描述。对于多变量模型，这是一个单变量的公式列表。见下面的例子。
参数：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 spGGT is called.
一个可选的数据框包含在模型中的变量。如果没有找到，数据，变量environment(formula)，通常是spGGT被称为环境。
参数：coords
an n x 2 matrix of the observation coordinates in R^2 (e.g., easting and northing).
n x 2R^2（例如，东向和北向）的观察坐标矩阵。
参数：run.control
a list with tags: n.samples the value of which is the number of MCMC iterations; sp.effects a logical value indicating if spatial random effects should be recovered; DIC a logical value indicating if marginalized DIC and associated statistics should be computed; DIC.start and DIC.end the MCMC iteration interval over which DIC should be calculated; and linpack a logical value indicating if LINPACK or LAPACK should be used for matrix algebra operations. The default is to use LINPACK.
一个列表标签：n.samples其价值是多少MCMC迭代; sp.effects一个逻辑值，表示如果空间随机效应应该被回收;DIC一个逻辑值，表示如果被边缘化DIC和相关的统计数据应计算; DIC.start和DIC.endMCMC迭代间隔，其中DIC应计算;和linpack一个逻辑值，表示如果LINPACK或LAPACK应使用矩阵代数运算。默认情况下是使用LINPACK。
参数：var.update.control
a list with each tag corresponding to a parameter name used to define the variance structure of the model. Valid list tags are K, Psi, phi, and nu. The value portion of each of these tags is another list with tags sample.order, starting, tuning, and prior. The positive scalar value of the optional tag sample.order determines the order in which parameters are updated by sequential Metropolis-Hastings steps. The sample.order value is relative to the sample order specified for the other parameters. The value of the starting tag defines the starting value of the parameter in the Metropolis-Hastings sampling. If the list element defines a vector or matrix of parameters, then the starting value should be set accordingly. If a scalar starting value is given for a vector or matrix of parameters then the scalar value is recycled. If the model defines K or Psi as a full cross-covariance matrix, opposed to a diagonal matrix or scalar, then the starting value should be the lower-triangle of the Cholesky square root of the desired starting matrix value. If K or Psi is a diagonal matrix, then the starting values should be passed as a vector or scalar. The value associated with the tuning tag defines the step size of the proposal used in the Metropolis-Hastings sampling. Again the tuning value is recycled if needed. If the parameter element is a matrix, i.e., K or Psi, then the diagonal elements in the tuning matrix correspond to the column major lower-triangle elements in the parameter matrix. The value of the prior tag is a valid prior object which describes the prior distribution of the parameter. Any parameter in the specified model can be fixed by using the FIXED prior with the starting value set to the desired fixed value. See example below.
与每个标记对应于用来定义的方差结构的模型的参数名称的列表。有效的列表标签是K，Psi，phi，nu。这些标签的价值部分是另一个标签列表sample.order，starting，tuning和prior。可选的标记sample.order的阳性标量值确定通过顺序大都市-赫斯廷斯步骤的参数进行更新的顺序进行。 sample.order这个值是相对于其他参数指定的样品订单。 starting标签的值定义在大都市黑斯廷斯采样参数的初始值。如果列表中的元素定义了一个向量或矩阵参数，那么应相应地设置初始值。如果一个标量初始值给出的矢量或矩阵的参数，然后在标量值被回收。如果模型定义K或Psi作为一个完整的互协方差矩阵，而不是一个对角矩阵或标量，那么初始值应该是三角的Cholesky平方根所需的开始矩阵值。如果K或Psi是一个对角矩阵，则初始值应该传递一个向量或标量。与tuning标签的值定义的步长在大都市黑斯廷斯采样的建议。同样的整定值，如果需要回收。如果该参数元素是一个矩阵，即，K或Psi，然后调谐矩阵中的对角线上的元素对应的列主要下部三角形参数矩阵中的元素。 prior标签的值是一个有效的prior描述对象的参数的先验分布。任何参数，可以固定在指定的模型通过使用FIXEDpriorstarting这个值设置到理想的固定值。见下面的例子。
参数：beta.update.control
a list with tags update, tuning, starting, prior, and sample.order. The value of update is either "GIBBS" or "MH". With the "MH" option, the value of the optional tag tuning is the upper-triangle of the Cholesky decomposed p x p tuning matrix, where p is the number of regression parameters. The tuning matrix controls the proposal step size in the Metropolis-Hastings sampling. The tuning values along the diagonal of this matrix correspond to the order of the regressors in the model formula. Off-diagonal elements are the covariance among the terms. The optional tag starting receives a vector of length p, whose elements define the starting values of the regressors used in the MCMC sampling. The order of the starting values in this vector corresponds to the order of the regressors in the model formula. The optional tag prior is a valid prior object that defines either a normal or flat distribution. If prior is not included in the beta.update.control list then the prior distribution on the regressors is assumed flat (i.e., proportional to 1.0). With the "MH" option, the value of the optional tag sample.order determines the order in which parameters are updated by sequential Metropolis-Hastings steps, see below for details.
标签的列表update，tuning，starting，prior和sample.order。 update是"GIBBS"或"MH"。 "MH"选项，可选的标记值tuning是上三角的Cholesky分解p x p调整矩阵，其中p是回归参数的数量。调整矩阵控制在大都市黑斯廷斯采样步长的建议。沿该矩阵的对角线的调谐值对应的回归量的模型公式中的顺序。非对角线元素是术语之间的协方差。可选的标记starting接收一个矢量的长度的p，其元素定义初始值的回归系数在MCMC采样。此向量中的初始值的顺序对应的顺序在模型公式中的回归量。可选的标记prior是有效的prior对象，定义正常或平面的分布。如果prior不包括在beta.update.control列表，然后回归量假定先验分布平坦的（即，比例为1.0）。随着"MH"选项，可选的标记sample.order的价值确定参数的更新顺序大都市黑斯廷斯步骤的顺序，请参阅下面的详细信息。
参数：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"。有关详细信息，请参见下文。
参数：verbose
if TRUE, model specification and progress of the sampler is printed to the screen. Otherwise, nothing is printed to the screen.
如果TRUE，型号规格和进步的取样打印到屏幕上。否则，一切都被打印到屏幕上。
参数：...
currently no additional arguments.
目前没有任何额外的参数。
Details
详细信息----------Details----------
Please refer to Section 3.1 in the vignette.
请参阅第3.1节中的小插曲。
值----------Value----------
An object of class spGGT, which is a list with the following tags:
类的一个对象spGGT，这是一个与下面的标签列表：
参数：cov.model
the covariance model specified by cov.model.
指定的cov.model的协方差模型。
参数：no.Psi
a logical value indicating if the parameter Psi was used.
一逻辑值，如果参数Psi使用。
参数：coords
the n x 2 matrix specified by coords.
n x 2的矩阵指定coords。
参数：X
the n x p matrix of regressors specified by formula.
n x p的回归系数矩阵指定的formula。
参数：Y
the the n x 1 response variable vector specified by formula.
在n x 1响应变量向量指定formula。
参数：p.samples
a matrix of posterior samples for the defined parameters. If K or Psi is specified as a full cross-covariance matrix then the corresponding samples in p.samples are column major lower-triangle elements of the matrix.
后的样本进行定义的参数的矩阵。如果K或Psi被指定为一个完整的互协方差矩阵，然后相应的样品中p.samples是按列的下三角矩阵的元素。
参数：acceptance
a matrix of the Metropolis-Hastings sampling acceptance rate. Row names correspond to the parameters updated with Metropolis-Hastings.
矩阵的大都市黑斯廷斯抽样验收率。行名对应的参数更新与大都市黑斯廷斯。
参数：DIC
a matrix that holds marginalized DIC and associated values is returned if DIC is true in the run.control list.
如果DIC是真实的run.control列表，一个矩阵，它拥有边缘化的DIC及相关值，则返回。
参数：sp.effects
a matrix that holds samples from the posterior distribution of the response specific spatial random effects. The rows of this matrix correspond to the n point observations and the columns are the posterior samples. If a multivariate model is fit to m response variables, the sp.effects matrix has mn rows. The random effects samples for points are held in rows 1:m, (m+1):2m, ..., ((i-1)m+1):im, ..., ((n-1)m+1):nm, where i = 1 ... n (e.g., the samples for the first point are in rows 1:m, second point in rows (m+1):2m, etc.).
一个矩阵保存样本的后验分布的响应特定的空间随机效应。这个矩阵的行对应的n点的观测和列后的样品。如果适合m的响应变量的多变量模型，sp.effects矩阵mn行。行1：米，第（m +1）：约2m，...，（第（i-1）m +1个）：肌肉注射，...，（第（n-1）米的随机效应样本点被保持在+1）：纳米，其中i = 1，... N（例如，为第一点的样品中的行1：米，行第（m +1）中的第二点：2米等）。
（作者）----------Author(s)----------
Andrew O. Finley <a href="mailto:finleya@msu.edu">finleya@msu.edu</a>, <br>
Sudipto Banerjee <a href="mailto:sudiptob@biostat.umn.edu">sudiptob@biostat.umn.edu</a>, <br>
Bradley P. Carlin <a href="mailto:brad@biostat.umn.edu">brad@biostat.umn.edu</a>.
参考文献----------References----------
参见----------See Also----------
prior, spPredict
prior，spPredict
实例----------Examples----------
## Not run: [＃不运行：]
data(FBC07.dat)
Y.1 <- FBC07.dat[1:100,"Y.1"]
Y.2 <- FBC07.dat[1:100,"Y.2"]
coords <- as.matrix(FBC07.dat[1:100,c("coord.X", "coord.Y")])
#############################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃]
## Univariate models[＃单变量模型]
#############################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃]
##Fit some competing models.[＃适合一些竞争机型。]
##Full model with nugget (Psi), partial sill (K),[＃模型金块（PSI），偏基台（K），]
##and spatial decay (Phi).[＃和空间衰减（PHI）。]
K.prior <- prior(dist="IG", shape=2, scale=5)
Psi.prior <- prior(dist="IG", shape=2, scale=5)
phi.prior <- prior(dist="UNIF", a=0.06, b=3)
var.update.control <-
list("K"=list(starting=5, tuning=0.5, prior=K.prior),
"Psi"=list(starting=5, tuning=0.5, prior=Psi.prior),
"phi"=list(starting=0.1, tuning=0.5, prior=phi.prior)
)
beta.control <- list(update="GIBBS", prior=prior(dist="FLAT"))
run.control <- list("n.samples"=1000, "sp.effects"=TRUE)
Fit.m1 <- spGGT(formula=Y.2~1, run.control=run.control,
coords=coords, var.update.control=var.update.control,
beta.update.control=beta.control,
cov.model="exponential")
plot(Fit.m1$p.samples)
summary(Fit.m1$p.samples)
##First reduced model with no nugget just partial sill and decay.[＃首先，没有金块只是局部的窗台和衰减模式减少。]
var.update.control <-
list("K"=list(starting=5, tuning=0.5, prior=K.prior),
"phi"=list(starting=0.1, tuning=0.5, prior=phi.prior)
)
Fit.m2 <- spGGT(formula=Y.2~1, run.control=run.control,
coords=coords, var.update.control=var.update.control,
beta.update.control=beta.control,
cov.model="exponential")
plot(Fit.m2$p.samples)
##Second reduced model with no nugget just partial sill and fixed[＃二没有金块的部分窗台和固定的简化模型]
##decay at the WRONG value, which forces the wrong estimate of K.[＃在错误的值衰减，这迫使错误的估计K.]
phi.prior <- prior(dist="FIXED")
var.update.control <-
list("K"=list(starting=5, tuning=0.5, prior=K.prior),
"phi"=list(starting=0.1, prior=phi.prior)
)
Fit.m3 <- spGGT(formula=Y.2~1, run.control=run.control,
coords=coords, var.update.control=var.update.control,
beta.update.control=beta.control,
cov.model="exponential")
plot(Fit.m3$p.samples)
##For this data, DIC cannot distinguish between the first two[＃对于这个数据，DIC无法区分第一]
##models but does show that the third model is clearly[＃模式，但确实表明，第三个模型是清楚]
##wrong. An empirical semivariogram is useful[＃错误的。经验半变异函数是很有用的]
##for recognizing that the full model is the correct choice[＃认识到完整的模型是正确的选择]
##over the second model.[＃在第二个模型。]
print(spDiag(Fit.m1))
print(spDiag(Fit.m2))
print(spDiag(Fit.m3))
##Use akima or MBA to produce interpolated surfaces of[＃使用Akima或MBA插值表面]
##estimated random spatial effects.[＃估计随机的空间效果。]
m1.surf <- mba.surf(cbind(coords, rowMeans(Fit.m1$sp.effects)),
no.X=100, no.Y=100, extend=TRUE)$xyz.est
image(m1.surf, xaxs="r", yaxs="r",
main="Y.2 random spatial effects")
points(coords, pch=19)
##Now the full model using the Matern spatial correlation function.[现在全模型，利用Matern的空间相关性函数。]
##Also, for fun, update the spatial decay parameters[＃此外，为了好玩，更新的空间衰减参数]
##phi and nu in a separate MH block. This provides finer control[＃phi和怒江在一个单独的MH块。这提供了更精细的控制]
##on the parameters' acceptance rate, but takes a bit longer to[＃参数的录取率，但需要花一点时间来]
##run.[＃运行。]
phi.prior <- prior(dist="UNIF", a=0.06, b=3)
nu.prior <- prior(dist="UNIF", a=0.01, b=2)
var.update.control <-
list("K"=list(starting=5, tuning=0.5, prior=K.prior),
"Psi"=list(starting=5, tuning=0.5, prior=Psi.prior),
"phi"=list(sample.order=1, starting=0.1,
tuning=0.5, prior=phi.prior),
"nu"=list(sample.order=1, starting=0.1,
tuning=0.5, prior=nu.prior)
)
run.control <- list("n.samples"=1000, "sp.effects"=FALSE)
Fit.m4 <-
spGGT(formula=Y.2~1, run.control=run.control,
coords=coords, var.update.control=var.update.control,
beta.update.control=beta.control,
cov.model="matern")
plot(Fit.m4$p.samples)
##Solve the Matern correlation function for the estimated[】＃解答的Matern的的相关函数的估计]
##95% CI for effective decay (e.g., distance at which the[＃95％CI为有效衰变（如距离]
##correlation drops to 0.05).[＃相关性下降到0.05）。]
phi.hat <- quantile(Fit.m4$p.samples[,"Phi"], c(0.50,0.025,0.975))
nu.hat <- quantile(Fit.m4$p.samples[,"Nu"], c(0.50,0.025,0.975))
matern <- function(cor, eff.d, phi, nu){
cor - (eff.d*phi)^nu/(2^(nu-1)*gamma(nu))*besselK(x=eff.d*phi, nu=nu)
}
##50[＃50]
print(uniroot(matern, interval=c(1, 50),
cor=0.05, phi=phi.hat[1], nu=nu.hat[1])$root)
print(uniroot(matern, interval=c(1, 50),
cor=0.05, phi=phi.hat[2], nu=nu.hat[2])$root)
print(uniroot(matern, interval=c(1, 50),
cor=0.05, phi=phi.hat[3], nu=nu.hat[3])$root)
##Finally, the full model again but this time update the Beta[＃最后，完整的模型，但是这一次更新的Beta]
##using MH. Using the default of sample.order=0, this sampling[＃MH。使用默认的sample.order = 0，本次抽查中]
##scheme is simply a single block MH proposal of all model parameters.[＃方案是一个简单的单块MH建议所有的模型参数。]
##Also, the first run uses the default tuning value for the Beta[＃另外，第一次运行时将使用默认的整定值的测试版]
##(i.e., chol(vcov(lm(Y~X))). Both chains use a normal prior on Beta;[＃（即哲（vcov（LM（Y~X）））。两家连锁店使用前正常测试;]
##however, the first chain allows for a vague prior, where as the[＃然而，第一链允许前一个模糊的，其中的]
##second is much too strict.[＃二是过于严格。]
var.update.control <-
list("K"=list(starting=5, tuning=0.5, prior=K.prior),
"Psi"=list(starting=5, tuning=0.5, prior=Psi.prior),
"phi"=list(starting=0.1, tuning=0.5, prior=phi.prior)
)
##Chain 1: vague prior normal variance of 1/0.001 = 1000[链1：1/0.001 = 1000含糊之前的正常变异]
beta.prior <- prior(dist="NORMAL", mu=5, precision=0.001)
beta.control <- list(update="MH", prior=beta.prior)
Fit.m5.a <-
spGGT(formula=Y.2~1, run.control=run.control,
coords=coords, var.update.control=var.update.control,
beta.update.control=beta.control,
cov.model="exponential")
##Chain 2: restrictive prior normal variance of 1/5 = 0.2[＃链2：限制性事先正常方差的1/5 = 0.2]
##(i.e., crazy strict)[＃（即，疯狂严格）]
beta.prior <- prior(dist="NORMAL", mu=5, precision=5)
beta.control <- list(update="MH", prior=beta.prior, tuning=0.75)
Fit.m5.b <-
spGGT(formula=Y.2~1, run.control=run.control,
coords=coords, var.update.control=var.update.control,
beta.update.control=beta.control,
cov.model="exponential")
plot(mcmc.list(Fit.m5.a$p.samples, Fit.m5.b$p.samples))
##Now just for fun, here is the full model but using the weakly[现在只是为了好玩，这里是完整的模型，但使用弱]
##informative half-Cauchy for the nugget (Psi) and partial sill (K).[＃翔实半柯西的金块（PSI）和部分的窗台（K）。]
K.prior <- prior(dist="HC", a=50)
Psi.prior <- prior(dist="HC", a=50)
phi.prior <- prior(dist="UNIF", a=0.06, b=3)
var.update.control <-
list("K"=list(starting=5, tuning=0.5, prior=K.prior),
"Psi"=list(starting=5, sample.order=1, tuning=2, prior=Psi.prior),
"phi"=list(starting=0.1, sample.order=2, tuning=1, prior=phi.prior)
)
beta.control <- list(update="GIBBS", prior=prior(dist="FLAT"))
run.control <- list("n.samples"=3000, "sp.effects"=TRUE)
Fit.m1 <- spGGT(formula=Y.2~1, run.control=run.control,
coords=coords, var.update.control=var.update.control,
beta.update.control=beta.control,
cov.model="exponential")
plot(Fit.m1$p.samples)
summary(Fit.m1$p.samples)
#############################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃]
## Multivariate models[＃多变量模型]
#############################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃]
##Full model with non-spatial cross-covariance[与非空间的互协方差模型＃]
##(Psi) matrix, spatial cross-covariance matrix (K),[＃（帕普西）矩阵，空间互协方差矩阵（K），]
##and response specific spatial decay parameter (Phi)[和响应特定的空间衰减参数（披）]
##(i.e., a non-separable model).[＃（即，非模型可分离）。]
##These chains really need to be run for several thousand[＃这些连锁店确实需要运行几千年]
##samples (e.g., 10,000). Also note that it is much easier to maintain[＃样品（例如，10000）。另外请注意，这是非常容易维护]
##a healthy acceptance rate for each parameter if they[＃一个健康的录取率为每一个参数，如果他们]
##are updated individually using the sample.order[＃单独更新使用的sample.order的]
##directive, as shown in the second example below.[＃指令，如下面的第二个例子中所示。]
K.prior <- prior(dist="IWISH", df=2, S=diag(c(3, 6)))
Psi.prior <- prior(dist="IWISH", df=2, S=diag(c(7, 5)))
phi.prior <- prior(dist="UNIF", a=0.06, b=3)
##The matrix starting values for K and Psi are actually[＃K和PSI的初始值的矩阵实际上是]
##passed as the sqrt of the desired starting values.[＃通过的sqrt所需的初始值。]
##This is only done if K or Psi are full cross-covariance[＃这是唯一的，如果K或PSI是互协方差]
##matrices and serve to insure that the true starting value[＃矩阵和服务，以确保真正的起始值]
##matrix is positive definite.[＃矩阵是正定的。]
K.starting <- matrix(c(2,-3, 0, 1), 2, 2)
Psi.starting <- diag(c(3, 2))
var.update.control <-
list("K"=list(starting=K.starting, tuning=diag(c(0.08, 0.3, 0.08)),
prior=K.prior),
"Psi"=list(starting=Psi.starting, tuning=diag(c(0.08, 0.3, 0.08)),
prior=Psi.prior),
"phi"=list(starting=0.1, tuning=0.5,
prior=list(phi.prior, phi.prior))
)
beta.control <- list(update="GIBBS", prior=prior(dist="FLAT"))
run.control <- list("n.samples"=2000, "sp.effects"=FALSE)
Fit.mv.m1 <-
spGGT(formula=list(Y.1~1, Y.2~1), run.control=run.control,
coords=coords, var.update.control=var.update.control,
beta.update.control=beta.control,
cov.model="exponential")
plot(Fit.mv.m1$p.samples)
##Calculate the mean K correlation matrix. Note that since K is[＃计算的K均值相关系数矩阵。请注意，由于K是]
##specified as a full cross-covariance matrix, the posterior[＃指定作为一个完整的互协方差矩阵，后]
##samples are returned for the lower triangle of the symmetric K,[＃样品返回的下三角对称的k，]
##in column major order.[＃列中的大订单。]
K <- matrix(0, 2, 2)
K[lower.tri(K, diag=TRUE)] <-
colMeans(Fit.mv.m1$p.samples[,c("K_1","K_2","K_3")])
K[upper.tri(K, diag=FALSE)] <- t(K)[upper.tri(K, diag=FALSE)]
print(cov2cor(K))
##Now using sample order. Note this takes ~3 times as long[＃现在使用的样品订单。请注意，这需要~3倍长]
##to run, but provides much finer control of acceptance rates.[＃运行，但录取率提供了更精细的控制。]
var.update.control <-
list("K"=list(sample.order=0, starting=K.starting,
tuning=diag(c(0.3, 0.5, 0.2)), prior=K.prior),
"Psi"=list(sample.order=1, starting=Psi.starting,
tuning=diag(c(0.2, 0.5, 0.2)), prior=Psi.prior),
"phi"=list(sample.order=2, starting=0.1,
tuning=0.5, prior=list(phi.prior, phi.prior))
)
run.control <- list("n.samples"=2000, "sp.effects"=TRUE)
Fit.mv.m2 <-
spGGT(formula=list(Y.1~1, Y.2~1), run.control=run.control,
coords=coords, var.update.control=var.update.control,
beta.update.control=beta.control,
cov.model="exponential")
plot(Fit.mv.m2$p.samples)
##Now again check out the random spatial effects. Recall,[＃现在又检查出随机的空间效果。回想一下，]
##these are stored in the sp.effects matrix, which is organized[＃这些存储在的sp.effects矩阵，这是组织]
##as m consecutive rows for each point.[＃m个连续行的每个点。]
rand.effect.Y.1 <-
Fit.mv.m2$sp.effects[seq(1, nrow(Fit.mv.m2$sp.effects), 2),]
rand.effect.Y.2 <-
Fit.mv.m2$sp.effects[seq(2, nrow(Fit.mv.m2$sp.effects), 2),]
rand.effect.Y.1.surf <-
mba.surf(cbind(coords, rowMeans(rand.effect.Y.1)),
no.X=100, no.Y=100, extend=TRUE)$xyz.est
rand.effect.Y.2.surf <-
mba.surf(cbind(coords, rowMeans(rand.effect.Y.2)),
no.X=100, no.Y=100, extend=TRUE)$xyz.est
##The estimated negative value of off-diagonal element in K is[＃估计负值的非对角线上元素在K]
##clearly seen in the map of random spatial effects.[＃清楚地看到在图中的随机的空间效果。]
par(mfrow=c(1,2))
image(rand.effect.Y.1.surf, xaxs="r", yaxs="r",
main="Y.1 random spatial effects")
contour(rand.effect.Y.1.surf, add=TRUE)
image(rand.effect.Y.2.surf, xaxs="r", yaxs="r",
main="Y.2 random spatial effects")
contour(rand.effect.Y.2.surf, add=TRUE)
##Now a simpler model. Specifically, the full model shows that[＃一个简单的模型。具体而言，完整的模型显示，]
##there is negligible non-spatial cross-covariance, therefore,[＃是可以忽略不计的非空间互协方差，因此，]
##Psi might just be modeled with a diagonal covariance matrix.[＃PSI可能只是一个对角协方差矩阵进行建模。]
##This is actually the model which gave rise to this synthetic[＃这实际上是产生了这种合成的模式，]
##data set.[＃数据集。]
K.prior <- prior(dist="IWISH", df=2, S=diag(c(3, 6)))
Psi.prior <- prior(dist="IG", shape=2, scale=5)
phi.prior <- prior(dist="UNIF", a=0.06, b=3)
var.update.control <-
list("K"=list(sample.order=0, starting=K.starting,
tuning=diag(c(0.3, 0.5, 0.2)), prior=K.prior),
"Psi"=list(sample.order=1, starting=5,
tuning=0.5, prior=list(Psi.prior, Psi.prior)),
"phi"=list(sample.order=2, starting=0.1,
tuning=0.5, prior=list(phi.prior, phi.prior))
)
run.control <- list("n.samples"=2000)
Fit.mv.m3 <-
spGGT(formula=list(Y.1~1, Y.2~1), run.control=run.control,
coords=coords, var.update.control=var.update.control,
beta.update.control=beta.control,
cov.model="exponential")
plot(Fit.mv.m3$p.samples)
summary(Fit.mv.m3$p.samples)
##Finally, a really simple model which we know makes no sense[＃最后，一个非常简单的模型，我们知道这是没有意义的]
##for this data set. Specifically, a single variance parameter[＃这组数据。具体而言，一个单一的方差参数]
##for the diagonal elements of K and Psi and common spatial decay,[编号K和Psi和共同空间衰减的对角线元素，]
##phi. Also, we use a spherical spatial covariance function, just[＃PHI。此外，我们使用了球形空间协方差函数，只]
##because we can.[＃因为我们能做到。]
K.prior <- prior(dist="IG", shape=2, scale=5)
Psi.prior <- prior(dist="IG", shape=2, scale=5)
phi.prior <- prior(dist="UNIF", a=0.06, b=3)
var.update.control <-
list("K"=list(starting=5, tuning=0.5, prior=K.prior),
"Psi"=list(starting=5, tuning=0.5, prior=Psi.prior),
"phi"=list(starting=0.1, tuning=0.5, prior=phi.prior)
)
Fit.mv.m4 <-
spGGT(formula=list(Y.1~1, Y.2~1), run.control=run.control,
coords=coords, var.update.control=var.update.control,
beta.update.control=beta.control,
cov.model="spherical")
plot(Fit.mv.m4$p.samples)
## End(Not run)[＃（不执行）]
转载请注明:出自 生物统计家园网(http://www.biostatistic.net)。
注：
注1：为了方便大家学习，本文档为生物统计家园网机器人LoveR翻译而成，仅供个人R语言学习参考使用，生物统计家园保留版权。
注2：由于是机器人自动翻译，难免有不准确之处，使用时仔细对照中、英文内容进行反复理解，可以帮助R语言的学习。
注3：如遇到不准确之处，请在本贴的后面进行回帖，我们会逐渐进行修订。
|
发表于 2012-9-30 14:24:07
