I am struggling to find the appropriate statistical test to analyze my data. I hope that my question will be understandable.

I have the following setup:

  • A porcine spine with three vertebral bodies (L1,L2,L3).

  • The spine was scanned on three different imaging modalities (Modality A,B,C)

  • On each of the modalities, different rings of fat were wrapped around the spine resulting in 5 different simulated sizes (size 1 to 5).
  • For each vertebral body of each of the sizes of each modality, I can measure the bone density (BD) as BD.L1, BD.L2, BD.L3

Here the first 10 rows of the table structure with some fictional values for the BD:

my.df <- structure(list(Modality = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L), .Label = c("A", "B", "C"), class = "factor"), 
    Size = structure(c(1L, 1L, 2L, 2L, 3L, 3L, 4L, 4L, 5L, 5L
    ), .Label = c("1", "2", "3", "4", "5"), class = "factor"), 
    Repeat = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
    2L), .Label = c("1", "2", "3"), class = "factor"), BD.L1 = c(1.3, 
    1.5, 2.2, 1.2, 1.8, 1.7, 0.7, 2.3, 2.5, 1.3), BD.L2 = c(1.2, 
    1.7, 1.6, 1.6, 1.1, 1.3, 1, 1.3, 1.2, 1.5), BD.L3 = c(1.6, 
    1, 1.8, 1.2, 1, 1.1, 1.6, 1.5, 1.6, 1.8)), row.names = c(NA, 
10L), class = "data.frame")

I would like to answer the following questions:

  1. Are there significant differences in bone density (BD) measurements among the three modalities for each phantom size?
  2. Are there significant differences in bone density (BD) measurements among the sizes within each modality?

Here the tricky part: for modality A all sizes were scanned twice (2 repeats) while for modalities B and C all sizes were scanned thrice (3 repeats).

Because the data points are very few, I thought to compare the BD measurements for each size not on a per-vertebra basis, but using the BD measurements of all three vertebra together for each modality and size.

Specific questions:

In regards to Analysis 1.) I was thinking about using the Friedman Test. However, I have unequal sample sizes (2 repeats for modality A) vs. (3 repeats for modality B). Which non-parametric test could I use here with unequal sample sizes?

In regards to Analysis 2.): Are the different sizes paired? If I add additional fat rings to the spine is it still considered the same or an independent sample. If independent is it correct to use Kruskal Wallis with Dunn post-hoc test to make comparisons among the five sizes?

I hope that my question is understandable.

Thank you very much!


For reproducibility a dataset representing the full data with fictive values has been added:


