Composite Likelihood Calculation for Spatial Ordinal Data without Replications (for implmentation)

Description

cl Calculate the negative composite log-likelihood value and score function for a particular subject given parameter value and other input variables.

Usage

cl(theta, y, X, dwdv, cmwdv, lt, wn, base, J, p)

Arguments

Value

cl returns a list: composite log-likelihood value and a vector of first-order partial derivatives for theta.

Composite Likelihood Calculation for Replciations of Spatial Ordinal Data (for illustration)

Description

cl.rord Calculate the negative composite log-likelihood value for replications of spatial ordinal data at given value of parameter value. Note that this function is not directly used in cle.rord but illustration only.

Usage

cl.rord(theta, response, covar, location, radius = 4)

Arguments

Value

cl.rord returns a list: negative composite log-likelihood, a vector of first-order partial derivatives for theta.

Examples

set.seed(1203)
n.subject <- 10
n.lat <- n.lon <- 10
n.site <- n.lat*n.lon

beta <- c(1,2,-1) # First 1 here is the intercept
midalpha <- c(1.15, 2.18) ; sigma2 <- 0.7 ; phi <- 0.8

true = c(midalpha,beta,sigma2,phi)

Xi = rnorm(n.subject,0,1) ; Xj <- rbinom(n.site,1,0.6)

 VV <- matrix(NA, nrow = n.subject*n.site, ncol = 3)

 for(i in 1:n.subject){ for(j in 1:n.site){
     VV[(i-1)*n.site+j,] <- c(1,Xi[i],Xj[j])
       }
 }

location = cbind(rep(seq(1,n.lat,length=n.lat),n.lat),rep(1:n.lon, each=n.lon))
sim.data <- sim.rord(n.subject, n.site, n.rep = 2, midalpha, beta, sigma2, phi, covar=VV, location)

cl.rord(theta=true,response=sim.data[[1]], covar=VV, location, radius = 4)

Composite Likelihood Calculation for Spatial Ordinal Data without Replications (for implmentation)

Description

cl_l Calculate the negative composite log-likelihood value for a particular subject given parameter value and other input variables

Usage

cl_l(theta, y, X, dwdv, cmwdv, lt, wn, base, J, p)

Arguments

Value

cl returns a list: composite log-likelihood value and a vector of first-order partial derivatives for theta.

Composite Likelihood Estimation for Replciations of Spatial Ordinal Data

Description

cle.rord Estimate parameters (including regression coefficient and cutoff) for replications of spatial ordinal data using pairwise likelihood approach.

Usage

cle.rord(
  response,
  covar,
  location,
  radius = 4,
  n.sim = 100,
  output = TRUE,
  SE = TRUE,
  parallel = FALSE,
  n.core = max(detectCores()/2, 1),
  ini.sp = c(0.5, 0.5),
  est.method = TRUE,
  maxiter = 100,
  rtol = 1e-06,
  factr = 1e+07
)

Arguments

Details

Given vector of ordinal responses, the design matrix, spatial location for sites, weight radius (for pair selection), and the prespecified number of simulation used for estimating the Godambe information matrix. Initial estimate is obtained by fitting model without spatial dependence (using MASS::polr()) and optional guess of spatial parameters. The function first estimates parameters of interest by either solving the gradient of composite log-likelihood using rootSolve::multiroot() or maximize the composite log-likelihood by optim(..., method="L-BFGS-B"). The asymptotic covariance matrix and standard error of parameters are then estimated by parametric boostrapping. Although the default root solving option is typically more efficient, it may encounter runtime error if negative value of \phi is evaluated (and L-BFGS-B approach should be used).

Value

cle.rord returns a list contains:

vec.par: a vector of estimator for \theta=(\alpha,\beta,\phi,\sigma^2);

vec.se: a vector of standard error for the estimator;

mat.asyvar: estimated asymptotic covariance matrix H^{-1}(\theta)J(\theta)H^{-1}(\theta) for the estimator;

mat.Hessian: Hessian matrix at the parameter estimate;

mat.J: Sensitivity matrix estimated by parametric boostrapping; and

CLIC: Composite likelihood information criterion (see help manual of clic() for detail).

Examples

set.seed(1228)
n.subject <- 20
n.lat <- n.lon <- 10
n.site <- n.lat*n.lon

beta <- c(1,2,-1) # First 1 here is the intercept
midalpha <- c(1.15, 2.18) ; phi <- 0.6 ; sigma2 <- 0.7

true <- c(midalpha,beta,phi,sigma2)

Xi <- rnorm(n.subject,0,1) ; Xj <- rbinom(n.site,1,0.6)

 VV <- matrix(NA, nrow = n.subject*n.site, ncol = 3)

 for(i in 1:n.subject){ for(j in 1:n.site){
     VV[(i-1)*n.site+j,] <- c(1,Xi[i],Xj[j])
       }
 }

location <- cbind(rep(seq(1,n.lat,length=n.lat),n.lat),rep(1:n.lon, each=n.lon))
sim.data <- sim.rord(n.subject, n.site, n.rep = 2, midalpha, beta, phi, sigma2, covar=VV, location)


options(digits=3)
result <- cle.rord(response=sim.data[[1]], covar=VV,
          location = location ,radius = 4, n.sim = 100, output = TRUE, parallel=TRUE, n.core =2)
