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 its 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)