Tom,
Perhaps the following code will do the trick :
ICC <- lm(x~j + i,data=iccdata)
edf <- df.residual(ICC)
ems <- deviance(ICC)/edf
bdf <- anova(ICC)[1,1]
bms <- anova(ICC)[1,2]/bdf
msb <- bms
jdf <- anova(ICC)[2,1]
jms <- anova(ICC)[2,2]/jdf
k <- jdf + 1
msw <- ((ems*edf)+(jms*jdf))/(edf+jdf)
wms <- msw
n <- bdf + 1
theta <- (msb - msw)/(k*msw)
winer11 <- theta / (1+theta)
winer1k <- (k*theta)/(1+k*theta)
sf11 <- (bms - wms)/(bms+(k-1)*wms)
sf21 <- (bms - ems)/((bms)+((k-1)*ems)+((k*(jms-ems))/n))
sf31 <- (bms - ems)/(bms+((k-1)*ems))
sf1k <- (bms - wms)/bms
sf1k <- (bms - ems)/(bms+((jms-ems)/n))
sf3k <- (bms - ems)/bms
Invoke a few print commands for the output desired and voila, Shrout and
Fleiss ICC's right from R. I might have a few typos in the above code
but the general picture works. I tested the critical parts. In
addition, the parceling of the anova output is rather clunky. I hope
others will show us a more efficient method of extracting the sums of
squares. Let me know if the results are consistent with your hand
calculations.
Cheers,
Patrick
NB: The general approach above was abstracted from Robert M. Hamer's
SAS macro.
__________________________________________________________________________________________
Patrick E. McKnight
Evaluation Group for Analysis of Data
University of Arizona
Department of Psychology
Tucson, AZ 85721
520-621-5463
pem at theriver.com
Original message....
Hello, R-folks:
Here is a statement I use to make a data frame:
iccdata <- data.frame(i=rep(1:10,rep(2,10)),j=rep(1:2,10),
x=c(0.35011,0.11989,0.13081,0.09919,0.16000,0.12000,0.00000,0.00000,
0.44023,0.32977,2.67081,2.63919,0.09050,0.03950,0.44019,0.30981,0.59000,
0.57000,4.03000,3.77000))
<snip>
Now, how do you use (lme in) R to compute the Intraclass correlation for
these data? I have written brute-force code to get the "right"
answer, but
I cannot get the proper mean squares from any regression model I have tried.
It seems like it should be trivial, but I cannot do it! I need this to work
in order to apply the technique to a much larger data set. I could use my
silly code, in SAS, but I want to see and understand the variance
components. After studying Pinheiro and Bates for a while, I naively tried
the following, which is not what I want.
>> ICC <- groupedData(x~j|i,data=iccdata)
>> ICC.lme <- lme(x~j,data=ICC)
>
Can you help me? Thanks in advance!
Tom
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._