Intraclass correlation example

# Load pacman(載入pacman, tidyverse,lme4套件)
pacman::p_load(tidyverse, lme4)

# input data (載入資料)
dta <- read.table("C:/Users/ASUS/Desktop/data/finger_ridges.txt", header = T)

# coerce pair variable to factor type (將pair變項轉為類別格式)
dta$Pair <- factor(dta$Pair)

# show first 6 lines(顯示資料前6列)
head(dta)

##   Pair Twin_1 Twin_2
## 1    1     71     71
## 2    2     79     82
## 3    3    105     99
## 4    4    115    114
## 5    5     76     70
## 6    6     83     82

# historical idea of intraclass correlation (計算ICC值)
cor(c(dta$Twin_1, dta$Twin_2), c(dta$Twin_2, dta$Twin_1))

## [1] 0.9625626

# turn data from wide format to long format (將資料由寬格式轉為長格式)
dtaL <- dta %>% 
        gather(key = Twin, value = Ridges, Twin_1, Twin_2) %>%
        mutate( Twin = as.factor(Twin)) 
 
# show mean and sd of ridges by twin pair  (列出每組雙胞胎指紋數的平均和標準差)     
dtaL %>% group_by(Pair) %>%
         summarize( ridge_ave = mean(Ridges), ridge_sd = sd(Ridges)) %>%
         arrange( desc(ridge_ave)) 

## `summarise()` ungrouping output (override with `.groups` argument)

## # A tibble: 12 x 3
##    Pair  ridge_ave ridge_sd
##    <fct>     <dbl>    <dbl>
##  1 4         114.     0.707
##  2 7         114.     0.707
##  3 9         114.     0.707
##  4 3         102      4.24 
##  5 10         92.5    2.12 
##  6 6          82.5    0.707
##  7 2          80.5    2.12 
##  8 11         79      5.66 
##  9 12         74      2.83 
## 10 5          73      4.24 
## 11 1          71      0    
## 12 8          50.5    9.19

# dot plot (繪圖 y為指紋數，x為每組雙胞胎)
ggplot(dtaL, aes(reorder(Pair, Ridges, mean), Ridges, color = Pair)) + 
 geom_point() + 
 geom_hline(yintercept = mean(dtaL$Ridges), linetype = "dashed") +
 labs(y = "Number of finger ridges", x = "Twin ID") +
 theme_bw() +
 theme(legend.position = "NONE")

# random effect anova for balanced designs (算出雙胞胎factor的殘差及SSW)
summary(m0 <- aov(Ridges ~ Error(Pair), data = dtaL))

## 
## Error: Pair
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 11   8990   817.3               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 12  171.5   14.29

# random effects anova (呈現random與fixed effect model結果)
summary(m1 <- lmer(Ridges ~ 1 + (1 | Pair), data = dtaL))

## Linear mixed model fit by REML ['lmerMod']
## Formula: Ridges ~ 1 + (1 | Pair)
##    Data: dtaL
## 
## REML criterion at convergence: 174.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.88917 -0.38618 -0.00834  0.37963  1.54959 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Pair     (Intercept) 401.51   20.04   
##  Residual              14.29    3.78   
## Number of obs: 24, groups:  Pair, 12
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   87.208      5.836   14.94

#呈現每組雙胞胎指紋平均值
ranef(m1)

## $Pair
##    (Intercept)
## 1   -15.924912
## 2    -6.591031
## 3    14.533018
## 4    26.814441
## 5   -13.959885
## 6    -4.626003
## 7    25.831927
## 8   -36.066447
## 9    25.831927
## 10    5.199136
## 11   -8.064801
## 12  -12.977371
## 
## with conditional variances for "Pair"

knitr::spin("C:\\Users\\ASUS\\Desktop\\finger_ridges.R", knit=FALSE)

## [1] "C:\\Users\\ASUS\\Desktop\\finger_ridges.Rmd"

# The end