In my last post, I outlined the ideas of statistical inference: that we have a model, and we have a philosophical scheme that tells us what maths to do to estimate the parameters of the model. But this is not the same as doing this in practice. The equations are often complex, and difficult or impossible to solve (anyone fancy trying to integrate in 1000 dimensions?). Because of this, a lot of methods have been devised for estimating the parameters, and particularly estimating their variability (which, from now on, I will write about as estimating the standard error).
For maximum likelihood, the problem of getting the estimates is an optimisation problem: find the set of parameters that give the largest likelihood. For the simple linear regression, this is easy, because the equations were worked out a couple of hundred years ago, by some French mathematicians trying to work out where Paris was. For more complex models, such as generalised linear models, there are algorithms that have been shown to efficiently iterate towards the ML estimates. In more complex models, more general or complicated search algorithms might be needed (for example, there is the EM algorithm which is useful for dealing with missing data. So useful that the missing data is sometimes invented, in order to use it).
These algorithms only give point estimates. We also want to know the standard error, and there are methods for doing this. Again, it may be that the equations can be derived directly, or can be found with a simple numerical search. For less standard problems, more complex algorithms can be used. Two popular ones are called the bootstrap and the jackknife. They both work by re-sampling the data, and using that variation to estimate the standard error. What is less well known is that they can also be used to estimate the bias in the point estimates.
Even more complex problems have even more complex solutions. For Bayesians, the most common method for fitting their models is a technique called MCMC. This estimates the parameters (actually their full distribution) by simulation. The technique is not simple to explain, but it is used because it is very flexible and for a lot of problems works well in practice. However, this is not the same as a Bayesian analysis: MCMC can be used by frequentists as well, and other complicated methods can also be used to fit Bayesian models (e.g. sequential importance sampling).
So, now we have a model fitted. There is still one part of the process that is often done: model selection. This is not always done, and is perhaps done more often than it should. It also creates a lot of heated debates, which may surprise some of you. That's because model selection is often disguised as hypothesis testing.
A general description of the role of model selection is that it is a way of choosing between different mathematical models for the data, to find the best one. For frequentists there are two principal ways of doing this: null hypothesis statistical testing and using information criteria to compare models.
Null hypothesis statistical testing (NHST) is the method everyone is taught to test hypotheses. What one does is set up a null hypothesis, for example that there is no difference between two groups, and then test whether this hypothesis is supported by the data. If it is, then the null hypothesis is accepted (strictly, it's not rejected), otherwise it is rejected. In one common formulation, whenever a null hypothesis is set up, a more general alternative hypothesis is also proposed: rejecting the null hypothesis means accepting the alternative hypothesis.
Where is the model selection in this? We can take the Harry Potter example, and test the null hypothesis that the change in book length over time is zero. The relevant part of the mathematical model is
mi = a + b Xi
where mi is the expected number of pages, and Xi is the year of publication. The null hypothesis is that b=0. But this is equivalent to this model
mi = a
And mi = a + b Xi is the alternative model for the alternative hypothesis. Hence, hypothesis testing becomes a problem of comparing two models. The models are compared by calculating the amount of variation in the data explained by the two models, and seeing if the larger model (i.e. the one with the effect of time) explains a significantly greater amount of variation.
The alternative method for model comparison uses information criteria. These are measures of model adequacy: they consist of a measure of model fit (the deviance: the smaller the deviance, the closer the fitted model is to the data), and a penalisation term for complexity (the more parameters in the model, the higher this penalisation). The different models being compared are fitted to the data, and the relevant information criterion (e.g. AIC or BIC) is calculated for each one. The model with the lowest criterion is declared the best. If there are several models with similar criteria, they are sometimes all examined, and one of them chosen to be the best, for other reasons (e.g. because it makes sense scientifically).
Both of these methods have their counterparts in Bayesian analysis: hypothesis testing can be carried out using Bayes factors, for instance, and there is an information criterion for Bayesians (DIC. Alas, not the Bayesian information criterion!).
So, there we have it. OF course, this is not all there is in statistical inference: for example, I have not dealt with model checking (i.e. asking whether the model actually fits the data well), and there are many details I have left out. But I hope this give a scheme that separates out the different parts of fitting statistical models, and reduces some of the confusion. To summarise, here is a list of the parts of the process, and examples of the actual methods used in each part:
Read more!
Wednesday, 1 August 2007
Statistical Modelling: The Bits Pt. 2
Monday, 30 July 2007
Statistical Modelling: The Bits Pt. 1
It's a sad fact of life that it's too short, so we don't have time to learn everything that we need to know. For most biologists (and probably most scientists, and indeed social scientists), one of the things they should learn to do science well is statistics. Of course, the poor dears are too busy learning about things like biochemistry or PCR, to have time to learn about the important stuff.
Because of this, there is a lot they don't know, and also a lot of confusion about the bits they do know. I therefore thought it was worth writing a few posts about some of these issues, to try and dispel at least some of the confusion.
Now we're below the fold, I should admit that I don't mind too much the lack of knowledge - there are large swathes of statistics that I don't know about either (I was too busy admiring the paper aeroplane building skills of my fellow biochemistry students). The confusion is more of a concern, and I think it is largely because biologists get little training in statistics as undergraduates, so don't have the grounding in the theory when they start doing research. And, yes, I am generalising - some biologists do become very good statisticians. But there is still plenty of confusion.
A lot of basic confusion comes from not understanding the different parts of using models to analyse data. For the moment I am ignoring things like model checking, and how to interpret the results. This is not because they are unimportant, but just because I want to write a blog entry, not a monograph. Anyway, There are, I think, 4 parts that can be separated out:
Point 1 is often passed over, and points 3 and 4 are often confused with point 2. I will deal with the first two points in this post, and treat the others later.
To help disentangle the different parts, it is useful to have an example. For this, we can use a simple linear regression, such as the one I used here. The mathematical model is self-explanatory. For the regression, it can be written as
mi = a + b Xi
Yi ~ N(mi, s2).
The first line is the equation for a straight line. Xi is the covariate (year, in the example). This is then related to an expected value, mi. The second line says that each data point, Yi, is normally distributed with mean mi and variance s2.
Note that the data are assumed to be drawn from some distribution. Hence, they are random. However, they depend on parameters, in this case the mean and variance. Here the mean is modelled further, using an equation, so it is deterministic. i.e. if we knew the parameters, we would know the true value of the mean.
More complex analyses also have a model underneath them: the usual ANOVA actually has a model that is almost the same. The mean can be a function of several parameters, and might not be the mean, but just a parameter that controls the likelihood. The variance might also be modelled. And the parameters that are modelled could themselves be functions of more parameters. And all these parameters might themselves be random. Yes, things can get complex.
The problem, then, is to estimate the parameters. There are several schemes for doing this. For the model above, the two commonly used schemes are the frequentist and the Bayesian approaches. Each of these gives us a way to estimate the parameters. They are, in essence, philosophical approaches to understanding what we mean by probability and randomness. In the frequentist scheme, what we observe is random, with fixed parameters, and hence any statistic we calculate (such as the slope of the straight line) is a random function of the data. In contrast, the Bayesian approach is to say that anything we are uncertain about is random. Hence, we treat the data as fixed (because we know it: we have observed it), and the parameters are random.
Both of these philosophical schemes are useful in that they then can be used to construct methods for estimating the parameters of the mathematical model. Both use the likelihood, i.e. the probability of the data given values of the parameters, albeit in different ways. For example, the frequentist scheme tells us that we should find the values of the parameters which give us the maximum value of the likelihood. Hence, this is often spoken about as a maximum likelihood approach (there is another philosophical scheme which leads to the same equations, but to me it always looks like asking Fisher to wave his fiducial wand to make everything mysteriously right. I assume I'm missing something important).
Estimating the parameters is not enough, though. We also want to know how reliable they are. It makes a difference to our estimate that for every year J.K. Rowling writes an extra 300 pages if the range of possible values is 290-310 or 0-600. In the former case, we have a pretty good idea how much time to book off work to read her latest book: in the latter, we have little idea about whether we need to have a headache or major heart surgery. We can summarise the reliability in a variety of ways: a popular one is the probability that the estimate is less than zero. Alternatives include giving a range of likely values (typically a range for which there is a 95% probability that the estimate is within it), or a statistic such as a standard deviation to summarise the variability in the estimate (in the frequentist scheme, this standard deviation is called the standard error).
In my next post, I will describe the model fitting and model selection, i.e. how these schemes are used in practice to make inferences. I suspect this will include a rant about p-values.
Before I stop to work myself up into a frenzy, one last point is worth mentioning. The approaches discussed above assume that the data are normally distributed. If we discard this assumption, we can use different approaches. One, called least squares, is based on optimising a property of the model and data (i.e. minimising the sum of the squares of the differences between the data and the expected values). It turns out that this is the same as maximising the likelihood, i.e. it is mathematically the same as the frequentist approach. I, therfore, use the frequentist interpretation of the equations, and it gets round a lot of fiddly problems. Another, called non-parametric statistics, does away with any distributional assumption, but also throws away part of the data: rather than use the data themselves, the analysis concentrates on using the ranks of the data (i.e whether each Yi is the smallest, second smallest etc.). In some cases, I think this can still have a likelihood interpretation, but I must admit that I haven't checked: it's just a hunch.
Read more!
