MATH5885 Longitudinal Data Analysis Exercise 1
Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit
MATH5885 Longitudinal Data Analysis
Week 2
Exercise 1
Read in the data, creating a new variable for log(FEV1/Height), and have a look at the first 5 rows:
> fev <- read .table(" . ./ . ./data/fev1 .txt", header = TRUE)
> fev$lfev1h <- fev$lfev1 - log(fev$height)
> head(fev, n = 5)
subjectid height age iheight iage lfev1 lfev1h
1 |
1 |
1.20 9.3415 |
1.2 9.3415 0.21511 0.03278844 |
2 |
1 |
1.28 10.3929 |
1.2 9.3415 0.37156 0.12469992 |
3 |
1 |
1.33 11.4524 |
1.2 9.3415 0.48858 0.20340106 |
4 |
1 |
1.42 12.4600 |
1.2 9.3415 0.75142 0.40076313 |
5 |
1 |
1.48 13.4182 |
1.2 9.3415 0.83291 0.44086791 |
Question 1
Extract the residuals from the cubic spline fit:
> fev$resid <- resid(smooth .spline(fev$age, fev$lfev1h))
We can produce a spaghetti plot of the residuals:
> plot(fev$age[fev$subjectid == 1], fev$resid[fev$subjectid == 1],
+
+
type="l", lty=1, col=1, xlim = range(fev$age), ylim = range(fev$resid), xlab="Age (years)", ylab="log(FEV1/Height) residual")
> points(fev$age[fev$subjectid == 1], fev$resid[fev$subjectid == 1], + type="p", pch=1, col = 1, cex = 0 .6)
> for (i in 2:max(fev$subjectid)) {
+ lines(fev$age[fev$subjectid == i], fev$resid[fev$subjectid == i], + type = "l", lty = 1, col = i) + points(fev$age[fev$subjectid == i], fev$resid[fev$subjectid == i], + type = "p", pch = 1, col = i, cex = 0 .6)
+ }
> text(fev$age[fev$subjectid == 197] + 0 .8,
+ fev$resid[fev$subjectid == 197], "id 197")
> title("Spaghetti Plot: Residuals log(FEV1/Height) minus cubic spline")
Spaghetti Plot: Residuals log(FEV1/Height) minus cubic splin
id 197 |
18
Age (years)
Question 2
First round the ages:
> fev$ageround <- round(fev$age)
In order to use the pairs function, we need to get the data into ‘wide’ format, with one row per individual and a column for each age. We can do this using the reshape 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,use="pairwise .complete .obs")) + 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)
+ }
> fev .resid .wide .reshape <- reshape (
+ fev[,c("subjectid", "resid", "ageround")],
+ direction = "wide",
+ v .names = "resid",
+ timevar = "ageround",
+ idvar = "subjectid")
> head(fev .resid .wide .reshape)
subjectid resid .9 resid .10 resid .11 resid .12 resid .13
1
8
16
1 -0 .18275144 -0 .15477209 -0 .15288553 -0 .03693050 -0 .07181568
2 NA NA NA NA -0 .10286747 3 0 .14248385 0 .19768037 0 .25395274 NA -0 .04114325
25 4 -0 .13857584 -0 .06967834 -0 .06919457 -0 .10017548 -0 .18301857 35 5 NA NA NA NA NA 42 6 0.06273254 0.11656534 0.06271131 0.04941064 0.00250609
resid .15 resid .16 resid .7 resid .8 resid .14 resid .17
1 -0 .13769410 -0 .18744768 NA NA NA NA 8 0 .01590316 0 .07632892 0 .167993829 0 .06849432 -0 .096306567 -0 .168731997 16 0 .02461597 -0 .02089331 0 .219756542 0 .04612666 -0 .004626005 NA 25 -0 .03473378 NA -0 .114311932 -0 .12945013 0 .008968504 NA 35 -0 .16260641 -0 .11006280 -0 .084387049 -0 .12803007 -0 .057921453 -0 .155794884 42 0.02490978 0.06420085 0.001955933 0.08270616 0.125522304 0.004705282
resid .18 resid .6 resid .19
1 NA NA NA 8 0 .03833939 NA NA
16 NA NA NA 25 0 .05364655 NA NA 35 -0 .25670811 NA NA
42 NA NA NA
> pairs(fev .resid .wide .reshape[,- 1], upper .panel = panel .cor )
−0.3 −0.2 −0.4 −0.3 −0.4 −0.3 −0.20
2022-08-06