mspeNERjack {saeMSPE}

Compute MSPE through Jackknife-based MSPE estimation method for Nested error regression model

Description

This function returns MSPE estimator with Jackknife-based MSPE estimation method for Nested error regression model.

Usage

mspeNERjack(ni, X, Y, Xmean, method = 1)

Arguments

ni

(vector). It represents the sample number for every small area.

X

(matrix). Stands for the available auxiliary values.

Y

(vector). It represents the response value for Nested error regression model.

Xmean

(matrix). Stands for the population mean of auxiliary values.

method

The MSPE estimation method to be used. See "Details".

Details

This bias-corrected jackknife MSPE estimator was proposed by J. Jiang and L. S. M. Wan, it covers a fairly general class of mixed models which includes gLMM, mixed logistic model and some of the widely used mixed linear models as special cases.

Default value for method is 1, method = 1 represents the MOM method , method = 2 and method = 3 represents ML and REML method, respectively.

Value

This function returns a list with components:

MSPE

(vector) MSPE estimates for NER model.

bhat

(vector) Estimates of the unknown regression coefficients.

sigvhat2

(numeric) Estimates of the area-specific variance component.

sigehat2

(numeric) Estimates of the random error variance component.

Author(s)

Peiwen Xiao, Xiaohui Liu, Yuzi Liu, Jiming Jiang, and Shaochu Liu

References

M. H. Quenouille. Approximate tests of correlation in time series. Journal of the Royal Statistical Society. Series B (Methodological), 11(1):68-84, 1949.

J. W. Tukey. Bias and confidence in not quite large samples. Annals of Mathematical Statistics, 29(2):614, 1958.

J. Jiang and L. S. M. Wan. A unified jackknife theory for empirical best prediction with m estimation. Annals of Statistics, 30(6):1782-1810, 2002.

Examples

### parameter setting 
Ni = 1000; sigmaX = 1.5; m = 5
beta = c(0.5, 1)
sigma_v2 = 0.8; sigma_e2 = 1
ni = sample(seq(1,10), m,replace = TRUE); n = sum(ni)
p = length(beta)
### population function
pop.model = function(Ni, sigmaX, beta, sigma_v2, sigma_e2, m){
  x = rnorm(m * Ni, 1, sqrt(sigmaX)); v = rnorm(m, 0, sqrt(sigma_v2)); y = numeric(m * Ni)
  theta = numeric(m); kk = 1
  for(i in 1 : m){
    sumx = 0
    for(j in 1:Ni){
      sumx = sumx + x[kk]
      y[kk] = beta[1] + beta[2] * x[kk] + v[i] + rnorm(1, 0, sqrt(sigma_e2))
      kk = kk + 1
    }
    meanx = sumx/Ni
    theta[i] = beta[1] + beta[2] * meanx + v[i]
  }
  group = rep(seq(m), each = Ni)
  x = cbind(rep(1, m*Ni), x)
  data = cbind(x, y, group)
  return(list(data = data, theta = theta))
} 
### sample function
sampleXY = function(Ni, ni, m, Population){
  Indx = c()
  for(i in 1:m){
    Indx = c(Indx, sample(c(((i - 1) * Ni + 1) : (i * Ni)), ni[i]))
  }
  Sample = Population[Indx, ]; Nonsample = Population[-Indx, ]
  return(list(Sample, Nonsample))
} 
### data generation process
Population = pop.model(Ni, sigmaX, beta, sigma_v2, sigma_e2, m)$data
XY = sampleXY(Ni, ni, m, Population)[[1]]
X = XY[, 1:p]
Y = XY[, p+1]
Xmean = matrix(NA, m, p)
for(tt in 1: m){
  Xmean[tt, ] = colMeans(Population[which(Population[,p+2] == tt), 1:p])
}
### mspe result
mspeNERjack(ni, X, Y, Xmean, method = 1)