In class Exercise3

1 load data

pacman::p_load(tidyverse, lme4)

# input data
dta <- read.table("finger_ridges.txt", header = T)

2 change pair variable to factor type

# coerce pair variable to factor type
dta$Pair <- factor(dta$Pair)
# show first 6 lines
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

3 intraclass correlation

# historical idea of intraclass correlation
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)) 
#寬變長

4 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)) 

## # 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

#算平均與1標準差 依照ridge_ave遞減排序
#可知第4對的雙胞胎手的Ridges較多

5 dot plot

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")

#每個觀察值以點圖型式呈現

6 random effect anova for balanced designs

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
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

#Estimated Var = (817.3 - 14.29) /2 = 401.51
#組內相關係數(ICC) = 401.51/(401.51+14.29) = 0.966

ranef(m1)

## $Pair
##    (Intercept)
## 1   -15.924913
## 2    -6.591031
## 3    14.533019
## 4    26.814443
## 5   -13.959886
## 6    -4.626003
## 7    25.831929
## 8   -36.066449
## 9    25.831929
## 10    5.199136
## 11   -8.064802
## 12  -12.977372
## 
## with conditional variances for "Pair"

#算12組雙胞胎的隨機效果