Pffft...Luke...I Am Your (Updating) Factor!: A Short Guide to Bayesian Statistics and Bayesian Experimental Design

Bayes Theorem – remember that mentioned way, way back in Lecture 2? No, it isn’t some new age way of predicting who you’ll be romantically involved with this winter, but there is a field of inference that comes from this theorem that plays into Bayesian statistics, and that is a subfield of statistics many scientists should be paying attention to.

Up until this point, we’ve basically been learning more frequentist statistics than Bayesian statistics (i.e., heavy on the linear regression, chi-squares, correlations, less so on multiple comparisons, etc.). This is evident by our HistoryStats projects: we’ve been looking at the lives and work of some of the founders of frequentists’ school of thought like Neyman, Pearson, and Wald. How do we best describe these  frequentist statisticians? Well, let’s take a simple, intuitive analogy described by this StackExchange forum. According to “user28” having a frequentist frame of mind is like hearing the phone go off, referring to a model upon which helps you identify the area of your home that phone is going off to make the inference on where the phone is. Having a Bayesian frame of mind means you may have that model in mind, but you also take into account places where you’ve mistakenly left the phone in the past. Simply, frequentists believe that data is a frequency, or a repeatable random sample, while Bayesians believe that data is observed from a real sample. Furthermore, frequentists believe that parameters are fixed, whereas Bayesians believe the parameters to be unknown but can be described by probabilities. (So…that would make Fisher’s maximum likelihoods a closet Bayesian statistic, wouldn’t it?).

TJ gave us some great examples of Bayes Theorem applied to real life, like the probabilities in clinical trials with cancer treatments. However, we never really got to see how Bayesian inference affects the experiment’s statistics and experimental design.

To understand the experimental design, we need to understand exactly how experimental design is updated or modified by Bayes Theorem, generally. Let’s say you are going to flip a coin 10 times and you suspect a probability distribution to describe these coin flips. Therefore, h would represent the probability of heads, and p(h) would represent the distribution settled on prior to any coin flips. Then the coin is flipped and way more heads come up than usual, say 8 heads. By using Bayesian inference, we need to update our prior belief about the coin – it’s now unfair. So our new beliefs may be modeled like p(h|f) where f is the number of heads experienced in those 10 flips. This abstraction is read as “what is the probability distribution of heads given the number of heads resulting from 10 tosses [in this case 8]?” This seems like a reasonable update as we pare down our hypotheses to fit our experimental data. Mathematically, the update would look like p(h|f) = u(h, f) x p(h) where u(h,f) is an updating factor written out as u(h, f) = (l(f|h))/l(f) where l(f|h) is a likelihood function or the probability we observed 8 heads given the parameters we modeled in the beginning. The denominator of the updating factor is just the likelihood of the data under no conditions. Because Bayesian statistics doesn’t believe parameters are fixed, they can have conditions added to them. Therefore, the likelihood of the data can be written as an integral l(f) = ∫l(f|h)p(h)dh (this is similar to a general expectation value). The denominator turns out to be a weighted average of likelihoods across all possible parameters. Or simply, a ratio that is able to tell you what parameter values are most likely.

Darth Vader: crafty with a lightsaber and some conditional probabilities.

How does this play out in the lab? Let’s take a hypothetical animal trial where dose concentrations of many drugs are tested on large amounts of animals to test their potencies. The lab wants to apply regression analyses to the different drugs based on the specimen they inject the drug into. For experimentation of one drug, the experimental design included six equally spaced doses given to ten mice each; so, 60 animals to test a range of concentrations for one drug. The investigators measured the number of surviving mice one week after drug administration. It turns out that about 90% of mice died at high concentrations of the drug, while 10-20% died at low concentrations of the drug. After each of the experiments, maximum likelihood estimations were used to estimate an LD50 value (or the dose at which the probability of mice dying is 50%). As it turns out, the investigators used results from the first few sets of experiments to predict a distribution for following experiments, in anticipation of constructing an updating factor, as described above. In total, if 50 drugs are tested with similar experimental design, the investigators can use these 50 LD50 values as a sample from a distribution of LD50 values.

Overall, these Bayesian inferences and the statistics are mathematically rooted in Bayes Theorem. This theorem relies on conditional probability. These conditional probabilities make the system easy to update and a noteworthy design for scientists to consider -- because writing grant proposals on frequentist assumptions can be dangerous when we try to predict a model for data without any prior knowledge of the system.  

Limits of Detection (and Credibility)

Methods of Effective Conjugation of Antigens to Nanoparticles as Non-Inflammatory Vaccine Carriers

Nanoparticles are promising carriers for recombinant vaccine antigens, although the physical properties that make some particles good vaccines and others not are still unknown. In this work[1], the authors conjugate model ovalbumin (OVA) antigen to polystyrene nanoparticles for use as a vaccine delivery system. They immunize mice and examine the in vivo responses in order to elucidate an immunological mechanism, however their poor statistical analysis of the resulting data leaves little in the way of conclusions.

