MATH5885 Longitudinal Data Analysis Take-home exercises
Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit
MATH5885 Longitudinal Data Analysis
Week 1, Lecture 1
Take-home exercises
Read in the dataset and extract the placebo group:
> vnames <- c("ID", paste0("Hour", 0:8), "treatment")
> resptrial <- read .table(" . ./ . ./data/resptrial .txt", + header = FALSE, col .names = vnames) > resptrial_placebo <- subset(resptrial, treatment == "p")
Have a look at the first 5 rows:
> resptrial_placebo[1:5, ]
ID Hour0 Hour1 Hour2 Hour3 Hour4 Hour5 Hour6 Hour7 Hour8 treatment
49 201 |
2 .14 |
2 .36 |
2 .36 |
2 .28 |
2 .35 |
2 .31 |
2 .62 |
2 .12 |
2 .42 |
p |
50 202 |
3 .37 |
3 .03 |
3 .02 |
3 .19 |
2 .98 |
3 .01 |
2 .75 |
2 .70 |
2 .84 |
p |
51 203 |
1 .88 |
1 .99 |
1 .62 |
1 .65 |
1 .68 |
1 .65 |
1 .85 |
1 .96 |
1 .30 |
p |
52 204 |
3 .10 |
3 .24 |
3 .37 |
3 .54 |
3 .31 |
2 .81 |
3 .58 |
3 .76 |
3 .05 |
p |
53 205 |
2 .91 |
3 .35 |
3 .92 |
3 .69 |
3 .97 |
3 .94 |
3 .63 |
2 .92 |
3 .31 |
p |
Create the pairs scatterplot. N.B. We have to first define the panel .cor function:
> panel .cor <- function(x, y, digits=2, prefix="", cex .cor) { + usr <- par("usr"); on .exit(par(usr))
+ par(usr = c(0, 1, 0, 1))
+ r <- abs(cor(x, y))
+ txt <- format(c(r, 0 . 123456789), digits=digits)[1]
+ txt <- paste(prefix, txt, sep="")
+ if(missing(cex .cor)) cex <- 0 .8/strwidth(txt)
+ text(0 .5, 0 .5, txt, cex = cex * r)
+ }
> pairs(resptrial_placebo[,2:10], upper .panel = panel .cor)
1.5 3.5
1.5 3.5
1.0 3.0
1.0 3.0
Maybe restrict it to the first few so that it’s a bit easier to see:
> pairs(resptrial_placebo[,2:5], upper .panel = panel .cor)
Hour0 |
|
|
|
|
1.5 2.5 3.5
0.82 |
Hour1 |
|
|
0.70 |
0.9 |
0.96 |
0.94 |
|
Time plot for all placebo patients:
> x <- 0:8
> y <- t(resptrial_placebo[,2:10])
> meany <- apply(y, 1, mean)
> matplot(x, y, type = "l", lty = 1, xlab = "Time (hours)", + ylab = "FEV1", xaxt = "n")
> axis(side = 1, at = x)
> matpoints(x, y, type = "p", pch = 1)
> title(main = "FEV1 profiles for respiratory trial placebo group") > lines(x, meany, type = "l", lwd = 3)
FEV1 profiles for respiratory trial placebo group
0 1 2 3 4 5 6 7 8
Time (hours)
Or just 5 randomly selected patients:
> set .seed(100)
> inds <- sample(1:nrow(resptrial_placebo), 5)
> y1 <- y[, inds]
> matplot(x, y1, type = "l", lty = 1, xlab = "Time (hours)", + ylab = "FEV1", xaxt = "n")
> axis(side = 1, at = x)
> matpoints(x, y1, type = "p", pch = 1)
> title(main = "FEV1 profiles for respiratory trial selected patients") > text(rep(0, 5), y1[1,] + 0 . 1, resptrial_placebo$ID[inds], col = 1:5)
FEV1 profiles for respiratory trial selected patients
224
220
210
206
217
|
|||||||||
|
|
|
|
|
|
|
|
|
|
0 1 2 3 4 5 6 7 8
Time (hours)
2022-08-06