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)