Parameter uncertainty

Just how confident are you in your parameter values, anyway?

You measure an observable over time and would like to extract a rate constant $k$ for the process. You are fitting the data using the code below. Each data point was measured independently of the others. The uncertainty in each measured point is normally distributed with a standard deviation of 0.05.

How robust are your estimates of $A$ and $k$ given the uncertainty in your measured points?

• Generate 1,000 simulated data sets where each experimental point is drawn from a normal distribution with a mean of d.obs and a standard deviation of 0.05.

• Generate a histogram of possible values of $A$ and $k$ from these simulations. (Hint: try fitting each simulated dataset independently...).

• What are the 95% confidence intervals on your estimate of $k$? The lower bound is the value of $k$ for which 2.5% of the histogram counts are below the value. The upper bound is the value of $k$ for which 2.5% of the histogram counts are above the value.

• You measure the same process under slightly different conditions. These data are stored in data/time_course_1.csv. Is there a statistically significant difference between $k$ from dataset 1 vs. 0?

Bonus:

These are a couple of challenge questions for students who are already comfortable fitting and want to extend this basic fitting example.

• You repeat the experiment again and explicitly measure the variance of each point (meaning some are now 0.1, others 0.3, etc.). These values are stored in data/time_course_2.csv in the obs_err column. Modify your sampling code so you incorporate these point-specific uncertainties.

• "Global fitting" is a powerful way to squeeze the most information out of experiments. Rather than fitting models to individual experimental replicates and then estimating parameter uncertainty, one can fit all experiments at once to extract a single estimate of the fit parameters. Implement a global fit function that simultaneously fits the data in data/time_course_2.csv, data/time_course_3.csv, and data/time_course_4.csv, then report your estimate and 95% confidence intervals of $k$. Do your confidence intervals change with more data?

Model selection

How do you choose which model to use?

The Akaike Information Criterion (AIC) helps you select between competing models. It penalizes models for the number of fittable parameters because adding parameters will almost always make a fit better.

$$AIC = -2ln(\hat{L}) + 2k$$

where $\hat{L}$ is the "maximum likelihood" and $k$ is (1 + number of fittable parameters in the model). $\hat{L}$ is proportional to the sum of the squared residuals ($SSR$), so we can write:

$$ AIC = nln(SSR) + 2k $$

where $n$ is the number of observations. If we want to compare $N$ different models, we first normalize $AIC$ for each model as $\Delta _{i}$.

$$\Delta _{i} = AIC_{i} - AIC_{min}$$

where $AIC_{min}$ is the best AIC. The support for a given model is given by the Akaike weight ($w$):

$$w_{i} = \frac{exp(-\Delta_{i}/2)}{\sum_{j=0}^{j<N} exp(-\Delta_{j}/2)}$$

The weight runs from 0.0 to 1.0, with 1.0 representing the strongest support. For a more thorough description, see here.

• Create a function called calc_aic that returns the Akaike weight for a list of models. calc_aic should take the following arguments:
  • list of ssr for each fit
  • list of num_parameters for each fit
  • number of observations you fit the models to

You can check your function with the following lists.

Real example: what can you learn from your data?

Sedimentation equilibrium experiments are used to study molecular assemblies. If spin a solution of molecules in a centrifuge very fast for a very long time, you get an equlibirum distribution of molecules along the length ($r$) of the tube. This is because centrigual force pulls them to the bottom of the tube, while diffusion tends to spread them out over the tube. If you measure the concentration of molecules along the length of the tube $c(r)$, you can then fit a model and back out the molecular weight(s) of the species in solution.

The equation describing a simple, non-interacting molecule is below. $c_{0}$ is the concentration of the molecule at $r=0$, $M$ is the molecular weight of the molecule, and $r$ is the position along the tube.

$$c(r) = c_{0}exp \Big [ M \Big ( \frac{r^{2}}{2} \Big ) \Big ]$$

You are interested in whether there is any dimer present (a dimer is a pair of interacting molecules). If there is dimer, our model gets more complicated. We have to add a term $\theta$ that describes the fraction of the molecules in dimer versus monomer:

$$c(r) = c_{0} \Big ( \theta exp \Big [ M \Big ( \frac{r^{2}}{2} \Big ) \Big ] + (1 - \theta )exp \Big [ 2M \Big ( \frac{r^{2}}{2} \Big ) \Big ] \Big )$$

(I've collapsed in a bunch of constants to simplify the problem. If you want a full discussion of the method, see here, particularly equation 4).

Main question: do these data provide evidence for a dimer?

• Implement a function and residuals function for each of these models. The first model should have the fittable parameters $c_{0}$ and $M$. The second model should have the fittable parameters $c_{0}$, $M$, and $\theta$.

• Fit both models to the data in data/sed_eq.csv. What are your estimates of $c_{0}$, $M$, and $\theta$? Are they the same between the two fits?

• Use your calc_aic function on these fits. Which model is supported? Can you conclude there is dimer present?

Gaussian

One common type of data is a sum of Gaussian functions. Among other places, these sort of data pop up in single-molecule work, fluorescence-activated cell sorting, spectroscopy and chromatography. You collect a dataset that looks like the following:

You find code to analyze this kind of data on the internet. Using the functions below, determine:

• how many gaussians you can extract from the data in lab-data/gaussian.csv.
• their means, standard deviations, and areas.

PS: you should play with your initial guesses PPS: negative area gaussians are not allowed