df <- data.frame(
  Modality = c(rep("A",30),rep("B",45),rep("C",45)),
  Size = factor(c(rep(rep(1:5,each=2),3),rep(rep(1:5,each=3),6)), levels=c(1,2,3,4,5),ordered=TRUE),
  Repeat = factor(c(rep(1:2,15),rep(rep(1:3,15),2))),
  Level = c(rep(c("L1","L2","L3"),each=10),rep(rep(c("L1","L2","L3"),each=15),2)),
  BD = c(runif(30,1,3),runif(45,2,4),runif(45,3,5))

'data.frame':   120 obs. of  5 variables:
 $ Modality: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ...
 $ Size    : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 1 1 2 2 3 3 4 4 5 5 ...
 $ Repeat  : Factor w/ 3 levels "1","2","3": 1 2 1 2 1 2 1 2 1 2 ...
 $ Level   : Factor w/ 3 levels "L1","L2","L3": 1 1 1 1 1 1 1 1 1 1 ...
 $ BD      : num  2.15 1.45 1.66 2.42 2.64 ...
  • 1
    $\begingroup$ A clarifying question: why are all variables measured categorically? For instance, BD seems like it could be measured continuously. Wouldn't continuous metrics, as appropriate, simplify your model? $\endgroup$
    – user78229
    Commented Apr 10, 2020 at 19:23
  • $\begingroup$ BD is on a continuous scale. The other variables are categorical (Modality and Repeat) and ordinal (Size). $\endgroup$ Commented Apr 11, 2020 at 8:51
  • 1
    $\begingroup$ Ok, my view is that a much more powerful model is possible by retaining quantitative information as quantitative, rather than collapsing them into categories. Potentially valuable information is lost with discretization. $\endgroup$
    – user78229
    Commented Apr 18, 2020 at 18:08
  • $\begingroup$ yes, but it had to be ensured that size was really quantitative in a physical context. $\endgroup$ Commented Apr 21, 2020 at 10:52

2 Answers 2


I am currently exploring the extent of this "aggregation problem" in my own work, for lack of a better term. One answer here suggests that you average the measurements in each modality. This is a form of aggregation and is an easily overlooked violation of assumptions, i.e. that the unobservable error terms of the model are IID, so I believe you are correct in being concerned with it. By averaging them, you are washing out the larger variance of the modality that has fewer measurements. You are also losing degrees of freedom/data points. It is better to include all measurements in the model to account for this difference in variance.

One concern I am noticing is whether phantom size should be continuous or ordinal. The answer to the question is in intro statistics books. Does phantom size have a zero? Considering they are fat layers, then I would think no fat layers could be a choice, so yes. Is twice a phantom size of 1 equivalent to a phantom size of 2? You'll have to answer that, but as far as radii of fat go, the answer is yes. These two yeses indicate a continuous variable. You'd also have to determine what's more reasonable, whether, 1) the area spanned by the fat around the spine affects bone density, or 2) the fat's distance from the center of the spine affects bone density. 1) implies a squared function of radius, whereas 2) implies a linear function. It depends on your research question. If both are true then you can include the linear and squared phantom sizes. Note if the only reason the phantom sizes are ordinal is because you have five phantoms to choose from, this is not a justification for ordinality given your (assumed) research question. You'll simply need to get a ruler and measure the radius of fat from the center of the spine minus the radius of the vertebrum for each phantom. If this paragraph is entirely unreasonable and phantom size should be ordinal, (I really believe it isn't), then I'm wrong and no need to read further, since I assume continuous phantom size. Though if I'm right, at this point I'd suggest changing "phantom size" to "fat layer radius".

Overall I agree with others that you can use a regression model. The response variable is bone density. Three indicator explanatory columns are used for each vertebrum, three more for modality, and one or two more for phantom size depending on if you include a squared term. Both 1) and 2) would be answered by including an interaction term of phantom size*modality (or two if squared term included). If you do not include the squared term, then you have a total of 3+3+1+3*1+3*1 = 13, or 3+3+2+3*2+3*2 = 20 terms if you include the squared term. If you can assume that the measurements from each vertebrum are IID (seems reasonable), then you don't need the 3 indicator variables for vertebra. That would reduce to 3+1+3*1 = 7, or 3+2+3*2 = 11 terms. You'll have to pick the appropriate model considering the IID vertebra question, linear vs squared question, and the number of data points you have. The rule of thumb I've learned is 10 to 1 data points to predictors. Fewer data points could still work depending on the effect size.

Note the above does not require any non-parametric methods. Non-parametric methods should only be used if your residuals do not end up following a normal distribution. The other answers seem to give good code for generating the model, just make sure phantom size is continuous. I see no reason to perform partial f tests, though. Significant interaction terms in the full model suffice in answering your two (equivalent) questions.

  • $\begingroup$ Thank you for your helpful comment. Do I understand correctly (I am a physician) that the model with interaction terms (lm(formula = BMD ~ Modality * Size + Vertebra) would be appropriate to answer both questions? $\endgroup$ Commented Apr 16, 2020 at 9:32
  • $\begingroup$ Yes. (Modality x Size) is equivalent to (Size x Modality), which appears to be your 1) and 2) questions. Again, this assumes Size is continuous, am I correct that this is a reasonable assumption? Since Modality is categorical, it would be split into 3 indicator columns, then each would be multiplied by Size to create 3 new columns. The coefficients from each of these columns should approximately match the slopes that you'd get by fitting a line in the 3 plots provided in your follow-up post. Let me know if they don't match. $\endgroup$
    – Rob G
    Commented Apr 17, 2020 at 4:20
  • $\begingroup$ Is there a way to calculate effect size of the (significant) factors of the multiple linear regression model? $\endgroup$ Commented Apr 21, 2020 at 10:52
  • $\begingroup$ I am actually still in doubt if the sample size of this experiment is enough to perform multiple linear regression because I think that a linear model is based on parametric assumptions. Would a bootstrap hypothesis testing be an alternative or is that not possible because we are looking at 3 groups? Thank you $\endgroup$ Commented May 5, 2020 at 6:59

