9 Plotting predicted and observed covariances

Once you have fitted a number of alternative models, it's a good idea to inspect how well the model actually fits the data. One way to do this is to compare the covariances between the gene frequencies predicted by the selected model with the covariances computed from the date directly.

The likelihood is computed on the basis of $n-1$ contrasts $\mathbf{y} = (y_1,y_2,\dots,y_{n-1})$ between the observed standardized gene frequencies. The covariance matrix of these contrasts is returned by function covpred. Two arguments are required; the name of the migration model FUN, the parameter vector theta.

Function covobs estimates the covariances from the observations directly using the estimator

$$\hat c_{ij}=\frac1{n_{loci}} \sum_{k=1}^{n_{loci}} \frac{(p_{i,k}-\bar p_k)(p_{j,k}-\bar p_k)}{\bar p_k(1-\bar p_k)}.$$ (4)

The only required argument is the matrix of gene frequencies p. The covariances can be estimated based on weighted gene frequency means at each loci, by providing a vector w of weights to covobs as a second optional argument.

Figure 2: Predicted and observed covariances plotted using the covplot function

The predicted and observed covariances can be compared through plotting them with the covplot function. Arguments are the computed predicted and observed covariance matrices. The following call

   > covplot(covpred(steppingstone,fit1$par),covobs(p))

produced the plot shown in Figure 2. Two additional optional arguments can be given; mfcol specifying the number of rows and columns in the plot, and file specifying the name of a postscript file (used as alternative output). The interpretation of these type of plots is discussed in Tufto et al. (1998).

Function covpred computes (an approximation) of the covariances under a given model, conditioned on the observed gene frequency means at each locus. The unconditional covariances may also be of interest. These can be computed by function courgeau, called with the migration matrix M and the vector of effective population sizes Ne as arguments. This function rewrites the matrix equation (Tufto et al., 1996, eq. 7) to a system of $n(n+1)/2$ equations (Tufto et al., 1996, eq. A.4 and A.5) and solves these. It may be noted that the returned unconditional covariances (which are exact to the extent that the order of the events in the life cycle can be ignored) can differ greatly from the conditional ones, especially if the long range rate of migration is low such that the gene frequency vector is in one of the states $(0,0,\dots,0)$ or $(1,1,\dots,1)$ for long periods of time.