Lab sessions built upon practicals provided by Steve Coad, SMS (QMUL).

Practical 1 - 5 October 2022 (Week 2)

This practical reminds you about the language and commands of R which you will practice through the statistical computing package RStudio.  This package should be available if you are using a university computer, or you may use your own laptop with your installation of the package. Another option to run your R commands is to use the site rdrr .io/snippets/.

Open RStudio and from the menu File>New  File>R  Script, open a new script in which you will type commands (alternatively, use Control+Shift+N). Remember to save this file for your records and for later practice.

The data

In this practical, we will refresh some of the techniques that you used in Statistical Modelling I, by fitting a familiar simple linear regression model.

Manatees are large, gentle sea creatures that live along the Florida coast. Many manatees are killed or injured by powerboats.  Below are data on powerboat regis- trations (x), in thousands, and the number of manatees killed by boats (y) in Florida in the years 1977 to 1987.

Year    1977    1978    1979    1980    1981    1982    1983    1984    1985    1986    1987

Do the data provide any evidence to conclude that there is a linear relationship between the number of powerboat registrations and the number of manatees killed by boats?

Entering and plotting the Data

First enter the x and y values as vectors in R:

x  <-  c(447,460,481,498,513,512,526,559,585,614,645)

y  <-  c(13,21,24,16,24,20,15,34,33,33,39)

Then produce a scatterplot of the data using

plot(x,y,main="Plot  of  y  against  x")

Does the relationship between y and x seem to be linear?

Fitting the Model

Fit a simple linear regression model to the data by

manatee  <-  lm(y  ~  x)

To see the details of the fitted model, we use

summary(manatee)

The values of the intercept βˆ0  and the slope βˆ1  can be seen from the output.  Add the fitted line to your scatterplot by using the following commands:

plot(x,y,main="Fitted  Line  Plot")

abline(manatee)

To test H0  : β1  = 0, we look at the analysis of variance table:

anova(manatee)

What are your conclusions?

Checking the Assumptions

To check the assumptions of the model, we examine the residual plots. We first plot the standardised residuals against x:

stres  <-  rstandard(manatee)

plot(x,stres,main="Standardised  Residuals  against  x")

Is there any reason to doubt the linearity of the model?   Next, we plot the standardised residuals against the fitted values:

fits  <-  fitted(manatee)

plot(fits,stres,main="Standardised  Residuals  against  Fitted Values")

Is there any reason to doubt that the variance is constant?  Finally, we look at the Q-Q plot:

qqnorm(stres,main="Q-Q  Plot")

qqline(stres)

Is there any reason to doubt the assumption that the errors are normally dis- tributed?

Finishing the Session

Remember to save your files and log out.

Practical 2 - 12th October 2022 (Week 3)

This practical reminds you how to read data from a file into R and how to fit a multiple linear regression model.

Diabetic Data

The file diabetic .csv on the module webpage contains data for 20 male insulin- dependent diabetics who had been on a high carbohydrate diet for six months. Compliance with the regime is thought to be related to age (x1 ), in years, body weight as a percentage of ‘ideal’ weight for height (x2 ) and the percentage of calories as protein (x3 ). The dependent variable is the percentage of total calories obtained from complex carbohydrates.

Load the data into R as follows.  First, you have to tell R where you have saved the data. This is known as your working directory. You set this by telling R where the data are.  For example, it is probably convenient to save the data onto the G:

drive, in which case you use

setwd("G:")

You can check that it is correct by

getwd()

Now read in the data by

diabetic  <-  read .csv("diabetic .csv")

The data will have been read into R, but they are stored in four columns, and we need to allocate the columns of the matrix to y , x1 , x2  and x3 . This is achieved by

y  <-  diabetic[,1]

x1  <-  diabetic[,2]

x2  <-  diabetic[,3]

x3  <-  diabetic[,4]

Check that data has been correctly read by looking at the first few rows using head(diabetic)

and examine the collection of pairwise scatterplots of data with

pairs(diabetic)

Model 1

Fit a multiple linear regression model to the data in which y is linearly related to

each of the explanatory variables by

diabetic1  <-  lm(y  ~  x1  +  x2  +  x3)

Obtain the fitted model using summary. Save the fitted values and the standardised residuals. Assess the assumptions of normality and constant variance of the random errors by examining suitable residual plots. What are your conclusions?

Model 2

Fit a multiple linear regression model to the data in which y is linearly related to x2  and x3  by

diabetic2  <-  lm(y  ~  x2  +  x3)

Obtain the fitted model and save values as before.  Assess the assumptions of nor- mality and constant variance of the random errors.  Compare the two model fits. Which one is best and why?

Model 3

From looking at the initial results (Model 1), it is clear that x3 is the most important term. Fit a univariate regression model explaining y in terms of x3 only and compare it with the previous two models.

Finishing the Session

Remember to save your files and log out.

Practical 3 - 19 October 2022 (Week 4)