result$vec.par
# alpha2  alpha3   beta0   beta1   beta2     phi sigma^2
# 1.249   2.319   1.169   1.990  -1.000   0.668   0.678

result$vec.se
# alpha2  alpha3   beta0   beta1   beta2     phi sigma^2
# 0.0704  0.1201  0.1370  0.2272  0.0767  0.0346  0.1050

Composite likelihood Information Criterion

Description

clic Calculating the Composite likelihood information criterion proposed by Varin and Vidoni (2005)

Usage

clic(logCL, mat.hessian, mat.J)

Arguments

Details

Varin and Vidoni (2005) proposed the information criterion in the form: -2*logCL(theta) + 2*trace(H^{-1}(\theta)J(\theta))

Value

CLIC: Composite likelihood information criterion proposed by Varin and Vidoni (2005)

clic: Composite likelihood information criterion proposed by Varin and Vidoni (2005)

References

Varin, C. and Vidoni, P. (2005) A note on composite likelihood inference and model selection. Biometrika 92: 519–528.

Composite Likelihood Calculation for Spatial Ordinal Data without Replications (for implmentation)

Description

Composite Likelihood Calculation for Spatial Ordinal Data without Replications (for implmentation)

Usage

## S3 method for class 'list'
merge(x, y = NULL, mergeUnnamed = TRUE)

Arguments

Value

merge.list returns a list: a merged list

Simulation of Replciations of Spatial Ordinal Data

Description

sim.rord Simulate replications of spatial ordinal data

Usage

sim.rord(
  n.subject,
  n.site,
  n.rep = 100,
  midalpha,
  beta,
  phi,
  sigma2,
  covar,
  location
)

Arguments

Value

sim.rord returns a list (length n.rep) of matrix (n.subject*n.site) with the underlying parameter as inputs.

Examples

set.seed(1203)
n.subject <- 100
n.lat <- n.lon <- 10
n.site <- n.lat*n.lon

beta <- c(1,2,-1) # First 1 here is the intercept
midalpha <- c(1.15, 2.18) ; phi <- 0.8 ; sigma2 <- 0.7

true <- c(midalpha,beta,sigma2,phi)

Xi <- rnorm(n.subject,0,1) ; Xj <- rbinom(n.site,1,0.6)

 VV <- matrix(NA, nrow = n.subject*n.site, ncol = 3)

 for(i in 1:n.subject){ for(j in 1:n.site){
     VV[(i-1)*n.site+j,] <- c(1,Xi[i],Xj[j])
       }
 }

location <- cbind(rep(seq(1,n.lat,length=n.lat),n.lat),rep(1:n.lon, each=n.lon))
sim.data <- sim.rord(n.subject, n.site, n.rep = 2, midalpha, beta, phi, sigma2, covar=VV, location)

length(sim.data)
head(sim.data[[1]])
dim(sim.data[[1]])
hist(sim.data[[1]])

Surrogate Residuals for Replciations of Spatial Ordinal Data

Description

surrogate.residual simulate the surrogate residual with the the given parameter value and covariate for model diagnostics.

Usage

surrogate.residual(
  response,
  covar,
  location,
  seed = NULL,
  midalpha,
  beta,
  sigma2,
  phi,
  burn.in = 20,
  output = TRUE
)

Arguments

Details

Given vector of observed responses, the design matrix, spatial location for sites and parameter value, raw surrogate residuals are simulated using an efficient Gibbs sampling, which can be used for model diagnostics. When the fitted model is correct, the raw surrogate residuals among subjects should follow multivariate normal with mean 0 and covariance Sigma. If the model is correct, residual plot should be close to a null plot or random scatter. For example, it can be used to check the potential missing in covariate, non-linearity of covariate and outliers. In particular for the example below, the residual plot shows that linearity of Xi is adequate for the model.

Value

surrogate.residual returns a (no. spatial site * no. subject) matrix contains raw surrogate residuals with element corresponds to the response matrix.

Examples

set.seed(1228)
n.subject <- 50
n.lat <- n.lon <- 10
n.site <- n.lat*n.lon

beta <- c(1,2,-1) # First 1 here is the intercept
midalpha <- c(1.15, 2.18) ; phi <- 0.6 ; sigma2 <- 0.7

true <- c(midalpha,beta,phi,sigma2)

Xi <- rnorm(n.subject,0,1) ; Xj <- rbinom(n.site,1,0.6)

 VV <- matrix(NA, nrow = n.subject*n.site, ncol = 3)

 for(i in 1:n.subject){ for(j in 1:n.site){
     VV[(i-1)*n.site+j,] <- c(1,Xi[i],Xj[j])
       }
 }

location <- cbind(rep(seq(1,n.lat,length=n.lat),n.lat),rep(1:n.lon, each=n.lon))
response <- sim.rord(n.subject, n.site, n.rep = 1,
midalpha, beta, phi, sigma2, covar=VV, location)[[1]]


# Example for linearity of covariate
sur.resid <- surrogate.residual(response, covar=VV, location, seed =1,
midalpha, beta, sigma2, phi,
burn.in=20, output = TRUE)

scatter.smooth(rep(Xi,each=n.site),c(sur.resid),
main="Surrogate residual against Xi", xlab="Xi", ylab="Surrogate residual",
lpars = list(col = "red", lwd = 3, lty = 2))

abline(h=0, col="blue")