rstar: Inference on a scalar function of interest by the r*...
In likelihoodAsy: Functions for Likelihood Asymptotics

Description Usage Arguments Details Value References See Also Examples

View source: R/likelihoodAsy.R

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

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

`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"`.

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.

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.

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.

rstar.ci.

# 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)