Read in the data (N.B. the first 20 rows are placebo, the next 20 are HRT) and have a look at the first 5 rows:

>  hrt  <- read .csv(file = " . ./ . ./data/hrt .txt")

>  hrt$treatment  <- c(rep(0, 20), rep(1, 20))

>  hrt[1:5,  ]

Question 1

Using the spaghetti plot function from the lecture (but modifying the apply call so that NAs are ignored):

>  spagplot  <- function(x,y,ylimit,xlabel,ylabel,heading) { + meany <- apply(y,1,mean, na .rm = TRUE)

+     matplot(x,y,type="l",lty=1,xlab=xlabel,

+                     ylab=ylabel,ylim=ylimit,xaxt="n",

+                     col="lightblue")

+      axis(side=1,at=x)

+     matpoints(x,y,type="p",pch=1,col="lightblue")

+      title(main=heading)

+      lines(x,meany,type="l",lwd=2,col="red")

+  }

> par(mfrow  =  c(1,2))

>  x  <- - 1:4

> pbo  <- t(subset(hrt, treatment == 0)[,1:6])

>  spagplot(x, pbo,  range(hrt[,1:6],  na .rm  =  TRUE),  "Time  (months)", +                    "Hamilton  Depression  Score",  "Placebo")

>  trt  <- t(subset(hrt, treatment == 1)[,1:6])

>  spagplot(x,  trt,  range(hrt[,1:6],  na .rm  =  TRUE),  "Time  (months)", +                    "Hamilton  Depression  Score",  "HRT")

Placebo

HRT

−1    0     1     2     3     4

Time (months)

−1    0     1     2     3     4

Time (months)

Question 2

Again using the function from the lecture:

> meanplot  <- function(x,means,xlabel,ylabel,heading, + legnames,legloc) {

+     matplot(x,means,type="l",lty=1,xlab=xlabel,ylab=ylabel, +                     xaxt="n",col=c("blue","red"))

+      axis(side=1,at=x)

+     matpoints(x,means,type="p",pch=1:2,col=c("blue","red")) +      title(main=heading)

+      legend(legloc,  legnames,  lty=1,  col=c("blue","red"),

+                   pch=1:2)

+  }

> means1  <- cbind(apply(pbo, 1, mean, na .rm = TRUE),

+                                 apply(trt,  1, mean,  na .rm  =  TRUE))

> par(mfrow  =  c(1,1))

> meanplot(x, means1,  "Time  (months)",  "Hamilton  Depression  Score",

+

+

"Mean profiles  for HRT data",

c("placebo","HRT"),"topleft")

Mean profiles for HRT data

−1                0                 1                 2                 3                 4

Time (months)

Question 3

The above plots indicate a tendency for the Hamilton Depression Score (HDS) to increase (i.e., worsen) as time progresses from baseline with the possibility that HRT slows the rate of increase somewhat: statistical tests of hypotheses are required.

The spaghetti plots suggest that the variance increases with increasing HDS and/or with the progression of time. They also suggest that there is a need to consider random effects at baseline in any model for the time courses because many traces tend to stay on a higher (or lower) trajectory for their whole time course. One possible model for consideration later would be one that includes a random intercept and a random trend component.

Question 4

Calculate the changes from the first baseline measurement:

>  hrt_diff  <- cbind(hrt[,2:6] - hrt[,1], treatment = hrt$treatment) >  hrt_diff[1:5,]

Spaghetti plots:

> par(mfrow  =  c(1,2))

>  x_diff  <- x[- 1]

> pbo_diff  <- t(subset(hrt_diff, treatment == 0)[,1:5]) >  spagplot(x_diff, pbo_diff,  range(hrt_diff[,1:5],  na .rm  =  TRUE),

+

+

"Time  (months)",

"Change  in  HDS  from  1st measurement",  "Placebo")

>  trt_diff  <- t(subset(hrt_diff, treatment == 1)[,1:5]) >  spagplot(x_diff,  trt_diff,  range(hrt_diff[,1:5],  na .rm  =  TRUE),

+

+

"Time  (months)",

"Change  in  HDS  from  1st measurement",  "HRT")

Placebo

HRT

0       1       2       3      4

0       1       2       3      4

Time (months)                                          Time (months)

Mean profile plots:

> means1_diff  <- cbind(apply(pbo_diff, 1, mean, na .rm = TRUE), + apply(trt_diff, 1, mean, na .rm = TRUE)) > par(mfrow  =  c(1,1))

> meanplot(x_diff, means1_diff,  "Time  (months)",

+

+

+

"Change  in  HDS  from  1st measurement",

"Mean  change profiles  for HRT data",

c("placebo","HRT"),"topleft")

Mean change profiles for HRT data

3

Time (months)

The above graphs showing changes from baseline confirm the observations made above for the complete time course including baseline measurements.

Question 5

Placebo:

>  cov_pbo  <- var(t(pbo), use = "pairwise .complete .obs") > print(round(cov_pbo,  2))

base1  base2  month1 month2 month3 month4

base1 8 .01 6 .22 6 .36 5 .58 7 .25 7 .75

base2 6 .22 6 .45 6 .44 5 .42 7 .78 8 .52

month1 6 .36 6 .44 7 .54 5 .26 8 .25 8 .58

month2 month3 month4