============================================================================================================ This up-to-date document is available at https://rpubs.com/sherloconan/1199857
Open Science Framework: https://osf.io/yv62c/
Disclaimer: I do not intend to teach or promote data manipulation or \(p\)-hacking. Please use this \(\texttt{R}\) markdown for discussion purposes only.
There are two types of data normalization used in the longitudinal data analysis:
\(\hspace{23pt}\frac{\text{Each Subject's Value}-({\color{red}{\text{Each Subject's Baseline at Time 0}}}-{\color{blue}{\text{Average Baseline at Time 0 Across All Subjects}}})}{{\color{blue}{\text{Average Baseline at Time 0 Across All Subjects}}}}\)
\(\hspace{23pt}\frac{\text{Each Subject's Value}}{{\color{red}{\text{Each Subject's Baseline at Time 0}}}}\)
Both methods will “normalize” the starting value to 1 for each subject. However the transformation does not make the data more normally distributed (Cousineau, 2019, p. 229).
The longitudinal data are from a simulation designed to study the effect of an intervention on the response measurements of \(n\) subjects, each of whom was repeatedly measured during and after the intervention. The data set (long format) consists of four columns:
ID - [factor] subject identifier; 15 levels;
Time - [numeric] time point, starting at 0;
Phase - [factor] treatment identifier; “2” as “intervention”, and “1” as “follow-up”;
Resp - [numeric] response measurements.
nSubj <- 15 # number of subjects
nTime <- 3 # number of time points in each phase
Trt <- c("intervention", "follow-up") # treatment phases
nTrt <- length(Trt)
set.seed(241) # try also 232
data <- data.frame("ID"=rep(paste0("s", 1:nSubj), each=nTime*nTrt),
"Time"=rep(0:(nTime*nTrt-1), times=nSubj),
"Phase"=rep(Trt, each=nTime, times=nSubj),
"Resp"=rnorm(nSubj*nTime*nTrt, mean=100, sd=10), # iid normal
stringsAsFactors=T)
# Note: no interaction effect in the true data-generating modeldatatable(data, filter="top")
\[\begin{equation} \tag{1} Y_{}=(\underset{\text{fixed intercept}}{b_0}+\underset{\text{random intercept}}{u_{0k}})+(\underset{\text{fixed slope}}{b_1}+\underset{\text{random slope}}{u_{1k}})\cdot\text{Time}_{}+b_2\cdot\text{Phase}_{}+b_3\cdot[\text{Phase}\times\text{Time}]_{}+\epsilon_{}, \end{equation}\]
for subject \(k=1,\dotsb,n\), \(\texttt{Phase}\in\{0,1\}\), and \(\texttt{Time}=0,1,\dotsb\).
|
\(p\)-Value (Raw) |
\(p\)-Value (Relative Change) |
\(p\)-Value (Subject-Centered) |
\(\textit{BF}_{01}\) (Raw) |
\(\textit{BF}_{01}\) (Relative Change) |
\(\textit{BF}_{01}\) (Subject-Centered) |
|
|---|---|---|---|---|---|---|
| \(\hat{b}_0\) | 0 | 0 | 0 | |||
| \(\hat{b}_2\) | 0.1261 | 0.0808 | 0.1100 | |||
| \(\hat{b}_1\) | 0.2300 | 0.2016 | 0.2210 | |||
| \(\hat{b}_3\) | 0.0856 | 0.0455 | 0.0725 | 1.67 | 2.01 | 1.06 |
# Fit the LME model to the raw data
fit_raw <- lme(fixed=Resp ~ Phase*Time,
random=~ Time | ID, data=data, # random intercept and slope for each subject
control=lmeControl(opt="optim", maxIter=1E3, msMaxIter=1E3))
summary(fit_raw)## Linear mixed-effects model fit by REML
## Data: data
## AIC BIC logLik
## 655.2483 674.8831 -319.6241
##
## Random effects:
## Formula: ~Time | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 2.5639459 (Intr)
## Time 0.1190702 -0.098
## Residual 8.8517581
##
## Fixed effects: Resp ~ Phase * Time
## Value Std.Error DF t-value p-value
## (Intercept) 107.36428 6.630840 72 16.191654 0.0000
## Phaseintervention -10.70776 6.919738 72 -1.547423 0.1261
## Time -1.95669 1.616395 72 -1.210526 0.2300
## Phaseintervention:Time 3.98394 2.285514 72 1.743128 0.0856
## Correlation:
## (Intr) Phsntr Time
## Phaseintervention -0.949
## Time -0.975 0.934
## Phaseintervention:Time 0.689 -0.826 -0.707
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.4214678 -0.6404955 0.1015791 0.4721396 2.1987894
##
## Number of Observations: 90
## Number of Groups: 15
Please note that relative changes cannot be compared across different subjects.
data <- data %>%
group_by(ID) %>%
mutate(Baseline = Resp[Time==0]) %>%
mutate(Relative_change = Resp / Baseline) %>%
select(-Baseline)# Now, fit the LME model to the normalized data
fit_relative <- lme(fixed=Relative_change ~ Phase*Time,
random=~ Time | ID, data=data,
control=lmeControl(opt="optim", maxIter=1E3, msMaxIter=1E3))
summary(fit_relative)## Linear mixed-effects model fit by REML
## Data: data
## AIC BIC logLik
## -113.1405 -93.50567 64.57023
##
## Random effects:
## Formula: ~Time | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.05375853 (Intr)
## Time 0.01423294 0.99
## Residual 0.08917850
##
## Fixed effects: Relative_change ~ Phase * Time
## Value Std.Error DF t-value p-value
## (Intercept) 1.1292647 0.06790354 72 16.630425 0.0000
## Phaseintervention -0.1234491 0.06971404 72 -1.770792 0.0808
## Time -0.0215091 0.01669127 72 -1.288641 0.2016
## Phaseintervention:Time 0.0468733 0.02302579 72 2.035689 0.0455
## Correlation:
## (Intr) Phsntr Time
## Phaseintervention -0.933
## Time -0.891 0.911
## Phaseintervention:Time 0.678 -0.826 -0.690
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.21501755 -0.67979964 0.08510331 0.54091821 2.20826421
##
## Number of Observations: 90
## Number of Groups: 15
pooled_mean_raw <- mean(subset(data, Time==0)$Resp) # mean raw baseline
data <- data %>%
group_by(ID) %>%
mutate(Baseline = Resp[Time==0]) %>%
mutate(Resp_centered = (Resp - Baseline + pooled_mean_raw) / pooled_mean_raw) %>% # subject-centering
select(-Baseline)fit_centered <- lme(fixed=Resp_centered ~ Phase*Time,
random=~ Time | ID, data=data,
control=lmeControl(opt="optim", maxIter=1E3, msMaxIter=1E3))
summary(fit_centered)## Linear mixed-effects model fit by REML
## Data: data
## AIC BIC logLik
## -116.8031 -97.16834 66.40156
##
## Random effects:
## Formula: ~Time | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.04785171 (Intr)
## Time 0.01409012 0.988
## Residual 0.08797127
##
## Fixed effects: Resp_centered ~ Phase * Time
## Value Std.Error DF t-value p-value
## (Intercept) 1.1156730 0.06672381 72 16.720764 0.0000
## Phaseintervention -0.1112694 0.06877031 72 -1.617986 0.1100
## Time -0.0203329 0.01646816 72 -1.234678 0.2210
## Phaseintervention:Time 0.0413990 0.02271409 72 1.822615 0.0725
## Correlation:
## (Intr) Phsntr Time
## Phaseintervention -0.937
## Time -0.899 0.911
## Phaseintervention:Time 0.681 -0.826 -0.690
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.26393017 -0.70223149 0.08104979 0.56120258 2.26098756
##
## Number of Observations: 90
## Number of Groups: 15
\(s_k\): Subject. \(\quad\alpha_i\): Time (treated as categorical). \(\quad\beta_j\): Phase. All have been standardized.
\[\begin{align} \tag{2} &\mathcal{M}_{1}:\ Y_{ijk}=\mu+\sigma_\epsilon\left(s_k+\alpha_i+\beta_j+(\alpha\beta)_{ij}+(\alpha s)_{ik}\right)+\epsilon_{ijk} \\ \text{versus}\qquad&\mathcal{M}_{0}:\ Y_{ijk}=\mu+\sigma_\epsilon\left(s_k+\alpha_i+\beta_j+(\alpha s)_{ik}\right)+\epsilon_{ijk},\qquad\text{for }\epsilon_{ijk}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,\sigma_\epsilon^2) \end{align}\]
\[\begin{equation} \tag{3} \pi(\mu,\sigma_\epsilon)\propto\frac{1}{\sigma_\epsilon} \end{equation}\]
\[\begin{equation} \tag{4} (\alpha_1^\star,\dotsb,\alpha_{a-1}^\star)=(\alpha_1,\dotsb,\alpha_a)\cdot\mathbf{Q},\quad\mathbf{I}_a-a^{-1}\mathbf{J}_a=\mathbf{Q}\cdot\mathbf{Q}^\top,\quad\dotsb \end{equation}\]
\[\begin{equation} \tag{5} \alpha_i^\star\mid g_A\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,g_A),\quad\beta_j^\star\mid g_B\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,g_B),\quad(\alpha\beta)_{ij}^\star\mid g_{AB}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,g_{AB}),\quad s_k\mid g_{s}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,g_s) \end{equation}\]
\[\begin{equation} \tag{6} g_A\sim\text{Scale-inv-}\chi^2(1,h_A^2),\quad g_B\sim\text{Scale-inv-}\chi^2(1,h_B^2),\quad g_{AB}\sim\text{Scale-inv-}\chi^2(1,h_{AB}^2),\quad g_s\sim\text{Scale-inv-}\chi^2(1,h_s^2) \end{equation}\]
By default, \(h=0.5\) for fixed effects and \(h=1\) for random effects (i.e., the standardized subject-specific random effects in this case).
The justification for modeling the mixed interaction \((\alpha s)_{ik}\) is not provided here. Please refer to Rouder et al. (2012, p. 365) for more details.
\(\mathbf{Q}\) is an \(a\times(a-1)\) matrix of the \(a-1\) eigenvectors of unit length corresponding to the nonzero eigenvalues of the left side term. For example, the projected main effect \(\alpha^\star=(\alpha_1-\alpha_2)/\sqrt{2}\) when \(a=2\). In the other direction (given \(\alpha_1+\alpha_2=0\)), \(\alpha_1=\alpha^\star/\sqrt{2}\) and \(\alpha_2=-\alpha^\star/\sqrt{2}\). The same reduced parametrization applies to \(\beta_j\).
As for the interaction effects, four variables are projected into one. \((\alpha\beta)^\star=((\alpha\beta)_{11}-(\alpha\beta)_{21}-(\alpha\beta)_{12}+(\alpha\beta)_{22})/2\). In the other direction (given \((\alpha\beta)_{11}+(\alpha\beta)_{12}=0\), \((\alpha\beta)_{21}+(\alpha\beta)_{22}=0\), and \((\alpha\beta)_{11}+(\alpha\beta)_{21}=0\)), \((\alpha\beta)_{11}=(\alpha\beta)^\star/2\), \((\alpha\beta)_{21}=-(\alpha\beta)^\star/2\), \((\alpha\beta)_{12}=-(\alpha\beta)^\star/2\), and \((\alpha\beta)_{22}=(\alpha\beta)^\star/2\).
df <- as.data.frame(data)
df$Time <- factor(df$Time) # non-parsimonious
set.seed(277)
full <- lmBF(Resp~Time*Phase+Time:ID+ID, df, whichRandom="ID", progress=F)
set.seed(277)
additive <- lmBF(Resp~Time+Phase+Time:ID+ID, df, whichRandom="ID", progress=F)
full / additive # BF10 = 0.597396 ±4.98%; BF01 = 1.673932 ±4.98%## Bayes factor analysis
## --------------
## [1] Time * Phase + Time:ID + ID : 0.597396 ±4.98%
##
## Against denominator:
## Resp ~ Time + Phase + Time:ID + ID
## ---
## Bayes factor type: BFlinearModel, JZSset.seed(277)
full_relative <- lmBF(Relative_change~Time*Phase+Time:ID+ID, df, whichRandom="ID", progress=F)
set.seed(277)
additive_relative <- lmBF(Relative_change~Time+Phase+Time:ID+ID, df, whichRandom="ID", progress=F)
full_relative / additive_relative # BF10 = 0.4970184 ±8.43%; BF01 = 2.011998 ±8.43%## Bayes factor analysis
## --------------
## [1] Time * Phase + Time:ID + ID : 0.4970184 ±8.43%
##
## Against denominator:
## Relative_change ~ Time + Phase + Time:ID + ID
## ---
## Bayes factor type: BFlinearModel, JZSset.seed(277)
full_centered <- lmBF(Resp_centered~Time*Phase+Time:ID+ID, df, whichRandom="ID", progress=F)
set.seed(277)
additive_centered <- lmBF(Resp_centered~Time+Phase+Time:ID+ID, df, whichRandom="ID", progress=F)
full_centered / additive_centered # BF10 = 0.9458405 ±9.02%; BF01 = 1.057261 ±9.02%## Bayes factor analysis
## --------------
## [1] Time * Phase + Time:ID + ID : 0.9458405 ±9.02%
##
## Against denominator:
## Resp_centered ~ Time + Phase + Time:ID + ID
## ---
## Bayes factor type: BFlinearModel, JZS
See also https://cran.r-project.org/web/packages/BayesFactor/vignettes/manual.html#generalTestBF
df <- as.data.frame(data)
set.seed(277)
generalTestBF(Resp~Time*Phase+Time:ID+ID, df, whichRandom="ID", neverExclude="^ID$", progress=F)## Warning in doNwaySampling(method, y, X, rscale, iterations, gMap, incCont, :
## Some NAs were removed from sampling results: 4407 in total.## Warning in doNwaySampling(method, y, X, rscale, iterations, gMap, incCont, :
## Some NAs were removed from sampling results: 4334 in total.## Error in optim(qs, Qg, gr = dQg, control = list(fnscale = -1), method = "BFGS", :
## initial value in 'vmmin' is not finite## Warning in doNwaySampling(method, y, X, rscale, iterations, gMap, incCont, :
## Some NAs were removed from sampling results: 4597 in total.## Bayes factor analysis
## --------------
## [1] Time + ID : 0.02608702 ±0.84%
## [2] Phase + ID : 0.03021746 ±6.07%
## [3] Time + Phase + ID : 0.009438913 ±2.02%
## [4] Time + Phase + Time:Phase + ID : 0.01214553 ±1.23%
## [5] Time + Time:ID + ID : NaNe+Inf ±NaN%
## [6] Time + Phase + Time:ID + ID : 1.177706e-07 ±99.98%
## [7] Time + Phase + Time:Phase + Time:ID + ID : NaNe+Inf ±0%
## [8] ID : 0.112958 ±0%
##
## Against denominator:
## Intercept only
## ---
## Bayes factor type: BFlinearModel, JZS
Normalization makes data dimensionless, meaning that the features are not dependent on the original units of measurement.
For example, \(\texttt{Google Trends}\) tracks popularity by representing each data point as its sampled search volume (multiplied by 100) divided by the highest sampled search volume within the specified filter settings and time range, a method known as baseline normalization. The data always contain 100 as the maximum value (baseline). Additionally, multiple sets of popularity data are comparable only when they share the same baseline.
Different normalization techniques can produce different results, and choosing the wrong method for a specific data set or analysis can negatively impact the outcomes. For example, Min-Max normalization is sensitive to outliers, which can skew the scaling process.
By normalizing data in different ways, researchers can inadvertently or intentionally find a version that yields significant results, even if those results are not truly reflective of underlying patterns. Read about the Gino case; post, video.
In transcriptomics, gene expression levels (logcounts) are often normalized and transformed to account for various technical biases and to make the data more comparable across different samples.
NormalizeData: https://satijalab.org/seurat/reference/normalizedata
SCTransform: https://satijalab.org/seurat/reference/sctransform
Similarly, in experimental psychology, it is common to re-parameterize effect size as \(t_i = \tau_i\ /\ \sigma_\epsilon\) so that the treatment effects are standardized relative to the standard deviation of the error and become dimensionless (Jeffreys, 1961; Rouder et al., 2012, p. 359; Wei, Nathoo, & Masson, 2023, p. 1762; but see also https://rpubs.com/sherloconan/1187595).
Scale invariance refers to the property of an algorithm or model to be unaffected by the scale of the input data. The value of the Bayes factor is unaffected by multiplicative changes in the unit of measure of the observations. For instance, if observations are in a unit of length, the Bayes factor is the same whether the measurement is in nanometers or light-years (Rouder et al., 2012, p. 360).
However, there is an ongoing debate and a round-table discussion about the standardization of effect sizes in van Doorn et al. (2023).
Many researchers refer to \(t_i\) as an obfuscating and decontextualized standardized effect size, which adds a layer of abstraction to its interpretation and reification, especially in mixed-effects models.
They argued that \(\hat{\sigma}_\epsilon\) can be contaminated by shrinkage. Specifically, if shrinkage is applied to the random slope estimate, then \(\hat{\sigma}_\epsilon\) is consequently contaminated (a negative relationship between \(\hat{\sigma}_\epsilon^2\cdot\hat{g}\) and \(\hat{\sigma}_\epsilon^2\); p. 153-154).
Also, consider the one-way between-subjects design with a heteroscedastic error covariance matrix and the unreliable Bayes factor estimates when there are many divergent transitions https://rpubs.com/sherloconan/1075410.
Real vs Vector parameter(s): one is working and the other is not (as expected). Why? Jun 9, 2022 • https://discourse.mc-stan.org/t/real-vs-vector-parameter-s-one-is-working-and-the-other-is-not-as-expected-why/27784
Writing the Jeffreys prior for a vector variable in Stan? Jul 26, 2022 • https://discourse.mc-stan.org/t/jeffreys-prior-for-a-vector-variable/28338 • https://stats.stackexchange.com/questions/535955/how-do-i-write-the-jeffreys-prior-for-error-variance-in-stan-p-mu-sigma-12
How do I write the Jeffreys prior for error variance (heteroscedasticity)? Jul 28, 2021 • https://discourse.mc-stan.org/t/how-do-i-write-the-jeffreys-prior-for-error-variance-heteroscedasticity/23697 • https://stats.stackexchange.com/questions/531632/is-my-stan-model-correct-the-jeffreys-prior-for-a-heteroscedastic-mixed-effects?noredirect=1#comment976543_531632
#########################
# Simulate Example Data #
#########################
nSubj <- 15 # number of subjects
nTime <- 3 # number of time points in each phase
Trt <- c("intervention", "follow-up") # treatment phases
nTrt <- length(Trt)
set.seed(241) # try also 232
data <- data.frame("ID"=rep(paste0("s", 1:nSubj), each=nTime*nTrt),
"Time"=rep(0:(nTime*nTrt-1), times=nSubj),
"Phase"=rep(Trt, each=nTime, times=nSubj),
"Resp"=rnorm(nSubj*nTime*nTrt, mean=100, sd=10), # iid normal
stringsAsFactors=T)
# Note: no interaction effect in the true data-generating model
View(data)
##############################
# Linear Mixed-Effects Model #
##############################
library(nlme) #v 3.1-165
# Fit the LME model to the raw data
fit_raw <- lme(fixed=Resp ~ Phase*Time,
random=~ Time | ID, data=data, # random intercept and slope for each subject
control=lmeControl(opt="optim", maxIter=1E3, msMaxIter=1E3))
summary(fit_raw)
# Fixed effects: Resp ~ Phase * Time
# Value Std.Error DF t-value p-value
# (Intercept) 107.36428 6.630840 72 16.191654 0.0000
# Phaseintervention -10.70776 6.919738 72 -1.547423 0.1261
# Time -1.95669 1.616395 72 -1.210526 0.2300
# Phaseintervention:Time 3.98394 2.285514 72 1.743128 0.0856 <-- not statistically significant
#################
# Normalization #
#################
library(dplyr) #v 1.1.4
# Normalize response by baseline value at Time 0 for each subject
data <- data %>%
group_by(ID) %>%
mutate(Baseline = Resp[Time==0]) %>%
mutate(Relative_change = Resp / Baseline) %>%
select(-Baseline)
subset(data, Time==0) # check the processed data
# Now, fit the LME model to the normalized data
fit_relative <- lme(fixed=Relative_change ~ Phase*Time,
random=~ Time | ID, data=data,
control=lmeControl(opt="optim", maxIter=1E3, msMaxIter=1E3))
summary(fit_relative)
# Fixed effects: Relative_change ~ Phase * Time
# Value Std.Error DF t-value p-value
# (Intercept) 1.1292647 0.06790354 72 16.630425 0.0000
# Phaseintervention -0.1234491 0.06971404 72 -1.770792 0.0808
# Time -0.0215091 0.01669127 72 -1.288641 0.2016
# Phaseintervention:Time 0.0468733 0.02302579 72 2.035689 0.0455 <-- statistically significant
## Disclaimer: I do not intend to teach or promote data manipulation or p-hacking.
## Please use the example R script for discussion purposes only.
pooled_mean_raw <- mean(subset(data, Time==0)$Resp) # mean raw baseline
data <- data %>%
group_by(ID) %>%
mutate(Baseline = Resp[Time==0]) %>%
mutate(Resp_centered = (Resp - Baseline + pooled_mean_raw) / pooled_mean_raw) %>% # subject-centering
select(-Baseline)
fit_centered <- lme(fixed=Resp_centered ~ Phase*Time,
random=~ Time | ID, data=data,
control=lmeControl(opt="optim", maxIter=1E3, msMaxIter=1E3))
summary(fit_centered)
# Fixed effects: Resp_centered ~ Phase * Time
# Value Std.Error DF t-value p-value
# (Intercept) 1.1156730 0.06672381 72 16.720764 0.0000
# Phaseintervention -0.1112694 0.06877031 72 -1.617986 0.1100
# Time -0.0203329 0.01646816 72 -1.234678 0.2210
# Phaseintervention:Time 0.0413990 0.02271409 72 1.822615 0.0725
#################
# Bayes Factors #
#################
library(BayesFactor) #v 0.9.12-4.7
df <- as.data.frame(data)
set.seed(277)
(bf <- generalTestBF(Resp~Time*Phase+Time:ID+ID, df, whichRandom="ID", neverExclude="^ID$", progress=F))
df$Time <- factor(df$Time) # non-parsimonious
set.seed(277)
full <- lmBF(Resp~Time*Phase+Time:ID+ID, df, whichRandom="ID", progress=F)
set.seed(277)
additive <- lmBF(Resp~Time+Phase+Time:ID+ID, df, whichRandom="ID", progress=F)
full / additive # BF10 = 0.597396 ±4.98%; BF01 = 1.673932 ±4.98%
set.seed(277)
full_relative <- lmBF(Relative_change~Time*Phase+Time:ID+ID, df, whichRandom="ID", progress=F)
set.seed(277)
additive_relative <- lmBF(Relative_change~Time+Phase+Time:ID+ID, df, whichRandom="ID", progress=F)
full_relative / additive_relative # BF10 = 0.4970184 ±8.43%; BF01 = 2.011998 ±8.43%
set.seed(277)
full_centered <- lmBF(Resp_centered~Time*Phase+Time:ID+ID, df, whichRandom="ID", progress=F)
set.seed(277)
additive_centered <- lmBF(Resp_centered~Time+Phase+Time:ID+ID, df, whichRandom="ID", progress=F)
full_centered / additive_centered # BF10 = 0.9458405 ±9.02%; BF01 = 1.057261 ±9.02%