wws: Wald, Wilks and Score tests
In smovie: Some Movies to Illustrate Concepts in Statistics

View source: R/wald_wilks_score.R

wws	R Documentation

Wald, Wilks and Score tests

Description

A movie to illustrate the nature of the Wald, Wilks and score likelihood-based test statistics, for a model with a scalar unknown parameter \theta. The user can change the value of the parameter under a simple null hypothesis and observe the effect on the test statistics and (approximate) p-values associated with the tests of this hypothesis against the general alternative. The user can specify their own log-likelihood or use one of two in-built examples.

Usage

wws(
  model = c("norm", "binom"),
  theta_range = NULL,
  ...,
  mult = 3,
  theta0 = if (!is.null(theta_range)) sum(c(0.25, 0.75) * theta_range) else NULL,
  panel_plot = TRUE,
  hscale = NA,
  vscale = hscale,
  delta_theta0 = if (!is.null(theta_range)) abs(diff(theta_range))/20 else NULL,
  theta_mle = NULL,
  loglik = NULL,
  alg_score = NULL,
  alg_obs_info = NULL,
  digits = 3
)

Arguments

`model`	A character scalar. Name of the the distribution on which one of two in-built examples are based. If `model = "norm"` then the setting is a random sample of size `n` from a normal distribution with unknown mean `mu` = `\theta` and known standard deviation `sigma`. If `model = "binom"` then the setting is a random sample from a Bernoulli distribution with unknown success probability `\theta`. The behaviour of these examples can be controlled using arguments supplied via `...`. In particular, the data can be supplied using `data`. If `model = "norm"` then `n`, `mu`, and `sigma` can also be chosen. The default cases for these examples are: `model = "norm"`: `n` = 10, `mu` = 0, `sigma` = 1 and `data` contains a sample of a sample of size `n` simulated, using `Normal`, from a normal distribution with mean `mu` and standard deviation `sigma`. `model = "binom"`: `data = c(7, 13)`, that is, 7 successes and 13 failures observed in 20 trials. For the purposes of this movie there must be at least one success and at least one failure.
`theta_range`	A numeric vector of length 2. The range of values of `\theta` over which to plot the log-likelihood. If `theta_range` is not supplied then the argument `mult` is used to set the range automatically.
`...`	Additional arguments to be passed to `loglik`, `alg_score` and `alg_obs_info` if `loglik` is supplied, or to functions functions relating to the in-built examples otherwise. See the description of `model` above for details.
`mult`	A positive numeric scalar. If `theta_range` is not supplied then an interval of width 2 x `mult` standard errors centred on `theta_mle` is used. If `model = "binom"` then `theta_range` is truncated to (0,1) if necessary.
`theta0`	A numeric scalar. The value of `\theta` under the null hypothesis to use at the start of the movie.
`panel_plot`	A logical parameter that determines whether the plot is placed inside the panel (`TRUE`) or in the standard graphics window (`FALSE`). If the plot is to be placed inside the panel then the tkrplot library is required.
`hscale`, `vscale`	Numeric scalars. Scaling parameters for the size of the plot when `panel_plot = TRUE`. The default values are 1.4 on Unix platforms and 2 on Windows platforms.
`delta_theta0`	A numeric scalar. The amount by which the value of `theta0` is increased (or decreased) after one click of the + (or -) button in the parameter window.
`theta_mle`	A numeric scalar. The user may use this to supply the value of the maximum likelihood estimate (MLE) of `\theta`. Otherwise, `optim` is used to search for the MLE, using `theta0` as the initial value and `theta_range` as bounds within which to search.
`loglik`	An R function, vectorised with respect to its first argument, that returns the value of the log-likelihood (up to an additive constant). The movie will not work if the observed information is not finite at the maximum likelihood estimate.
`alg_score`	A R function that returns the score function, that is, the derivative of `loglik` with respect to `\theta`.
`alg_obs_info`	A R function that returns the observed information that is, the negated second derivative of `loglik` with respect to `\theta`.
`digits`	An integer indicating the number of significant digits to be used in the displayed values of the test statistics and p-values. See `signif`.

Details

The Wald, Wilks (or likelihood ratio) and Score tests are asymptotically equivalent tests of a simple hypothesis that a parameter of interest \theta is equal to a particular value \theta_0. The test statistics are all based on the log-likelihood l(\theta for \theta but they differ in the way that they measure the distance between the maximum likelihood estimate (MLE) of \theta and \theta_0. The Wilks statistic is the amount by which the log-likelihood evaluated \theta_0 is smaller than the log-likelihood evaluated at the MLE. The Walk statistics is based on the absolute difference between the MLE and \theta_0. The score test is based on the gradient of the log-likelihood (the score function) at \theta_0. For details see Azzalini (1996).

This movie illustrates the differences between the test statistics for simple models with a single scalar parameter. In the (default) normal example the three test statistics coincide. This is not true in general, as shown by the other in-built example (distn = "binom").

A user-supplied log-likelihood can be provided via loglik.

Value

Nothing is returned, only the animation is produced.

References

Azzalini, A. (1996) Statistical Inference Based on the Likelihood, Chapman & Hall / CRC, London.

Examples

# N(theta, 1) example, test statistics equivalent
wws(theta0 = 0.8)

# binomial(20, theta) example, test statistics similar
wws(theta0 = 0.5, model = "binom")

# binomial(20, theta) example, test statistic rather different
# for theta0 distant from theta_mle
wws(theta0 = 0.9, model = "binom", data = c(19, 1), theta_range = c(0.1, 0.99))

# binomial(2000, theta) example, test statistics very similar
wws(theta0 = 0.5, model = "binom", data = c(1000, 1000))

set.seed(47)
x <- rnorm(10)
wws(theta0 = 0.2, model = "norm", theta_range = c(-1, 1))

# Log-likelihood for a binomial experiment (up to an additive constant)
bin_loglik <- function(p, n_success, n_failure) {
  return(n_success * log(p) + n_failure * log(1 - p))
}

wws(loglik = bin_loglik, theta0 = 0.5, theta_range = c(0.1, 0.7),
    theta_mle = 7 / 20, n_success = 7, n_failure = 13)

bin_alg_score <- function(p, n_success, n_failure) {
  return(n_success / p - n_failure / (1 - p))
}
bin_alg_obs_info <- function(p, n_success, n_failure) {
  return(n_success / p ^ 2 + n_failure / (1 - p) ^ 2)
}
wws(loglik = bin_loglik, theta0 = 0.5, theta_range = c(0.1, 0.7),
    theta_mle = 7 / 20, n_success = 7, n_failure = 13,
    alg_score = bin_alg_score, alg_obs_info = bin_alg_obs_info)

smovie documentation built on May 29, 2024, 10:28 a.m.