The beta-binomial distribution is a combination of the binomial and the beta distribution. For a test with n items, which are scored dichotomously with 0 (= false) and 1 (= correct), let X be the number of correct responses. The probability function of the beta-binomial distribution is given by: \[P(X = k | n, \alpha, \beta) = \binom{n}{k} \frac{B(k + \alpha, n - k + \beta)}{B(\alpha, \beta)}\] where \(B(\cdot,\cdot)\) is the beta function, and \(\alpha\) and \(\beta\) are shape parameters of the beta distribution.

We model the shape parameters \(\alpha\) and \(\beta\) of the beta-binomial distribution as functions of age (or another explanatory variable) using polynomial regression. Specifically: \[\log(\alpha(age)) = a_0 + a_1age + a_2age^2 + ... + a_mage^m\] \[\log(\beta(age)) = b_0 + b_1age + b_2age^2 + ... + b_sage^s\] where m and s are the degrees of the polynomials for \(\alpha\) and \(\beta\) respectively. We use the logarithm of \(\alpha\) and \(\beta\) to ensure that they are positive, since the beta-binomial distribution is only defined for positive \(\alpha\) and \(\beta\). This transformation also helps in stabilizing the variance and improving the optimization process. The mean \(\mu\) and variance \(\sigma^2\) of the beta-binomial distribution can be derived from \(\alpha\) and \(\beta\) as follows: \[\mu = \frac{n\alpha}{\alpha + \beta}\] \[\sigma^2 = \frac{n\alpha\beta(\alpha + \beta + n)}{(\alpha + \beta)^2(\alpha + \beta + 1)}\] To estimate the parameters (\(a_0, ..., a_m, b_0, ..., b_s\)), we use maximum likelihood estimation. The log-likelihood function for N observations is:

\[L(a, b | X, Age) = \sum_{i=1}^N \log[P(X_i | n, \alpha(Age_i), \beta(Age_i))]\]

where \(X_i\) is the score and \(Age_i\) is the age for the i-th observation.

The data fitting is performed using numerical optimization techniques, specifically the L-BFGS-B (Limited-memory BFGS) algorithm of the 'optim' function. The latter is a quasi-Newton method for solving large nonlinear optimization problems with simple bounds. By approximating the Hessian matrix, it simultaneously determines the coefficients of the regression equations for \(\alpha\) and \(\beta\) that maximize the log-likelihood, thus providing the best fit to the observed data.