Study on risks factor for low birth weights from Baystate medical center in Springfield (MA, USA) during 1986 [Hosmer & Lemeshow (1989), Applied Logistic Regression]

- low: low birth weight indicator coded as: 1 (low birth weight, i.e. \(\leq 2.5kg\)) and 0 (normal birth weight)
- age: mother’s age (in years) at the time of birth
- lwt: mother’s weight (in pounds) at the time of the last menstrual period
- race: mother’s race coded as: 1 (white), 2 (black) and 3 (other)
- smoke: mother’s smoking status coded as 1 (smoker) and 0 (non-smoker)
- ptl: number of previous premature labours
- ht: hypertension history indicator coded as: 1(yes) and 0 (no)
- ui: uterine irritability indicator coded as: 1 (yes) and 0 (no)
- ftv: number of physician visits during the first trimester
- bwt: infant’s birth weight (in grams)

1. launch Rstudio
2. open a new .Rmd file from the default template
3. compile it with the "Knit" button
4. discuss with your neighbor to understand each part of the script

1. import the dataset in birthweight.txt (you can use the "Import Dataset" button from Rstudio
2. describe the data (use nice table outputs in Rmarkdown and ggplot2 graphics)

• program a function that add the logarithm of two numbers
• use browser() to debug your function when evaluating at negative arguments

Statistics: summarizing information from experimental observations and quantifying the associated uncertainty

Always start with the research/scientific question !

• ideally we would start from our scientific question, figure out what data we need to answer the question
• usually the question is not very clear & the data have already been collected

Statistical Inference: we use a simple Generative Probabilistic model that could have generated the observations (Machine Learning sometimes reject this paradigm – cf. L. Breiman)

The likelihood is a fundamental building block of Biostatistics:

• its maximization yields estimators that have good asymptocic (when \(n\) the number of observations increases indefinitely) properties – no bias and smallest variance,
• it also allows to perform significance testing of model parameters in many settings.

The likelihood function quantifies how likely it is that a given (set of) observation(s) has been generated by our hypothesized Generative Probabilistic model.

The likelihood function is equal to the model joint probability distribution computed for the observations, and thus only a function of the model parameters.

The idea of the MLE is to to optimize the likelihood function given the observations, by finding the model parameters that would give these observations the highest probability of being generated under the model.

seems like a reasonable and intuitive idea !

Reverend Thomas Bayes proposed an alternative framework for statistical inference (actually before the “frequentist” method). It also relies on a probabilistic model through the likelihood function, but has different philosophical grounds than the frequentist.

To be continued…

Maximizing the likelihood can easily be done analytically for simple linear models.

However: non-linear likelihoods are hard (sometimes impossible) to optimize analytically !

\(\Rightarrow\) numerical optimization

An algorithm to find values for which a function is zero.

Applied to the derivative of the likelihood function, this will identify the MLE
given that the log-likelihood is a concave function

Generally, we maximize the log-likelihood instead of the likelihood.

• same thing (\(\log\) is a monotonically increasing function)
• products \(\prod_{i=1}^n\) (e.g. likelihoods) become sums \(\sum_{i=1}^n\)
  sums are easier to compute & easier to derivate
• more precision for small positives numbers (e.g. small probabilities)

This is not taught often but it can come very handy if you are writing your own statistical/optimisation program:

\[ \log \sum_{i=1}^n e^{x_i} = c + \log \sum_{i=1}^n e^{x_i-c} \]

Disclaimer: not useful today

1. start from an initial point \(x_0\)
2. tangential linear approximation of \(f\) in \(f(x_0)\)
3. get the next point \(x_1\) as the solution
   …

1. Repeat until \(x_{n+1}\approx x_n\)

About the linear approximation at the \(n + 1\) step: - goes through \(f(x_n)\) - has slope \(f'(x_n)\)

So it has the following equation: \(y= f'(x_n)(x-x_n) + f(x_n)\). Thus we find \(x_{n+1}\) by setting \(y=0\), which gives us \(x_{n+1} = x_n -\frac{f(x_n)}{f'(x_n)}\)

#install.packages("animation")
library("animation")
newton.method(FUN = function(x) (x - 2)^2 - 1, init = 9.5, 
              rg = c(-1, 10), tol = 0.001, interact = FALSE, 
              col.lp = c("orange", "red3", "dodgerblue1"), 
              lwd=1.5)

1. propose a generative probabilistic model

2. define the parameter of interest

3. program the associated likelihood function

4. maximize this likelihood analytically

Let’s Program a Newton-Raphson algorithm to maximise this likelihood numerically

1. write two funtions that computes the first and second derivatives of the log-likelihoodrespectively

2. write a Newton-Raphson function with 5 arguments (the first deriative of the function to maximize, its second derivative, the initial starting point, the tolerance, the maximum number of iterations)

3. use all three functions to compute the MLE of the low birthweight prevalence

1. Try several initial values, different tolerance, and experiment with different stopping rules – i.e. convergence criterion – (you can combine them). Compare the results. Does it always work ?
2. Compute the p-value from the Wald test, the score test and the Likelihood Ratio Test respectively, and compare to the glm function output

1. Repeat the exercise to assess wether the mother’s smoking status impacts the probability of low birth weight