rstar {likelihoodAsy}

Inference on a scalar function of interest by the r* statistic

Description

This function evaluates the r* statistic for testing of a scalar function of interest.

Usage

rstar(data, thetainit, floglik, fscore=NULL, fpsi, psival, datagen, R=1000, seed=NULL, 
      trace=TRUE, ronly=FALSE,  psidesc=NULL, constr.opt="solnp")

Arguments

data

The data as a list. All the elements required to compute the log likelihood function at a given parameter value should be included in this list.

thetainit

A numerical vector containing the initial value for the parameter of the model. It will be used as starting point in the numerical optimization of the log likelihood function.

floglik

A function which returns the log likelihood function at a given parameter value. In particular, for a certain parameter value contained in a numerical vector theta, a call floglik(theta, data) should return a scalar numerical value, the log likelihood function at theta.

fscore

An optional function which returns the score function at a given parameter value. It must return a numerical vector of the same length of thetainit. For a certain parameter value contained in a numerical vector theta, a call fscore(theta, data) should return the gradient of the log likelihood function at theta. Default is NULL, implying that numerical differentiation will be employed.

fpsi

A function which specifies the parameter of interest. A call fpsi(theta) should return a scalar value.

psival

A numerical scalar value containing the value of the parameter of interest under testing.

datagen

A function which simulates a data set. A call datagen(theta, data) will generate a copy of the data list, with the values of the response variable replaced by a set of values simulated from the parametric statistical model assumed for the response variable.

R

The number of Monte Carlo replicates used for computing the r* statistic. A positive integer, default is 1000.

seed

Optional positive integer, the random seed for the Monte Carlo computation. Default is NULL.

trace

Logical. When set to TRUE will cause some information on the computation to be printed. Default is FALSE.

ronly

Logical. If set to TRUE the computation of the r* statistic will be skipped, and only the value of the signed likelihood ratio test statistic r will be returned by the procedure, without any Monte Carlo computation. Default is FALSE.

psidesc

An optional character string describing the nature of the parameter of interest. Default is NULL.

constr.opt

Constrained optimizer used for maximizing the log likelihood function under the null hypothesis. Possible values are "solnp" or "alabama", with the former employing the solnp function from package Rsolnp and the latter the constrOptim.nl from the package alabama. Defauls is "solnp".

Details

The function computes the r* statistic proposed by Skovgaard (1996) for accurate computation of the asymptotic distribution of the signed likelihood ratio test for a scalar function of interest.

The function requires the user to provide three functions defining the log likelihood function, the scalar parametric function of interest, and a function for generating a data set from the assumed statistical model. A further function returning the gradient of the log likelihood is not required, but if provided it will speed up the computation.

When ronly = TRUE the function returns the value of the signed likelihood ratio test statistic r only.

The function handles also one-parameter models.

Value

The returned value is an object of class "rstar", containing the following components:

r

The observed value the signed likelihodo ratio test statistic r for testing fpsi(theta)=psival.

NP, INF

The Nuisance Parameter adjustment (NP) and the Information adjustment (INF) from the decomposition of the r*-r adjustment. The former is not computed for one-parameter models. Neither one is computed when ronly = TRUE.

rs

The observed value of the r* statistic. Not computed when ronly = TRUE.

theta.hat

The maximum likelihood estimate of the parameter theta, the argument of the floglik, fscore, datagen and fpsi functions.

info.hat

The observed information matrix evaluated at theta.hat. Not computed when ronly = TRUE.

se.theta.hat

The estimated standard error of theta.hat. Not computed when ronly = TRUE.

psi.hat

The parameter of interest evaluated at theta.hat.

se.psi.hat

The estimated standard error for the parameter of interest. Not computed when ronly = TRUE.

theta.hyp

The constrained estimate of the parameter, under the null hypothesis fpsi(theta)=psival.

psi.hyp

The value under testing, equal to psival.

seed

Random seed used for Monte Carlo trials. Not returned when ronly = TRUE.

psidesc

A character string describing the nature of the parameter of interest.

R

Number of Monte Carlo replicates used for computing the r* statistic. Not returned when ronly = TRUE.

There are print and summary methods for this class.

References

The method implemented in this function was proposed in

Skovgaard, I.M. (1996) An explicit large-deviation approximation to one-parameter tests. Bernoulli, 2, 145–165.

For a general review

Severini, T.A. (2000). Likelihood Methods in Statistics. Oxford University Press.

See Also

rstar.ci.

Examples

# Autoregressive model of order 1
# We use the lh data from MASS
library(MASS)
data(lh)
dat.y <- list(y=as.numeric(lh))
# First let us define the function returning the log likelihood function
# We employ careful parameterizations for the correlation and variance to
# avoid numerical problems
likAR1 <- function(theta, data)
{ 
  y <- data$y
  mu <- theta[1]
  phi <- theta[2] ### phi is log(sigma) 
  sigma2 <- exp(phi*2)
  z <- theta[3]   ### z is Fisher'z transform for rho
  rho <- (exp(2*z)-1) / (1 + exp(2*z))
  n <- length(y)
  Gamma1 <- diag(1+c(0,rep(rho^2,n-2),0))
  for(i in 2:n)
    Gamma1[i,i-1]<- Gamma1[i-1,i] <- -rho 
  lik <- -n/2 * log(sigma2) + 0.5 * log(1-rho^2) -1/(2*sigma2) * 
        mahalanobis(y, rep(mu,n), Gamma1, inverted = TRUE)
  return(lik)
}
# We need a function for simulating a data set
genDataAR1 <- function(theta, data)  
{
  out <- data
  mu <- theta[1]
  sigma <- exp(theta[2])
  z <- theta[3]
  rho <- (exp(2*z)-1) / (1 + exp(2*z))
  n <- length(data$y)
  y <- rep(0,n)
  y[1] <- rnorm(1,mu,s=sigma*sqrt(1/(1-rho^2)))
  for(i in 2:n)
    y[i] <- mu + rho * (y[i-1]-mu) + rnorm(1) * sigma 
  out$y <- y 
  return(out)
}
# For inference on the mean parameter we need a function returning the first component of theta
psifcn.mu <- function(theta) theta[1]
# Now we can call the function
rs.mu <- rstar(dat.y, c(0,0,0), likAR1, fpsi=psifcn.mu, psival=2, datagen=genDataAR1, R=1000, 
               trace=TRUE, psidesc="mean parameter")
summary(rs.mu)