8
$\begingroup$

I have 7 patients with blood measurements from two different tissues and at two different oxygenation states. We first took measurements from tissue A and B (Normoxia) and then we administrated oxygen (Hyperoxia) and took measurement from tissue A and tissue B. I would like to see the effect of oxygen administration as well as if there is a difference between the two tissues.

Initially, I performed a paired t-test between tissue A and tissue B in Normoxia state and a second paired t-test between tissue A and tissue B in Hyperoxia. With this way I can see if there is a difference between tissue types in each oxygenation state (Normoxia vs Hyperoxia) but I cannot see the effect of hyperoxia. Should I perform instead a One-Way ANOVA between the 4 measurements I have?

I have added my data below.

           Normoxia            Hyperoxia
     Tissue A  Tissue B     Tissue A Tissue B
1    0.0014639 0.0013601    0.001512 0.001464
2    0.0013740 0.0013480    0.001638 0.001518
3    0.0014739 0.0014101    0.001745 0.001559
4    0.0015490 0.0015230    0.001585 0.001531
5    0.0013931 0.0013829    0.001498 0.001400
6    0.0015456 0.0014864    0.001655 0.001599
7    0.0015980 0.0015920    0.001764 0.001672

image of data

$\endgroup$
7
  • 3
    $\begingroup$ Is each row an independent sample? Or is the design further complicated by sampling multiple sites within a tissue specimen for multiple tissue samples? $\endgroup$
    – AdamO
    Commented Jun 18 at 13:24
  • $\begingroup$ Each row is an independent sample $\endgroup$
    – user161260
    Commented Jun 18 at 13:41
  • $\begingroup$ When you say "a difference between the two tissues" are you referring to a baseline difference or any difference in the effect of hyperoxia on the tissues? (I will note that a graph of the data is often much more informative as a "picture" than a table of values.) $\endgroup$ Commented Jun 18 at 21:16
  • $\begingroup$ The standard analysis for this sort of balanced repeated measures data setup is that described below by @dariober. It gives exactly the same effect as a paired analysis but with four repeated measures instead of two and allows you to compare hyperoxia and normoxia. $\endgroup$ Commented Jun 19 at 0:26
  • $\begingroup$ There is no advantage is multilevels models as recommended by others on this thread because (1) all the comparisons of interest are within-patient, so only one level is relevant and (2) for a balanced design like this there is no between-patient information to be recovered from random effects. So the patient blocking approach has just as much power as more elaborate approaches while being simpler and more robust. $\endgroup$ Commented Jun 19 at 8:30

6 Answers 6

6
$\begingroup$

I feel I'm making things too easy here... In addition to estimating the effect of tissues and the effect of oxygen level, you want to test whether there is an interaction between tissue and oxygen. That is, you ask whether tissue A responds differently than B to oxygen.

In R code this could be fitted as:

fit <- lm(y ~ patient + tissue * treatmeant)
summary(fit) # or use package emmeans for more sophisticated output

(Posting the data as text instead of picture would help people playing with it)


Some simulated data:

pt <- sprintf('pt%03d', 1:7)
tissue <- c('A', 'B')
ox <- c('N', 'H')

dat <- data.frame(
    pt = rep(pt, each=4),
    tissue = rep(rep(tissue, 2), length(pt)),
    ox = rep(rep(ox, each=2), length(pt))
)

Add response variable with effect for tissue and oxygen but no interaction:

set.seed(1234)
dat$y <- (as.numeric(as.factor(dat$tissue)) + as.numeric(as.factor(dat$ox))) + rnorm(n=nrow(dat), sd=0.1)

enter image description here

summary(lm(y ~ pt + tissue * ox, data=dat))

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.9777     0.0618   31.99  < 2e-16 ***
ptb           0.0501     0.0732    0.69     0.50    
ptc          -0.0185     0.0732   -0.25     0.80    
ptd           0.0582     0.0732    0.80     0.44    
pte           0.0587     0.0732    0.80     0.43    
ptf           0.0463     0.0732    0.63     0.53    
ptg          -0.0100     0.0732   -0.14     0.89    
tissueB       0.9652     0.0553   17.45  1.0e-12 ***
oxN           0.9503     0.0553   17.18  1.3e-12 ***
tissueB:oxN   0.0391     0.0782    0.50     0.62    

Now a response with interaction:

set.seed(1234)
dat$y <- (as.numeric(as.factor(dat$tissue)) * as.numeric(as.factor(dat$ox))) + rnorm(n=nrow(dat), sd=0.1)

enter image description here

