Chapter 2 Introduction

Why do the conventional methods fail in the presence of correlated data?

In basic neuroscience research, data dependency due to clustering or repeated measurements is probably the norm rather than the exception. Unfortunately, the most widely used methods under these situations are still the t-test and ANOVA, which do not take dependence into account, thus leading to incorrect statistical inference and misleading conclusions.

We will use a data set collected from our own work to assess the degree of data dependency due to clustering (animal effects) and to illustrate the consequences of ignoring the dependent structure. In this example, we measured the change in pCREB immunoreactivity of 1,200 putative excitatory neurons in mouse visual cortex at different time points: collected at baseline (saline), 24, 48, 72 hours, and 1 week following ketamine treatment, from 24 mice. See Grieco et al. (2020) for more details. Figure 2.1 shows that the changes in pCREB immunoreactivity tend to be clustered, i.e., measurements from the same animal tend to be more similar to each other than measurements from different animals.

knitr::include_graphics("Fig1S.jpg")
Normalized pCREB staining intensity values from 1,200 neurons (Example 1). The values in each cluster were from one animal. In total, pCREB values were measured for 1,200 neurons from 24 mice at five conditions: saline (7 mice), 24h (6 mice), 48h (3 mice), 72h (3 mice), 1week (5 mice) after treatment.

Figure 2.1: Normalized pCREB staining intensity values from 1,200 neurons (Example 1). The values in each cluster were from one animal. In total, pCREB values were measured for 1,200 neurons from 24 mice at five conditions: saline (7 mice), 24h (6 mice), 48h (3 mice), 72h (3 mice), 1week (5 mice) after treatment.

We compute the intra-class correlation (ICC) to quantify the magnitude of dependency within animals using the software R, a free and open source software (CRAN) (R Development Core Team, 2020). One major advantage of R over other open source or commercial software is that R is widely adopted and continuously reassessed for accuracy, and has a rich collection of user-contributed packages (over 15,000), thus supporting a programing environment for developers and access to cutting-edge statistical methods. In this tutorial, we will use the following R packages: lme4 (Bates et al. (2014)), __ (Jose Pinheiro et al. (2007)), icc (Wolak and Wolak (2015)), pbkrtest (Halekoh and Højsgaard (2014)), brms (Bürkner (2017b), Bürkner (2017a)), lmerTest (Kuznetsova, Brockhoff, and Christensen (2017)), emmeans (Lenth et al. (2019)), car (Fox and Weisberg (2018)) , and sjPlot (Lüdecke (2018)). If they have not been installed onto your computer, you will need to install them by removing the “#” symbol and copy one line at a time to your R console. The “#” symbol is used for commenting out code in R. The installation of a package to a computer only needs to be done once. However, the libraries for data analysis need to be loaded each time you start the R software. We recommend you only load a library when it is needed.

library(lme4) #if the library has been downloaded to your local computer, your can use install.pckages("lme4") to download it
library(nlme)
library(ICC)
library(brms)
library(pbkrtest)
library(emmeans)
library(car)
library(sjPlot)

We start with reading the pCREB data (Example 1) into R. Because the data file is comma separated, we use the function “read.csv” to read it. The option “head=T” reads the first row as the column names. Most R packages of LME require the “long,” also known as “vertical” format, in which data are organized in a rectangular data matrix, i.e., each row of the dataset contains only the values for one observation. The columns contain necessary information about this observation such as the experimental condition, treatment, cell ID, and animal ID. In this example, the data are stored in a 1,200-by-3 matrix, with the first column being the pCREB immunoreactivity values, the second column being the treatment labels, and the last column being the animal identification numbers. The treatment information is in the second column and it is coded as labels 1 through 5: 1 for baseline (saline), 2-5 for 24, 48, 72 hours, and 1 week after ketamine treatment, respectively. By default, the treatment information is read into numerical values. To convert it to a categorical variable, we apply the “as.factor” function to the treatment variable.

# The following lines of code read the Example 1 data
Ex1 = read.csv("https://www.ics.uci.edu/~zhaoxia/Data/BeyondTandANOVA/Example1.txt", head=T)
 
# checking the dimensions of the dataset
dim(Ex1)
## [1] 1200    3
# checking the names of each column
names(Ex1)
## [1] "res"           "treatment_idx" "midx"
# a frequency table for the treatment variable
table(Ex1$treatment_idx)
## 
##   1   2   3   4   5 
## 357 309 139 150 245
#Do not forget to factor the treatment IDs
Ex1$treatment_idx = as.factor(Ex1$treatment_idx)

Next, we examine the magnitude of clustering due to animal effects by computing the ICC for each treatment group.

### load the ICC library
library(ICC) #load the library to conduct ICC analysis with its function ICCbare
 