Answer about analysis 1

If you want to perform a Friedman test, you have to use 15 groups of 3 observations (one group for each combination of spine and fat size, and one observation for modality), so that the samples are paired. Repeated measurements must be reduced to one single observation by averaging them.

The alternative is to use not a non-parametric test, but a linear model, which seems a confortable solution, as it helps with both your problems at one time.

Answer about analysis 2

If you average repeated measurements, you can create five paired samples of 9 observations (same averages as before, but this time divided in samples according to the size). Also, it seems to me that you should avoid Kruskal Wallis test, as it assumes independence between the samples, so it is less powerful than tests that take pairing into account. On the other hand, Kruskal Wallis doesn't requires pairing, so you can use it without averaging repeated measurements; I wouldn't recommend going down this way though, blocking for nuisance factors is generally more effective than having a not too larger sample.

So, you could use a Friedman test again. However, I don't understand why you would even be interested in a post-hoc test between sizes, how comes that a specific size of fat around the bone would be specially interesting? Wouldn't you be satisfied by learning that fat wraps have an effect overall?

In fact, as the size is clearly an ordinal variable (so, basically a numeric one), I would not go for a Friedman test either, but for a linear model. This seems by any mean the best option: it lets you use both blocking factors and repeated measurements, and it allows better insight into data. If you are worried about scientific correctness, just check the diagnostic plots before looking at the tests results: if the effect of the size is not linear, add non-linear effects; if the residuals have variance depending on factors... well, let's hope not, because in that case you'll have to use a different model.

Code for the linear model

Following the discussion on the comments here there is the code you could use for doing analysis n. 2:

# developing the data frame to use it for linear modeling
BD= c(as.matrix(my.df[, 4:6]))
Vertebra= paste('L', gl(3, nrow(my.df)), sep= '')
df= cbind(my.df[, 1:3], Vertebra, BD)
df$Size= as.numeric(as.character(df$Size))

# linear models
m_wo_mod= lm(DB~Size+Vertebra, data= df)
m_wo_size= lm(DB~Modality+Vertebra, data= df)
m_complete= lm(DB~Modality*Size+Vertebra, data= df)

# diagnostic plots
plot(df$Size, m_wo_mod$residuals, col= df$Modality)

# test question 1
anova(m_wo_mod, m_complete)
# test question 2
anova(m_wo_size, m_complete)
  • $\begingroup$ In regards to Answer 1: Do I have enough data points to perform a linear regression? I am interested to test if there are significant differences in BD among the three modalities. Would it be appropriate to take the average BD values from the repeats for all three vertebrae and 5 sizes together (15 values per modality) and compare the modalities using Friedman (then they would be paired)? Or is that inappropriate because I average across different vertebrae? $\endgroup$ Commented Apr 11, 2020 at 9:40
  • $\begingroup$ In regards to Answer 2: Your idea to assess if size has an impact overall is interesting. However, I think this has to be done for each modality separately. This results in the same question as for Answer 1: Is my sample size large enough for a linear regression? $\endgroup$ Commented Apr 11, 2020 at 9:42
  • $\begingroup$ linear models have a hard lower limit for sample size, but it is incredibly low. in practice, when we have already collected the sample, we want the most powerful method, which is the one effective with the least data points. the linear model is more powerful than any non-parametric method, and this holds even more true in your case, as it doesn't require any loss of information from data. $\endgroup$
    – carlo
    Commented Apr 11, 2020 at 12:56
  • 1
    $\begingroup$ you can definitely test effect of size for each modality separately, but you have to divide your sample, which is not that big. I would recommend that only on a linear model, after assessing overall effect of size, when estimating an interaction between size and modality. ask if you need the code. $\endgroup$
    – carlo
    Commented Apr 11, 2020 at 13:14
  • 1
    $\begingroup$ Thank you. Did you see the full dataframe structure I have updated at the end of the question? This is the complete structure of the data with fictive data $\endgroup$ Commented Apr 12, 2020 at 20:40

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.