Before I get to that, another issue I had with their data was an odd regression curve used to determine concentration from absorbance (Figure 3). In the adjacent figure, the authors used what appear to be serial dilutions of nanoparticles to generate a regression between optical density (OD) at 248 nm and the number of nanoparticles per mL. While their R² value is especially high, as would be expected for a serial dilution, there are no error bars on the points or indications of the number of replicates performed. More seriously, however, is the regression equation they use. In drawing a correlation between the number of nanoparticles in units of 10¹³/mL and the optical density, both numbers with one significant figure, it seemed odd that their regression equation would have an intercept with 4 significant figures, leading me to question the value of their equation and its corresponding R² value. Additionally, the fact that serial dilutions of a solution would have linearly decreasing absorbance values is not necessarily novel or informative, and this figure could have been put in the supplemental information.

More statistically egregious, however, is their table of multiplexed, cytokine-bead array results, expressed as pg/mL (Tables 3 and 4):

Looking only at the first column, one wonders what the difference between 1.4 ± 1.9, or 0 ± 0.21 and “Not detectable” is. Given the absurd concentrations seen in naïve serum and serum from mice injected with NPs, it’s clear that the data is very widely skewed, in which case actually seeing the data points in a graph would be much more helpful. It is also suspicious that these values were the only data in the paper displayed in tabular form, whereas all others were either histograms or bar graphs of some sort. No mention of the number of mice used for this analysis was found in the figure caption, or the results or methods sections either. Not only are their conclusions on immunological mechanism weak, but their poor statistical analysis calls the significance of other data they present into question as well.

[1] Xiang SD, Wilson K, Day S, Fuchsberger M, Plebanski M. Methods of effective conjugation of antigens to nanoparticles as non-inflammatory vaccine carriers. Methods 2013;60:232-41.

Tiny bird headphones and testing model fit: An example use of the extra sum of squares F test

One of the main avenues of research in our lab is quantitative analysis of behavior. One behavior we are particularly interested in is sensorimotor error correction in songbirds. “Sensorimotor error correction” refers to the process by which sensory feedback (such as the auditory feedback of hearing oneself sing) is used to correct a motor behavior (such as singing).

To induce a sensory “error”, our lab fits birds with sets of miniature headphones. While the bird sings, a microphone in the cage records them, and sound processing software will artificially “shift” the pitch of the song up by a couple semitones. This pitch-shifted version song is then played back to the bird through the headphones, virtually in real time. To compensate for the “error” it hears in the auditory feedback, the bird will start singing at a lower pitch (note that if you artificially shift the pitch down, the bird will shift its pitch up).

Previous work had shown that birds will learn to compensate for pitch shift at a faster rate if the shift was small. For very large pitch shifts, they barely learn at all. It’s important to note here that while each bird has its own individual song, they don’t sing it exactly the same way every time. The pitch of a particular note will vary from rendition to rendition, in a normally distributed manner. A technician in our lab hypothesized that it wasn’t the raw amount of sensory error from the pitch shift that influenced learning rate, it was the overlap between the error and the distribution of pitches the bird typically sang at that mattered most.

To test this, the technician decided to use an extra sum of squares f-test.  He created two models (not discussed in detail here to avoid going too egregiously over the word limit, see the paper for more), one of which included parameters for both error size and overlap with the prior distribution, and one that just included error. Then, he took birds with a range of pitch distributions (young adult birds have more variability, and thus a wider distribution of pitches they sing at, while older birds have a narrower distribution). He then tested those birds with a variety of different error sizes via the headphones pitch shift.

The extra sum of squares F test is a way of comparing the fit of two nested models. “Nested” models are models which are identical, but one has additional parameter(s). The extra sum of squares test asks whether the additional parameters reduce the residual error or not. In the case of my labmate, he wanted to know whether there would be less residual error in the error + prior distribution model than the error-only model.

GraphPad’s help page offers another great example of nesting, which may be more intuitive to most biologists:

“If you asked Prism to test whether parameters are different between treatments, then the models are nested. You are comparing a model where Prism finds separate best-fit values for some parameters vs. a model where those parameters are shared among data sets. The second case (sharing) is a simpler version (fewer parameters) than the first case (individual parameters).”

The change in residual sum of squares is divided by the additional degrees of freedom for the extra variables, giving us a mean square. The mean square is then compared to the residual mean square from the full model. An F-test allows us to determine the likelihood of our result, assuming the null is true.

In my labmate’s experiment, the more complex model that included prior distribution overlap significantly reduced the residual error. Check out the full paper here. Here is a longer explanation of extra sum of squares F tests.

Use Prism sample data to understand model fitting

Having trouble fitting a model to your data? Don’t worry because Prism is there to help. The program makes it easy to walk you through the basics of model fitting, and even includes sample data sets for you to practice on. Prism offers the following analyses that can fit lines and curves:

• • Linear regression
• • Deming linear regression (use when both X and Y variables are subject to error)
• • Nonlinear regression
• • Spline and Lowess (for curve fitting without selecting a model)
• • Interpolating from a standard curve

