1.3 Motivating example

The first data set contains the birthweight and gestational age (time from conception to birth) of 24 babies born in hospital. The data is held in R as birthweight. A scatter plot of this data set is presented below.

birthweight = read.table(’birthweight.dat’)
attach(birthweight)
plot(age, weight, pch=16)
abline(a = -1485.0, b = 115.5)

These observations suggest a linear trend of birthweight with increasing gestational age, together with some random scatter around the overall trend. A simple linear model is a clear candidate for modelling this data. The fitted least-squares line is show on the plot.

Let $Y_i$ denote the birthweight, for $i=1,2,\mathrm{\dots },n$, and $x_i$ denotes the gestational age for baby $i$. The linear model is:

$Y_i={\mu }_i+{ϵ}_i\mathit{\hspace{1em}}\mathrm{with}\mathit{\hspace{1em}}{\mu }_i=a+b{x}_i,\mathrm{and}\mathit{\hspace{1em}}{ϵ}_i\sim N(0,{\sigma }^2)$

This can be equivalently be represented as,

$Y_i\sim N({\mu }_i,{\sigma }^2)\mathit{\hspace{1em}}\mathrm{with}\mathit{\hspace{1em}}{\mu }_i=a+b{x}_i.$

Using the lm command in R we can easily find the fitted mean line as $\widehat{\mu }=-1485.0+115.5x$.

 
Exercise 1.3

1. 1.

   Interpret this finding in context.

2. 2.

   Estimate a suitable value of $\sigma$.

 

Between 1983 to 1986, data was collected on the number of deaths from AIDS in Australia in consecutive three-month periods. The data is held in R as aids and shown below.

The dotted line on the plot is the least-squares regression line estimated from the data. At first sight, this linear model does not appear to fit too badly, but we notice that the fitted values are negative for periods 1 and 2. This is a bad aspect of the model since we know that negative values are not possible. Also, the observations are counts, not continuous, and the variance seems to increase with the mean.

A reasonable model for this data might be that the number of deaths $Y_i$ at each time $t_i$ is Poisson(${\mu }_i$), with the means ${\mu }_i$ increasing in time. To make sure that these means do not go negative, we can model ${\mu }_i=\exp \{a+b{t}_i\}$ which is always positive. This model can be written as:

$Y_i\sim \mathrm{Pois}({\mu }_i)\mathit{\hspace{1em}}\mathrm{with}\mathit{\hspace{1em}}\log ({\mu }_i)=a+b{t}_i.$

Fitting this model to the data by maximum likelihood, we obtain estimates $\widehat{a}=0.340$ and $\widehat{b}=0.257$. The corresponding fitted curve, $\widehat{\mu }=\exp \{\widehat{a}+\widehat{b}t\}$ is shown in the plot as the solid line, which provide a better description of the data.

aids = read.table(’aids.dat’)
attach(aids)
plot(time, number)
mu = exp(0.340 + 0.257*time)
line(time, mu, col="blue")

 
Exercise 1.4
Give an interpretation of the model

$Y_i\sim \mathrm{Pois}({\mu }_i)\mathit{\hspace{1em}}\mathrm{with}\mathit{\hspace{1em}}\log ({\mu }_i)=0.340+0.257t_i.$

 

This model is not a causal model, but merely a description of how the AIDS epidemic is growing.

The Poisson model for the AIDS data can be defined in general terms as:

	$Y_i\sim G({\mu }_i)\mathit{\hspace{1em}}\mathrm{where}\mathit{\hspace{1em}}𝔼[Y_i]={\mu }_i,$
	$g({\mu }_i)={\eta }_i\mathit{\hspace{1em}}\mathrm{and}{\eta }_i={\beta }^{\prime }{𝐱}_i$

where we have to identify:

• •

  $i$ index, here 3 month period.

• •

  $Y_i$ response, here number of AIDS deaths.

• •

  ${𝐱}_i$ covariates, here time period.

• •

  $\beta$ coefficients, here $\beta =(a,b)$.

• •

  $G$ distribution, here Poisson.

• •

  $g$ link function, here $\log$.

In a phase I clinical trial to find the effective dose of a new drug where patients were randomly assigned to receive different doses of the drug. The table below shows the number $z_i$ of patients responding positively to the drug for each dose $x_i$.

