8 Likelihood ratio tests

Code is also provided for doing likelihood ratio tests between different nested, and perhaps also, non-nested alternative models. This is done, by brute force, by function lrtest by simulating bootstrap data from $H_0$ and computing the likelihood ratio by fitting both $H_0$ and $H_1$ numerically to each bootstrap data set.

Let's say we want to use the island model as our null hypothesis. Function islandmodel simply returns a diagonal matrix with nonzero elements equal to $1-m$:

> islandmodel(.1,10)
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]  0.9  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
 [2,]  0.0  0.9  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
 [3,]  0.0  0.0  0.9  0.0  0.0  0.0  0.0  0.0  0.0   0.0
 [4,]  0.0  0.0  0.0  0.9  0.0  0.0  0.0  0.0  0.0   0.0
 [5,]  0.0  0.0  0.0  0.0  0.9  0.0  0.0  0.0  0.0   0.0
 [6,]  0.0  0.0  0.0  0.0  0.0  0.9  0.0  0.0  0.0   0.0
 [7,]  0.0  0.0  0.0  0.0  0.0  0.0  0.9  0.0  0.0   0.0
 [8,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.9  0.0   0.0
 [9,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.9   0.0
[10,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.9

Note that this model is nested within the steppingstone model.

lrtest takes six arguments; H0 and H1 which specifies the names of the functions specifying migration pattern under the null and the alternative hypothesis, the gene frequency data matrix p, optionally a vector theta specifying from which parameter values under $H_0$ to simulate from, and nboots specifying the number of bootstrap replicates of the likelihood ratio to compute. After some time (maybe some hours), lrtest, as a side-effect, displays a summary of the test results:

   > test <- lrtest(islandmodel,steppingstone,p,nboots=50)
   Likelihood-ratio test summary

   Observed value of the test statistic:  26.8696459394272
   Number of bootstrap samples:  10
   Significance probability:  0

In addition, a list containing a number of objects is returned: The $H0 and $H1 components contains the lists returned by fitmodel when $H_0$ and $H_1$ are fitted to the observed data. Thus, for our example data set, the maximum likelihood estimates under the islandmodel is stored in

   > test$H0$par
   [1] 0.04816419

As we see, under the fitted (false) islandmodel, the rate of migration from the `mainland' has to be as large as 0.04 to account for the small between-population differentiation generated by the high rate of migration between neighbouring population in the true (steppingstone) model.

Figure 1: Histogram of simulated log likelihood ratios

Note that, unless the optional argument theta is given, lrtest simulates the distribution of the likelihood ratio from the maximum likelihood estimates under $H_0$. However, if there is a large amount of bias in the estimates (see e.g. Tufto et al., 1998), it may be adviceable to first do some sort of bias correction and supply lrtest with the bias corrected estimate through the theta argument.

The simulated (log) likelihood ratios are stored in the $lr component of the list returned by lrtest. A histogram of these can be produced with S-Plus' hist function (see Figure 1). As can be seen from the figure, the distribution is certainly far from $\chi^2$.

Nothing prevents the user from carrying out likelihood ratio tests between non-nested models, although the interpretation of such tests may be questionable from a theoretical point of view (see e.g. Vuong, 1989). It should be kept in mind that such tests can be carried out in two directions. It is therefore up to the user to decide which model represents the null hypothesis and which represent the alternative.