This practical introduces you to generalised linear models by modelling some binary response probabilities.

Beetle Data

In bioassays, the response may vary with a covariate x termed the dose.  Data for a typical example involving a binary response are given in the file beetle .csv on the module webpage.  Here, a certain number of beetles (r) are exposed to various concentrations of gaseous carbon disulphide, in milligrammes per litre, for five hours and the number of beetles killed (y) is recorded. The dose is the base 10 logarithm of the concentration.

The three columns of beetle .csv contain x, r and y in that order.   Create variables x, r and y from the columns of the data.  r and y.  The proportion of

beetles killed at each of the doses may be calculated using

p  <-  y/r

Now plot these proportions against the doses by

plot(x,p,main="Plot  of p  against  x")

The purpose of this practical is to identify a model which provides a good description of the data.

Logistic Model

Fit a logistic model to the data in which the probability of a beetle being killed πi at the ith dose xi  satisfies

log ( ) = β0 + β1 xi

by

beetle1  <-  glm(p  ~  x,family=binomial(link=logit),weights=r)

Note that, since a logit link is the default one for the binomial distribution, it is

sufficient to put family=binomial as the second argument in this case.

Obtain the fitted model using summary. The standard errors of the intercept βˆ0 and the slope βˆ1  are the square roots of the diagonal elements of the inverse of the estimated Fisher information matrix.  Under the null hypothesis that the logistic model describes the data, the residual deviance has an approximate χ6(2)  distribution, since there are eight doses and two parameters in the model. Does this model provide a good fit?

Probit Model

Fit a probit model to the data in which Ti  satisfies

Φ − 1 (Ti ) = β0 + β1 xi ,

where Φ denotes the standard normal distribution function, by              beetle2  <-  glm(p  ~  x,family=binomial(link=probit),weights=r)

Obtain the fitted model as before. Compare the two model fits. Are they similar?

Extreme Value Model

Fit an extreme value model to the data in which Ti  satisfies

log{−log(1 − Ti )} = β0 + β1 xi

by

beetle3  <-  glm(p  ~  x,family=binomial(link=cloglog),weights=r)

Here, cloglog means that a complementary log-log link is being used.  Obtain the fitted model as before. Which of the three models provides the best description of the data?

Finishing the Session

Remember to save your files and log out.

Practical 4 - 26 October 2022 (Week 5)

This practical shows you how to fit a Poisson regression model to count data.

Count Data

The data below are counts (y) observed at various values of a covariate (x).

y     2     3   6   7   8   9    10    12    15

x   -1   -1   0   0   0   0     1     1     1

Enter the x and y values as vectors in R. Call them x and y.  Then produce a scatterplot of the data. Is the variance constant?

Linear Model

Fit a linear model to the data in which the mean response µi for covariate xi satisfies µi  = β0 + β1 xi

by

count1  <-  glm(y  ~  x,family=poisson(link=identity))

Note that, since an identity link is not the default one for the Poisson distribution, it is necessary to put family=poisson(link=identity) as the second argument in this case.

Obtain the fitted model using summary. Under the null hypothesis that the linear model describes the data, the residual deviance has an approximate χ7(2)  distribution, since there are nine observations and two parameters in the model. Does this model provide a good fit?

Log-Linear Model

Fit a log-linear model to the data in which µi  satisfies

log(µi ) = β0 + β1 xi

by

count2  <-  glm(y  ~  x,family=poisson(link=log))

Obtain the fitted model as before.   Which of the two models provides the best description of the data?

Finishing the Session

Remember to save your files and log out.

Practical 5 - 2 November 2022 (Week 6)

This practical shows you how to assess whether a link function is appropriate and how to check the fit of a model.

Beetle Data Revisited

In Practical 3, you fitted three models to the beetle data. Each of these was based on a different link function.  The purpose of this practical is to see which of these link functions is best, and then to examine the fitted values of the chosen model.

Assessing a Link Function

To assess whether a logit link function might be appropriate, plot the empirical logits

log ( )

against the doses. Note that the halves are used to avoid problems when yi  = 0 or Ti . The plot is produced by

el  <-  log((y  +  0 .5)/(r  -  y  +  0 .5))

plot(x,el,main="Plot  of  Empirical  Logits  against  x")

and should be approximately linear.

The possibility of a probit link function can be examined by plotting the empirical

probits

Φ − 1  ( )

against the doses. This time, the plot is produced by

ep  <-  qnorm((y  +  0 .5)/(r  +  1))

plot(x,ep,main="Plot  of  Empirical  Probits  against  x")

Is this plot similar to the first one?

Similarly, to see whether the complementary log-log link function is adequate, plot the empirical log-log values

log { −log ( )}

against the doses. Here, the plot is produced by

ell  <-  log(-log((r  -  y  +  0 .5)/(r  +  1)))

plot(x,ell,main="Plot  of  Empirical  Log-Log Values  against  x") Which of the three link functions appears to be best?