summary(lm(y ~ pt + tissue * ox, data=dat))

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.9777     0.0618   15.81  5.3e-12 ***
ptb           0.0501     0.0732    0.69     0.50    
ptc          -0.0185     0.0732   -0.25     0.80    
ptd           0.0582     0.0732    0.80     0.44    
pte           0.0587     0.0732    0.80     0.43    
ptf           0.0463     0.0732    0.63     0.53    
ptg          -0.0100     0.0732   -0.14     0.89    
tissueB       0.9652     0.0553   17.45  1.0e-12 ***
oxN           0.9503     0.0553   17.18  1.3e-12 ***
tissueB:oxN   1.0391     0.0782   13.29  9.6e-11 ***

Code for plot

library(ggplot2)

gg <- ggplot(data=dat, aes(x=tissue, y=y, colour=ox, group=paste(pt, ox))) +
    geom_point() +
    geom_line() +
    theme_light()
$\endgroup$
2
  • 1
    $\begingroup$ +1. The model ~patient + group where group has four levels would also be effective and perhaps more flexible. BTW, you've left the ~ out from the lm model formula. $\endgroup$ Commented Jun 19 at 0:28
  • $\begingroup$ @GordonSmyth lm formula fixed, thanks. $\endgroup$
    – dariober
    Commented Jun 19 at 7:51
4
$\begingroup$

Neither of those. The first one, as you note, doesn't let you look at the difference in conditions, which you probably want to do. The second one ignores the "matching" or repeated measures, and thus violates the assumption of independent errors.

I would recommend a multilevel model, with "treatment" and "location" as fixed effects and "person" as a random effect. RM-ANOVA is another possibility, but that assumes sphericity, which may not be reasonable and which will be hard to examine with such a small. N.

$\endgroup$
2
  • $\begingroup$ Could you please recommend a multilevel model to use? I am familiar with R and GraphPad Prism. $\endgroup$
    – user161260
    Commented Jun 18 at 13:42
  • $\begingroup$ R has lots of resources for this. This page is a good starting place $\endgroup$
    – Peter Flom
    Commented Jun 18 at 15:09
3
$\begingroup$

ANOVA can be tricked into solving this kind of problem by calculating the paired differences for $k$-independent samples, $k$ being the number of independent samples. Just like with the paired t-test, the corresponding paired ANOVA would test whether the intercept and all the regressors are simultaneously equal to 0 - linear regression is yet another method to achieve this.

The problem with multiple t-tests is that it ignores multiplicity unless you adjust the $p$-value - like with this relevant XKCD.

There are many ways to adjust $p$-values such as the Bonferroni, Benjamini-Hochberg, or the Tukey HSD methods. On the other hand, multiple t-tests are favorable because they account for heteroscedasticity more efficiently than ANOVA.

$\endgroup$
2
$\begingroup$