### conduct ICC analysis by organizing all the information into a data frame
icc.analysis=data.frame(n=rep(0,5), icc=rep(0,5), design.effect=rep(0,5),
effective.n=rep(0,5), M=rep(0,5), cells=rep(0,5))
 
for(i in 1:5){
    trt= Ex1[Ex1$treatment_idx==i,]
    trt$midx=factor(trt$midx)
    icc=ICCbare(factor(trt$midx), trt$res) #ICCbare is a function in the ICC package
    icc.analysis$cells[i]=dim(trt)[1]
    M=dim(trt)[1]/length(unique( trt$midx))
    def=1 + icc*(M-1)
    icc.analysis$n[i]=length(unique( trt$midx))
    icc.analysis$icc[i]=icc
    icc.analysis$design.effect[i]=def
    icc.analysis$effective.n[i]=dim(trt)[1]/def
    icc.analysis$M[i]=M
}
icc.analysis
##   n        icc design.effect effective.n        M cells
## 1 7 0.62094868     32.047434   11.139737 51.00000   357
## 2 6 0.33006327     17.668195   17.489053 51.50000   309
## 3 3 0.01780304      1.807071   76.920039 46.33333   139
## 4 3 0.62810904     31.777343    4.720344 50.00000   150
## 5 5 0.53694579     26.773398    9.150874 49.00000   245

The results are organized in the following table:

tmp=t(icc.analysis[,c(6,2)])
row.names(tmp)= c("# of cells", "ICC")
knitr::kable(
  tmp,
  col.names = c("Saline (7 mice)","24h (6 mice)","48h (3 mice)","72h (3 mice)","1wk (5 mice)")
)
Saline (7 mice) 24h (6 mice) 48h (3 mice) 72h (3 mice) 1wk (5 mice)
# of cells 357.0000000 309.0000000 139.000000 150.000000 245.0000000
ICC 0.6209487 0.3300633 0.017803 0.628109 0.5369458

The ICC indicates that the dependency due to clustering is substantial. Therefore, the 1,200 neurons should not be treated as 1,200 independent cells. When dependence is not adequately accounted for, the type I error rate can be much higher than the pre-chosen level of significance. To see how serious this problem is, we examine the false positives based on the dependence structure observed in our own study. In the simulation script we wrote “simulation_TypeIErrorRate.R” (see download section above), we generated 1000 data sets, each of which follows the same ICC structure and assumes NO difference between the five conditions. Surprisingly, the type I error rate when treating 1,200 neurons as independent observations is over 90% at the significance level of \(\alpha=0.05\).

### run the simulation script. It might take 1-2 minutes
source("https://www.ics.uci.edu/~zhaoxia/Data/BeyondTandANOVA/simulation_TypeIErrorRate.R")

This is a situation for which the number of observational units is much larger than the number of experimental units. We will show how to use a linear mixed-effects model to correctly analyze the data in the next section.

References

Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2014. “Fitting Linear Mixed-Effects Models Using Lme4.” arXiv Preprint arXiv:1406.5823.
Bürkner, Paul-Christian. 2017a. “Advanced Bayesian Multilevel Modeling with the r Package Brms.” arXiv Preprint arXiv:1705.11123.
———. 2017b. “Brms: An r Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software 80: 1–28.
Fox, John, and Sanford Weisberg. 2018. An r Companion to Applied Regression. Sage publications.
Grieco, Steven F, Xin Qiao, Xiaoting Zheng, Yongjun Liu, Lujia Chen, Hai Zhang, Zhaoxia Yu, et al. 2020. “Subanesthetic Ketamine Reactivates Adult Cortical Plasticity to Restore Vision from Amblyopia.” Current Biology 30 (18): 3591–3603.
Halekoh, Ulrich, and Søren Højsgaard. 2014. “A Kenward-Roger Approximation and Parametric Bootstrap Methods for Tests in Linear Mixed Models–the r Package Pbkrtest.” Journal of Statistical Software 59: 1–32.
Kuznetsova, Alexandra, Per B Brockhoff, and Rune HB Christensen. 2017. “lmerTest Package: Tests in Linear Mixed Effects Models.” Journal of Statistical Software 82: 1–26.
Lenth, R, H Singmann, J Love, P Buerkner, and M Herve. 2019. “Estimated Marginal Means, Aka Least-Squares Means. R Package Version 1.3. 2.”
Lüdecke, Daniel. 2018. “sjPlot: Data Visualization for Statistics in Social Science.” R Package Version 2 (1).
Pinheiro, Jose, Douglas Bates, Saikat DebRoy, Deepayan Sarkar, and R Core Team. 2007. “Linear and Nonlinear Mixed Effects Models.” R Package Version 3 (57): 1–89.
Wolak, Matthew, and Maintainer Matthew Wolak. 2015. “Package ‘ICC’.” Facilitating Estimation of the Intraclass Correlation Coefficient.