a. Follow the Week 2 R-session and obtain the plot showing each subjects data and straight-line fit. Use lmList to obtain the 40 slopes for the straight-line fits. Compare the five-number summary of rates of change for the "X" measurement with that obtained for the perfectly measured "Xi" measurements in the posted R-session. --------------------- > week1X = read.table(file="D:\\drr12\\stat222\\week2\\mythdata_X", header = T) > head(week1X) subject X.1 X.3 X.5 W 1 1 37.51632 51.35238 59.44765 15.97247 2 2 45.12749 52.81792 61.64658 15.37724 3 3 35.14619 56.82575 66.15056 11.47902 4 4 44.12592 49.18999 64.57075 16.88944 5 5 52.74255 66.55824 70.48820 19.17834 6 6 30.42937 49.95363 64.29086 11.81822 > cor(week1X) X.1 X.3 X.5 W X.1 1.0000000 0.7158619 0.4906650 0.7162989 X.3 0.7158619 1.0000000 0.7405647 0.6522400 X.5 0.4906650 0.7405647 1.0000000 0.6192915 W 0.7162989 0.6522400 0.6192915 1.0000000 > attach(week1X) > diff = X.5 - X.1 > cor(diff, X.1) [1] -0.1993601 > cor(diff, W) [1] 0.1584012 > cor.test(diff, W) Pearson's product-moment correlation data: diff and W t = 0.9889, df = 38, p-value = 0.3289 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: -0.1610551 0.4478139 sample estimates: cor 0.1584012 > # note cor(W, theta) = 0 > library(reshape) Loading required package: plyr Attaching package: ‘reshape’ The following object(s) are masked from ‘package:plyr’: rename, round_any Warning messages: 1: package ‘reshape’ was built under R version 2.14.2 2: package ‘plyr’ was built under R version 2.14.2 > detach(week1X) > # For the time being I want to ignore the "W" column and just deal with the time-ordered observations > > dat.wide = subset(week1X, select = c(subject, X.1 , X.3, X.5 )) > # rename the observation headers by the numerical value of timepoint > names(dat.wide) = c("subject", "1", "3", "5") > dat.molten = melt(dat.wide, id.vars = c("subject")) > head(dat.molten) subject variable value 1 1 1 37.51632 2 2 1 45.12749 3 3 1 35.14619 4 4 1 44.12592 5 5 1 52.74255 6 6 1 30.42937 > dat.long = cast(dat.molten, subject + variable ~ .) > names(dat.long) = c("subject", "time", "X") > head(dat.long) subject time X 1 1 1 37.51632 2 1 3 51.35238 3 1 5 59.44765 4 2 1 45.12749 5 2 3 52.81792 6 2 5 61.64658 > # write out this "long" dataset (clears out any baggage from the 'melt' process) > write.table(dat.long, file = "D:\\drr12\\stat222\\week2\\xlong.dat", quote = FALSE) > week1X.long = read.table(file="D:\\drr12\\stat222\\week2\\xlong.dat", header = T) > library(lme4) Loading required package: Matrix Loading required package: lattice Attaching package: ‘Matrix’ The following object(s) are masked from ‘package:reshape’: expand The following object(s) are masked from ‘package:base’: det Attaching package: ‘lme4’ The following object(s) are masked from ‘package:stats’: AIC, BIC > xList = lmList(X ~ time |subject, data = week1X.long) #fit straight-line to each subject > rate = coef(xList)[2] > fivenum(rate) [1] 0.5290025 3.9306175 5.1907113 6.1846750 10.6300700 > # from file XiRsession (linked) > # fivenum(rate[,1]) #5-number summary for Xi data > #[1] 1.095690 3.989401 5.033039 6.336699 8.991205 > # so as you would expect extremes wider for "X" data > # make plot of each subject's data and fit > xyplot(X ~ time | subject, type=c("p", "r"), index.cond=function(x,y) + + {coef(lm(y ~ x))[1]}, data=week1X.long) >