Unless I am missing something, the answer is in fact rather simple, and depends completely on the independence of the samples; is the tissue sample under "Normoxia, Tissua A, row 1", the exact same sample as under "Hyperoxia, Tissue A, row 1" (and similarly for the other data points), or are these 2 independent tissue samples (coming both from Tissue type A, but independent)? I interpret the OP's answer to a previous commnent as the samples being independent?
If they are indeed independent, then a 2-way ANOVA (1 numerical DV, 2 categorical IV's) will do the job.
If the tissue samples are not independent (in fact, paired), then compute the paired differences, and run a one-way ANOVA on the paired difference (factor being tissue type), which in fact becomes just a 2-sample t-test.
Fwiw, in the first case, your treatment is significant (p=0.003), but not the tissue type (p=.064...) (assuming "traditional" .05 significance). In the 2nd case, tissue is not significant again (p=0.26). But the treatment is significant (p=0.000).

$\endgroup$
2
  • $\begingroup$ The samples are independent so I plan to use a 2-way ANOVA. Shall I also run a post-hoc Tukey analysis to see which combination of variables are significant? $\endgroup$
    – user161260
    Commented Jun 20 at 8:24
  • $\begingroup$ The 2-way ANOVA alone will only tell you that there are significant differences between the means; we know this just from graphing the data. Howver, the linear model from the 2-way ANOVA will also tell you which factor is significant, and whether there are interactions (i.e. a tissue type has an larger response, on average, to the treatment). So do the 2-way ANOVA (basically factorial ANOVA) and look at the linear model, and the adjusted p-values for the various terms. $\endgroup$
    – jginestet
    Commented Jun 20 at 16:14
1
$\begingroup$

The Problem

You want to see the effect of oxygen administration and whether that difference is affected by tissue type. This can be done with a linear model using an interaction, or by using a factorial ANOVA. These are two sides of the same coin; they will give you identical answers. [1], [2]

The challenge is how to obtain fair $p$-values, because these aren't 28 independent observations, but rather 4 measurements from each of 7 patients.

Two Solutions

One type of model that can do that is a mixed model (this is the suggestion made by Peter Flom).$^\dagger$

In R, this looks something like this:

require("lmerTest")
LMM <- lmer(y ~ Tissue * State + (1 | Patient), data = DF)

As has been pointed out by others already, since you happen to have a balanced design (all treatments combinations occur with the same frequency, no missing observations), there is also another way to model this, by including a fixed effect for patient in the model.

LM <- lm(y ~ Patient + Tissue * State, data = DF)

The Difference Between the Solutions

Now let's compare their output:

# Mixed model with a random effect for patient
summary(LMM)

# Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
# Formula: y ~ Tissue * State + (1 | Patient)
#    Data: DF
# 
# REML criterion at convergence: -385.8
# 
# Scaled residuals: 
#      Min       1Q   Median       3Q      Max 
# -1.32589 -0.68938  0.01221  0.57604  1.93157 
# 
# Random effects:
#  Groups   Name        Variance  Std.Dev. 
#  Patient  (Intercept) 6.136e-09 7.834e-05
#  Residual             2.413e-09 4.912e-05
# Number of obs: 28, groups:  Patient, 7
# 
# Fixed effects:
#                          Estimate Std. Error         df t value Pr(>|t|)    
# (Intercept)             1.485e-03  3.495e-05  9.428e+00  42.503 4.31e-12 ***
# TissueB                -4.214e-05  2.625e-05  1.800e+01  -1.605    0.126    
# StateHypoxia            1.428e-04  2.625e-05  1.800e+01   5.438 3.63e-05 ***
# TissueB:StateHyperoxia -5.129e-05  3.713e-05  1.800e+01  -1.381    0.184    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Correlation of Fixed Effects:
#             (Intr) TissuB SttHyp
# TissueB     -0.376              
# StateHypoxi -0.376  0.500       
# TssB:SttHyp  0.266 -0.707 -0.707

# Ordinary linear model with a fixed effect for patient
summary(LM)

# Call:
# lm(formula = y ~ Patient + Tissue * State, data = DF)
# 
# Residuals:
#        Min         1Q     Median         3Q        Max 
# -6.729e-05 -2.848e-05  1.143e-06  2.089e-05  9.271e-05 
# 
# Coefficients:
#                          Estimate Std. Error t value Pr(>|t|)    
# (Intercept)             1.412e-03  2.935e-05  48.120  < 2e-16 ***
# Patient2                1.950e-05  3.473e-05   0.561  0.58142    
# Patient3                9.700e-05  3.473e-05   2.793  0.01202 *  
# Patient4                9.700e-05  3.473e-05   2.793  0.01202 *  
# Patient5               -3.150e-05  3.473e-05  -0.907  0.37641    
# Patient6                1.215e-04  3.473e-05   3.498  0.00257 ** 
# Patient7                2.065e-04  3.473e-05   5.946 1.26e-05 ***
# TissueB                -4.214e-05  2.625e-05  -1.605  0.12586    
# StateHypoxia            1.428e-04  2.625e-05   5.438 3.63e-05 ***
# TissueB:StateHyperoxia -5.129e-05  3.713e-05  -1.381  0.18411    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 4.912e-05 on 18 degrees of freedom
# Multiple R-squared:  0.8715,  Adjusted R-squared:  0.8073 
# F-statistic: 13.57 on 9 and 18 DF,  p-value: 2.516e-06

There are quite some lines to scroll through, but the relevant parts to answer your original question are:

  • The $p$-value for the marginal effect of hyperoxia;
  • The $p$-value for the interaction term.

If the latter is significant, then that means there is no overall effect of hyperoxia, but it depends on which tissue type you consider. For both models though, the interaction term is not significant ($p = 0.184$).

This means we can interpret the former as the overall effect of oxygen administration. Apparently, it has a significant effect on the outcome ($p = 3.63 \cdot 10^{-5}$), increasing whatever it is you measured by about $0.00014$ on average.

Which Model is Better?

Here is where I disagree with some of the other suggestions. The fixed effect model is neither the standard way of analyzing such data, nor is it necessarily simpler.

I think that both are valid ways to analyze these data, each with advantages and disadvantages.

Some considerations to decide which to use:

  • Generally speaking, a simpler model is better, because not everyone has an extensive knowledge of statistics. You should not use something complex if it does not have a benefit, because all you're doing is making it harder for readers to interpret your results and/or judge the validity of your approach.
  • That said, mixed models are a fairly standard method of analysis in clinical science, and I would argue the output is easier to interpret:
    • The fixed effects coefficients tab of the mixed model is much smaller.
    • The differences between these specific 7 individuals is unlikely to be of interest.
  • The mixed model allows you to add an interesting piece of information to the conclusion: The estimated variance in the outcome between patients ($\hat{\sigma}^2_\texttt{Patient} = 6.136 \cdot 10 ^ {-9}$). This allows you to easily compare the difference in oxygenation state, or tissue type, to the overall spread between patients.
  • The mixed model has a different interpretation of the intercept.
    • In the ordinary linear model, it is the normoxic, tissue A value of the outcome for patient 1.
    • In the mixed model, it is the normoxic, tissue A value of the outcome for the average patient.
  • Mixed models require a minimum number of groups to consistently estimate a (non-zero) variance component. This is very nicely illustrated here. With $n = 7$ patients, you are on the lower, but safe side. Still, this is good to keep in mind, should you ever work with an even smaller number of patients.
  • If any of the measurements fail for whatever reason, or should be removed (e.g., outlying because it was not stored correctly), the design is no longer balanced, and the ordinary linear model no longer provides fair $p$-values.
  • Even if you don't end up using them now, mixed models are worthwhile learning about for other reasons
    • They not only allow just the intercept, but also slopes to vary among patients, which is often a realistic scenario (e.g., patients react differently to dosage).
    • They do so in a much more efficient way than ordinary linear models, because you don't need $1$ degree of freedom for every individual estimate, but instead, you are estimating a single variance component.
    • Mixed models allow more complex study designs as well, like dependency structures with crossed and nested random effects. There may be better resources out there, but I have an introductory video on that here.

Data

PS: Here is your data in a format you can use to run the example analyses in this and other answers:

DF <- data.frame(
  y = c(0.0014639, 0.0013740, 0.0014739, 0.0015490, 0.0013931, 0.0015456, 0.0015980,       
        0.0013601, 0.0013480, 0.0014101, 0.0015230, 0.0013829, 0.0014864, 0.0015920,       
        0.001512,  0.001638,  0.001745,  0.001585,  0.001498,  0.001655,  0.001764,        
        0.001464,  0.001518,  0.001559,  0.001531,  0.001400,  0.001599,  0.001672),
  Tissue = factor(rep(c("A", "B"), each = 7, times = 2)),
  State  = factor(rep(c("Normoxia", "Hyperoxia"), each = 14), 
                  levels = c("Normoxia", "Hyperoxia")),
  Patient = factor(rep(1:7, times = 4))
)

$^\dagger$: Mixed, multilevel, hierarchical models are some of the names commonly used for this type of model.

$\endgroup$
5
  • $\begingroup$ Thanks a lot for the detailed explanation. One thing I didn't understand is why should I include a fixed effect for patient in the linear model? If I want to see the effect of oxygen administration and whether that difference is affected by tissue type a model like this wouldn't work? LM <-lm(y ~ Tissue * State, data = DF) $\endgroup$
    – user161260
    Commented Jun 19 at 10:21
  • $\begingroup$ One more thing I've noticed and doesn't make sense to me is that if I restructure the data (add first Hyperoxia values and then Normoxia) the results change. TissueB is significant and StateHypoxia. Why is that happening? y = c(0.001512, 0.001638, 0.001745, 0.001585, 0.001498, 0.001655, 0.001764, 0.001464, 0.001518, 0.001559, 0.001531, 0.001400, 0.001599, 0.001672, 0.0014639, 0.0013740, 0.0014739, 0.0015490, 0.0013931, 0.0015456, 0.0015980, 0.0013601, 0.0013480, 0.0014101, 0.0015230, 0.0013829, 0.0014864, 0.0015920) $\endgroup$
    – user161260
    Commented Jun 19 at 11:51
  • 2
    $\begingroup$ The fixed model has to estimate 10 parameters, but there are only 28 observations. Don't we have to worry about overfitting here? $\endgroup$
    – Schmuddi
    Commented Jun 19 at 12:47
  • $\begingroup$ @Schmuddi I only have an informal answer - we are interested in the average effect of tissue and oxygen on the response y and for that we have 7 replicates (patients), which is ok. If we had thousands of patients it would be even better even if there would be thousands of parameters. The effect of each patient on y is indeed poorly estimated but we are not interested in that. If you are familiar with R you can use the simulation code in my answer and replace "7" with some other number. $\endgroup$
    – dariober
    Commented Jun 20 at 8:04
  • $\begingroup$ Great answer- In my experience, though, mixed models are more temperamental to fit than fixed effects ones. See for example stats.stackexchange.com/questions/509892/… $\endgroup$
    – dariober
    Commented Jun 20 at 8:08
0
$\begingroup$

I am ignoring all issues regarding the specific details of the problem that were already addressed above, and just point that (basic) ANOVA assumes that the standard deviations of the two groups are equal, while t-test can handle different standard deviations.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.