# Validating a model for a set of DNA sequences

Consider a set of DNA sequences i.e. strings composed from the four letters A, C, G, T.

We model these sequences as a random vector, $\mathbf{X}=(X_1,\dots, X_n)$. Each member of the set of sequences is a realization of this random vector.

In this case the $X_i$ only take the values A, C, G, T. More generally let $r$ be the number of letters.

We define a family of probability distributions for X. Our aim is to find the probability distribution that is the best fit for the data.

Our assumption is that the variables $X_i$ can be divided into disjoint subsets, such that variables in each subset are dependent among themselves, and the subsets are mutually independent.

It suffices to define the distribution on one of these subsets. This distribution is defined as follows. Let $\mathbf{S}=(X_{j_1},\dots, X_{j_k})$ be a sub-vector of $\mathbf{X}$ for some subset $\{j_1,\dots, j_k\} \in \{1,2\dots, n\}$. Then $\mathbf{S}$ takes $r^k$ possible values.

Then we define a categorical probability distribution on $\mathbf{S}$.

$P(S=(s_{1}, s_{2}, \dots, s_{k})) = p_{s_{1}s_{2}\dots s_{k}}$

So the distribution of $\mathbf{S}$ is defined by $r^k$ probabilities.

Observe that the collection of random variables $\{X_{j_1},\dots, X_{j_k}\}$ with this distribution is not independent in general.

We use MDL to find a suitable model for the data, in combination with a genetic search algorithm. This gives good results when applied to simulated data, with the caveat that it requires a sufficiently large data set to converge to the correct model.

However, the important question is how well this model validates against real data.

I believe this can be divided into the question of how much explanatory power this model has when applied to real data, as well as what predictive power it has.

I thought about the explanatory power issue, and I'm not sure what can be usefully done. Maybe this model just has too many parameters from the perspective of model fitting. However, if anyone has suggestions, please let me know.

As regards prediction, I think I can use cross-validation as well as permutation resampling.

For cross-validation, I can divide the data into K portions, and set aside in turn each of the portions to use for prediction, and select the model with the other K-1 portions. So, the question is then how one should do prediction.

Since the aim here is to check whether the predicted model for the validation set is the correct one. Usually, one needs a decision rule, which (I think) in this case would be of the form likelihood > value in this case, since the higher the likelihood, the more likely the data comes from the model. I'm not sure how to formalize this, or if one can justify it from some theory. I've found precisely zero references of anyone doing this, which suggests I may be doing something wrong. I've looked at some decision theory stuff but most results seem to assume conditions are met that are not true in this case. I'm also not sure how to set the cutuff. One possible approach would be to use the empirical distribution of the training set and set some cutoff based on that.

A related alternative approach would be to generate ROC curves, and then generate some summary statistics from that, like Area Under Curves.

For permulation resampling, I can take the data set and randomly permute it across the strings. So the different values of $X_i$ in the different strings will be exchanged among themselves. This will produce a data set which has the same marginal distributions as the original data set. but where the individual variables are uncorrelated.

One can then check to see how well the model distinguishes this permuted data set with the original data set. It seems to me this if there is clear separation, i.e. the likelihoods of the permuted data set are well below the likelihoods of the original data set, then that is evidence that the original data does have correlation structure. I'm not sure how to quantify this separation.