This data is held in R as clintrial. Below is a plot of the observed proportion of positive responses, $y_i=z_i/m_i$, at each dose.

The fitted least-squares line of best fit is shown as a dashed-line. This goes outside the defined $(0,1)$ interval for proportions. We therefore need to take account of the non-normal distribution of the data.

For each dose, a binomial model is appropriate, for example:

$Z_i\sim \mathrm{Binomial}(m_i,{\mu }_i)$

where ${\mu }_i$ is the probability that a patient will respond at does $x_i$. Here, the expectation is $𝔼[Z_i]=m_i{\mu }_i$. To match with the previous examples, we consider the the transformation $Y_i=Z_i/m_i$ to define the proportion of positive responses for each dose level. The expectation of the proportion is therefore $𝔼(Y_i)=𝔼(Z_i)/m_i={\mu }_i$. Taking the proportion as our response variable, it follows that $Y_i$ is a binomial-proportion model:

$Y_i\sim \mathrm{Binoprop}(m_i,{\mu }_i).$

Now we need to model ${\mu }_i$ as a function of $x_i$. Assuming a linear relationship such as ${\mu }_i=a+b{x}_i$ may conflict with the interpretation of $\mu$ as a probability. We must then constrain the mean to the required interval. One such mapping is the logistic function:

$h(x)={\displaystyle \frac{\exp \{x\}}{1+\exp \{x\}}}$

A proposed model for the clinical trial data could then be:

$Y_i\sim \mathrm{Bionprop}(m_i,{\mu }_i),\mathrm{with}\mathit{\hspace{1em}}\mathrm{logit}({\mu }_i)=a+b{x}_i$

The logit function is the inverse function to the logistic function.

 
Exercise 1.5
Evaluate the logit function.

 

The solid line in the plot above presents the curve of best fit with maximum likelihood estimates $\widehat{a}=-60.1$ and $\widehat{b}=33.9$.

clintrial = read.table(’clintrial.dat’)
attach(clintrial)
plot(dose, propn)
eta = -60.1 + 33.9 * dose
mu = exp(eta)/(1+exp(eta))
lines(dose, mu, col="blue")

 
Exercise 1.6
Give a contextual interpretation to this model, using the plot.

 

 
Exercise 1.7
The Binoprop model follows the general structure:

	$Y_i$	$\sim G({\mu }_i)\mathit{\hspace{1em}}\mathrm{where}\mathit{\hspace{1em}}𝔼(Y_i)={\mu }_i,$
	$g({\mu }_i)$	$={\eta }_i\mathit{\hspace{1em}}\mathrm{and}{\eta }_i={\beta }^{\prime }{x}_i.$

where each component is:

 

Twenty-six dishes of bacterial cultures were given different amounts, $x_i$, of a drug. For each, it was recorded whether they responded to the drug, $Y_i=1$, or not, $Y_i=0$. This data is held in R as bacteria. Below is a plot of the bacteria data at each dose level.

It is particularly apparent in this case that a normal linear model is not appropriate. The line of best fit again extends beyond the set of valid probability values. The response is categorical not continuous, as it can only take one of two values. In no sense can we say that the residuals from the fit are approximately normal distributed about the fitted line!

A more sensible model might assume that at each dose, $x_i$, there is a probability ${\mu }_i$ that a bacterial culture dish will respond so that $Y_i\sim \mathrm{Bernoulli}({\mu }_i)$. This is a special case of the binomial model discussed in the previous section. It makes sense to model ${\mu }_i$ in the same way:

	$Y_i$	$\sim \mathrm{Bernoulli}({\mu }_i)\mathit{\hspace{1em}}\mathrm{with}$
	$\mathrm{logit}({\mu }_i)$	$=a+b{x}_i.$

The curve of best fit using the maximum likelihood estimates $\widehat{a}=-4.111$ and $\widehat{b}=3.581$ is represented in the above plot with a solid line.

bacteria = read.table(’bacteria.dat’)
attach(bacteria)
plot(dose, response)
eta = -4.111 + 3.581 * dose
mu = exp(eta)/(1+exp(eta))
lines(dose, mu, col="blue")