For simplicity’s sake I’ll walk you through how to take advantage of linear and nonlinear regression model fitting using Prism. 

Linear Regression

Linear regression is used when you can describe the data using the equation y=mx+b. Knowing x, you can predict y using the goodness-of-fit line predicted by Prism. To start you can use the sample data Linear Regression – predicting the slope. You should get a data table that looks like this:

Minutes	Control			Treated
1.0	34.	29.	28.	31.	29.	44.
2.0	38.	49.	53.	61.		89.
3.0	57.		55.	78.	99.	77.
4.0	65.	65.	50.	93.	111.	109.
5.0	76.	91.	84.		109.	141.
6.0	79.	93.	98.	134.	145.	129.
7.0	100.	107.	89.	156.	134.	167.
8.0	105.	123.	119.	167.		180.
9.0	121.	143.	134.	178.	192.	175.
10.0	135.	156.		198.	203.	234.

Next it is time to analyze the data: click on the analyze button and then select linear regression.

	Control	Treated
Best-fit values
Slope	12.42 ± 0.5679	17.96 ± 0.8684
Y-intercept when X=0.0	17.42 ± 3.471	28.48 ± 5.446
X-intercept when Y=0.0	-1.402	-1.585
1/slope	0.08050	0.05568
95% Confidence Intervals
Slope	11.26 to 13.59	16.17 to 19.75
Y-intercept when X=0.0	10.28 to 24.55	17.26 to 39.69
X-intercept when Y=0.0	-2.163 to -0.7630	-2.432 to -0.8818
Goodness of Fit
R square	0.9485	0.9448
Sy.x	8.440	13.13
Is slope significantly non-zero?
F	478.5	427.8
DFn, DFd	1.000, 26.00	1.000, 25.00
P value	< 0.0001	< 0.0001
Deviation from zero?	Significant	Significant
Data
Number of X values	10	10
Maximum number of Y replicates	3	3
Total number of values	28	27
Number of missing values	2	3

When determining whether or not the model fits your data, take a look at the R² value. This is called the coefficient of determination and provides you with an idea of how well the best-fit line fits the data. The closer the R² value is to 1, the better the fit.  In this example, we see that R² is equal to .9485 and .9448 for the control and treatment group respectfully. The high fit of this model is further more confirmed by the graph in which there is little deviation of the actual from the expected values of Y according to the model.

Nonlinear Regression

Nonlinear regression is used when the predicted relationship between x and y is not as simple as a linear line. These models use different functions to derive a goodness-of-fit line and may depend on multiple independent variables. We’ll use an enzyme kinetics model to see how nonlinear regression models can be used to fit your data. Use the sample data provided by Prism titles Enzyme Kinetics – Michaelis-Menten. The data table should look like this:

[Substrate]	Enzyme Activity
2.	265.	241.	195.
4.	521.	487.	505.
6.	662.	805.	754.
8.	885.	901.	898.
10.	884.	850.
12.	852.		914.
14.	932.	1110.	851.
16.	987.	954.	999.
18.	984.	961.	1105.
20.	954.	1021.	987.

When analyzing the data, use the Michaelis-Menton model. The results of the Prism analysis are as follows:

	Enzyme Activity
Michaelis-Menten
Best-fit values
Vmax	1353
Km	5.886
Std. Error
Vmax	75.93
Km	0.9498
95% Confidence Intervals
Vmax	1197 to 1509
Km	3.933 to 7.839
Goodness of Fit
Degrees of Freedom	26
R square	0.9041
Absolute Sum of Squares	170343
Sy.x	80.94
Constraints
Km	Km > 0.0
Number of points
Analyzed	28

Notice, again that the R² value is close to 1, and that when you look at the generated graph that the actual values deviate very little from the expected. To show you that this model fits the data best, let’s see what happens if we had tried to fit a linear regression model to the data.

Best-fit values
Slope	35.79 ± 4.456
Y-intercept when X=0.0	408.6 ± 55.70
X-intercept when Y=0.0	-11.42
1/slope	0.02794
95% Confidence Intervals
Slope	26.63 to 44.95
Y-intercept when X=0.0	294.1 to 523.1
X-intercept when Y=0.0	-19.32 to -6.649
Goodness of Fit
R square	0.7128
Sy.x	140.1
Is slope significantly non-zero?
F	64.53
DFn, DFd	1.000, 26.00
P value	< 0.0001
Deviation from zero?	Significant
Data
Number of X values	10
Maximum number of Y replicates	3
Total number of values	28
Number of missing values	2

As you can see, the R² value significantly decreases, and when you compare the linear to the nonlinear goodness of fit line, you see that there is much more deviation from the actual observed values from those that are predicted. That is why fitting the data to the model is so important, because if we use a model with poor fitness we are unlikely to make accurate predictions about the dependent variables from independent variables.

A more thorough walk through, as well as model fitting examples using the other analyzes provided by Prism, click here.