Description Usage Arguments Details Value References Examples
func.cle.ord
performs composite likelihood estimation of parameters and their standard errors in a spatial ordered probit model by maximizing its composite log-likelihood.
1 2 | func.cle.ord(vec.yobs, mat.X, mat.lattice, radius, n.cat, n.sim = 100,
parallel = TRUE, n.core = max(detectCores()/2, 1), output = TRUE)
|
vec.yobs |
a vector of observed responses for all N sites. |
mat.X |
regression (design) matrix, including intercepts. |
mat.lattice |
a data matrix containing geographical information of sites. The ith row constitutes a set of geographical coordinates. |
radius |
weight radius. |
n.cat |
number of categories. |
n.sim |
number of simulations used for calculate the Godambe matrix (default: 100). |
parallel |
logical flag indicates using parallel processing (default: |
n.core |
number of physical cores used for parallel processing (when |
output |
logical flag indicates whether printing out result (default: |
Given the design matrix, the vector of observed responses, spatial lattice data, weight radius, number of categories, and the prespecified number of simulated vectors of responses used in estimating the Godambe information, this function assumes initial values of cutoff points and β as the estimates from the standard ordered probit regression with independent responses. After initial reparameterization, it first estimates parameters of interest by maximizing the composite log-likelihood using optim
, then computes the reparameterized sample covariance matrix and the set of standard errors, and finally reverse the reparameterization to obtain estimates corresponding to the original parameterization.
func.cle.ord
returns a list containing:
vec.par
: a vector of estimator for θ=(cutoff,β,σ^2,ρ);
vec.se
: a vector of standard error for the estimator;
mat.asyvar
: estimated asymptotic covariance matrix H^{-1}(θ)J(θ)H^{-1}(θ) for the estimator; and
vec.comp
: a vector of computational time for parameter and standard error estimation.
CLIC
: Composite likelihood information criterion proposed by Varin and Vidoni (2005), i.e. -2*logCL(θ) + 2*trace(H^{-1}(θ)J(θ))
Feng, Xiaoping, Zhu, Jun, Lin, Pei-Sheng, and Steen-Adams, Michelle M. (2014) Composite likelihood Estimation for Models of Spatial Ordinal Data and Spatial Proportional Data with Zero/One values. Environmetrics 25(8): 571–583.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 | # Example of n.cat = 3 (Spatial ordinal regression)
# True parameter
vec.cutoff <- 2; vec.beta <- c(1, 2, 1, 0, -1); sigmasq <- 0.8; rho <- 0.6; radius <- 5
vec.par <- c(vec.cutoff, vec.beta, sigmasq, rho)
# Coordinate matrix
n.cat <- 3; n.lati <- 30; n.long <- 30
n.site <- n.lati * n.long
mat.lattice <- cbind(rep(1:n.lati, n.long), rep(1:n.long, each=n.lati))
mat.dist <- as.matrix(dist(mat.lattice, upper=TRUE, diag=TRUE))
mat.cov <- sigmasq * rho^mat.dist
set.seed(1228)
# Generate regression (design) matrix with intercept
mat.X <- cbind(rep(1, n.site),scale(matrix(rnorm(n.site*(length(vec.beta)-1)),nrow=n.site)))
vec.Z <- t(chol(mat.cov)) %*% rnorm(n.site) + mat.X %*% vec.beta
vec.epsilon <- diag(sqrt(1-sigmasq), n.site) %*% rnorm(n.site)
vec.ylat <- as.numeric(vec.Z + vec.epsilon)
# Convert to the vector of observation
vec.yobs <- func.obs.ord(vec.ylat, vec.alpha=c(-Inf,0,vec.cutoff,Inf))
# With parallel computing
## Not run:
ord.example <- func.cle.ord(vec.yobs, mat.X, mat.lattice, radius, n.cat,
n.sim=100, parallel = TRUE, n.core = 2)
round(ord.example$vec.par,4)
# alpha1 beta0 beta1 beta2 beta3 beta4 sigma^2 rho
# 1.8395 0.9550 1.9690 0.9565 0.0349 -1.0398 0.8200 0.5578
round(ord.example$vec.se,4)
# alpha1 beta0 beta1 beta2 beta3 beta4 sigma^2 rho
# 0.1602 0.1222 0.1463 0.0916 0.0485 0.0889 0.1935 0.1267
## End(Not run)
# Without parallel computing
## Not run:
ord.example2 <- func.cle.ord(vec.yobs, mat.X, mat.lattice, radius,
n.cat, n.sim=100, parallel = FALSE)
## End(Not run)
# Example for n.cat = 2 (i.e. Spatial probit regression)
# True parameter
vec.beta <- c(1, 2, 1, 0, -1); sigmasq <- 0.5; rho <- 0.6; radius <- 5
vec.par <- c(vec.beta, sigmasq, rho)
# Coordinate matrix
n.cat <- 2 ; n.lati <- n.long <- 40
n.site <- n.lati * n.long
mat.lattice <- cbind(rep(1:n.lati, n.long), rep(1:n.long, each=n.lati))
mat.dist <- as.matrix(dist(mat.lattice, upper=TRUE, diag=TRUE))
mat.cov <- sigmasq * rho^mat.dist
set.seed(123)
# Generate regression (design) matrix with intercept
mat.X <- cbind(rep(1, n.site),scale(matrix(rnorm(n.site*(length(vec.beta)-1)),nrow=n.site)))
vec.Z <- t(chol(mat.cov)) %*% rnorm(n.site) + mat.X %*% vec.beta
vec.epsilon <- diag(sqrt(1-sigmasq), n.site) %*% rnorm(n.site)
vec.ylat <- as.numeric(vec.Z + vec.epsilon)
# Convert to the vector of observation
vec.yobs <- func.obs.ord(vec.ylat, vec.alpha=c(-Inf,0,Inf))
## Not run:
probit.example <- func.cle.ord(vec.yobs, mat.X, mat.lattice, radius, n.cat,
n.sim=100, parallel = TRUE, n.core = 4)
round(probit.example$vec.par,4)
# beta0 beta1 beta2 beta3 beta4 sigma^2 rho
# 1.0427 2.2250 1.0422 0.0156 -1.1489 0.4402 0.6636
round(probit.example$vec.se,4)
# beta0 beta1 beta2 beta3 beta4 sigma^2 rho
# 0.1198 0.1413 0.0863 0.0523 0.0935 0.1600 0.1263
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.