Hello Imen Zribi. Readers will be better able to advise you if you provide a bit more info:

Thanks for clarifying.

Article Practice Statistics Notes Analysis of continuous data from s...

Bruce Weaver Thank you for your answer I updated my question. I added the mean±SD. Can you please have a look at them?

You don't mention which software you are using. If using R, the default for

oneway.test

the procedure for oneway ANOVA, is to use Welch's (1951) procedure that does not make this assumption. But this may not be your best approach. What is your response variable? Is it bounded by zero? Also, is there a reason your sample size is so small? Is the larger sd for condition 1 due to a single big outlier? I guess with the sample size it is tough to tell.

B. L. Welch (1951). On the comparison of several mean values: an alternative approach. Biometrika, 38, 330–336. doi: 10.2307/2332579.

As Daniel Wright noted, that SD for treatment 1 is an eye-opener. In hindsight, I should have just asked for the raw data. Can you share it with us? And to repeat the first question in my earlier post (echoed by Daniel), what is the dependent variable? Thanks.

Another question is whether you are just interested in whether there are differences among treatments, or interested in where those differences are? Neither way be practical to look at given the sample size, but even more so the later.

Imen Zribi ,

Thanks you for supplying the summary of the data.

I agree that with the standard deviations reported, even with only a few observations for each group, that assuming homogeneity may not be the best approach.

There are quite a few options to analyze a one-way design with heterogeneity in groups. There's a very useful list here: https://stats.stackexchange.com/questions/91872/alternatives-to-one-way-anova-for-heteroskedastic-data/91881#91881 .

##########################

Just for fun, I'm going to assume that the data in each group are symmetric about the mean, and conduct a few analyses.

But in reality, consider Bruce Weaver's request to look at the actual raw data, and Daniel Wright's question about if the large variance in Group 1 is caused by a single outlying observation.

* * *

In R, which could be run here without installing software: https://rdrr.io/snippets/ .

A standard anova with White adjustment to account for heteroscedasticity. Welch's anova. Kruskal-Wallis test with Dunn test post-hoc. Fisher-Pitman permutation test.

And then I transformed the dependent variable (YTrans). I used Elving inverse scores transformation, which here makes the data better meet anova assumptions. A standard anova with White adjustment with pairwise comparison of estimated marginal means.

There's a lot of discussion about how useful or desirable transforming data is.

The upshot is that --- for my hypothetical data --- I have some evidence to support the idea that group 1 is different in response from group 3. But different methods of analysis lead to different interpretations.

Finally, I used Gamma regression. See ( https://en.wikipedia.org/wiki/Gamma_distribution ). Here, the results suggest more strongly that Treatment 1 is different from the others.

######### Create data ##########

Treatment1 = c(20.90, 20.90+11.63, 20.90-11.63)

Treatment2 = c(7.29, 7.29+1.37, 7.29-1.37)

Treatment3 = c(4.17, 4.17+1.67, 4.17-1.67)

Treatment4 = c(7.54, 7.54+1.74, 7.54-1.74)

Y = c(Treatment1, Treatment2, Treatment3, Treatment4)

Treatment = factor(rep(c("1","2","3","4"), each=3))

library(rcompanion)

YTrans = blom(Y)

Data = data.frame(Treatment, Y, YTrans)

########## Summarize and plot data ##########

library(FSA)

Summarize(Y ~ Treatment)

# Treatment n mean sd min Q1 median Q3 max

# 1 1 3 20.90 11.63 9.27 15.085 20.90 26.715 32.53

# 2 2 3 7.29 1.37 5.92 6.605 7.29 7.975 8.66

# 3 3 3 4.17 1.67 2.50 3.335 4.17 5.005 5.84

# 4 4 3 7.54 1.74 5.80 6.670 7.54 8.410 9.28

library(ggplot2)

ggplot(Data, aes(x=Treatment, y=Y)) +

geom_dotplot(binaxis='y', stackdir='center')

########### Analyze un-transformed data ##########

model = lm (Y ~ Treatment)

library(car)

Anova(model, white.adjust=TRUE)

# Analysis of Deviance Table (Type II tests)

# Df F Pr(>F)

# Treatment 3 2.779 0.1101

# Residuals 8

oneway.test(Y ~ Treatment)

# One-way analysis of means (not assuming equal variances)

# F = 3.1755, num df = 3.0000, denom df = 4.2639, p-value = 0.1395

kruskal.test(Y ~ Treatment)

library(FSA)

dunnTest(Y ~ Treatment)

# Kruskal-Wallis rank sum test

# Kruskal-Wallis chi-squared = 8.0256, df = 3, p-value = 0.04548

#

# Comparison Z P.unadj P.adj

# 1 1 - 2 1.4719601 0.141031641 0.5641266

# 2 1 - 3 2.8306926 0.004644734 0.0278684

# 3 2 - 3 1.3587324 0.174231388 0.5226942

# 4 1 - 4 1.3587324 0.174231388 0.3484628

# 5 2 - 4 -0.1132277 0.909850033 0.9098500

# 6 3 - 4 -1.4719601 0.141031641 0.7051582

library(coin)

oneway_test(Y ~ Treatment)

# Asymptotic K-Sample Fisher-Pitman Permutation Test

# chi-squared = 6.9911, df = 3, p-value = 0.07218

########### Analyze transformed data ##########

ggplot(Data, aes(x=Treatment, y=YTrans)) +

geom_dotplot(binaxis='y', stackdir='center')

model2 = lm (YTrans ~ Treatment)

library(car)

Anova(model2, white.adjust=TRUE)

library(emmeans)

marginal = emmeans(model2, ~ Treatment)

pairs(marginal)

# Analysis of Deviance Table (Type II tests)

# Df F Pr(>F)

# Treatment 3 5.2566 0.02697 *

# Residuals 8

# contrast estimate SE df t.ratio p.value

# 1 - 2 1.1357 0.475 8 2.392 0.1560

# 1 - 3 2.2029 0.475 8 4.640 0.0072

# 1 - 4 1.0672 0.475 8 2.248 0.1901

# 2 - 3 1.0672 0.475 8 2.248 0.1901

# 2 - 4 -0.0685 0.475 8 -0.144 0.9988

# 3 - 4 -1.1357 0.475 8 -2.392 0.1560

#

# P value adjustment: tukey method for comparing a family of 4 estimate

########### Gamma regression ##########

model3 = glm(Y ~ Treatment, family=Gamma())

Anova(model3)

marginal3 = emmeans(model3, ~ Treatment)

summary(marginal3, type="response")

pairs(marginal3)

# Analysis of Deviance Table (Type II tests)

# LR Chisq Df Pr(>Chisq)

# Treatment 31.189 3 7.755e-07 ***

#

# Treatment response SE df asymp.LCL asymp.UCL

# 1 20.90 4.51 Inf 14.69 36.21

# 2 7.29 1.57 Inf 5.12 12.63

# 3 4.17 0.90 Inf 2.93 7.23

# 4 7.54 1.63 Inf 5.30 13.06

#

# contrast estimate SE df z.ratio p.value

# 1 - 2 -0.08933 0.0313 Inf -2.850 0.0227

# 1 - 3 -0.19196 0.0528 Inf -3.638 0.0016

# 1 - 4 -0.08478 0.0304 Inf -2.787 0.0273

# 2 - 3 -0.10263 0.0596 Inf -1.722 0.3122

# 2 - 4 0.00455 0.0412 Inf 0.110 0.9995

# 3 - 4 0.10718 0.0591 Inf 1.813 0.2672

# Note: contrasts are still on the inverse scale

# P value adjustment: tukey method for comparing a family of 4 estimates

Let's all remember that it is the population variances, not the sample variances, that are assumed to be homogeneous. And with n = 3 in each group, one could easily see considerable variance in the sample variances, even if all population variances were exactly equal. Rather than relying on what such small samples can tell you, you ought to rely on what you know about the variable from previous experience, or from the literature. What are the 4 treatments? Have you or others examined them before? Is the treatment with the large SD expected to increase variability of the scores?

Returning to a question Daniel asked earlier, why do you have only n = 3 per group? Having such a small sample size per group does not really lend itself to any kind of serious analysis.

Thank you Bruce Weaver Salvatore S. Mangiafico Daniel Wright for answering my question. My response variable is non-structural carbohydrates concentration in below-ground parts. Indeed the larger SD for condition 1 is due to the presence of a big outlier. Here is the raw data :

Treatment 1 : 34.39 ; 13.75 ; 14.75

Treatment 2 : 6.12 ; 8.80 ; 6.97

Treatment 3 : 6.08 ; 3.48 ; 2.96

Treatment 4 : 9.55 ; 6.45 ; 6.63

PS: I am using SPSS to run the statistical analysis. Thank you again for the advice and recommendations.

IIRC, concentrations are often log-transformed. Do people in your area of research log-transform non-structural carbohydrate concentration?

Imen Zribi , are you interested if there are differences among treatments (and these treatments in particular or some population of treatment) or are you interested in differences between particular pairs of these four treatments (or other contrasts, and if so which ones)?

Log-transformation sounds like a reasonable approach. The transformed data appears fairly reasonable for ordinary least squares analysis (including common anova). Otherwise, the Kruskal-Wallis or Gamma regression approach may be useful.

I think before advising on a statistical procedure it is worth knowing why such a small sample is needed. Further, given the size of the SDs, assuming something like that was used in the power analysis, the minimum mean difference to detect would need to be very large.

Daniel Wright Salvatore S. Mangiafico Bruce Weaver Thank you again for your time and help.

@Bruce Weaver Yes indeed some studies have mentioned that in order to meet the homogeneity assumption data were log10(x + 1) or square root transformed.

Daniel Wright I am interested to see the difference between these treatments. I guess I have to clarify that Treatment 1 is the control and treatment 2, 3 and 4 are actually the same treatment but with different gradients. Treatment 2: level1 /treatment 3: level2 / treatment 4: level 3.

Imen Zribi , if the treatments could be thought of as None, Low, Medium, and High treatments, you might want to consider this in the analysis. If you are using a common anova approach, there is a way to fit polynomial contrasts so that you can see if there is e.g. a linear trend or quadratic trend across treatments. In your case, there doesn't seem to be any kind of trend within the three non-control treatments. But considering the four treatments together, there is a suggestion that there is a quadratic response. That is, there is a hint that the level of Treatment 3 may be associated with the lowest level of the response variable. This kind of response can be common in agricultural research. Like you apply a nutrient and that increases yield, but at too-high of a rate it decreases yield. ( Example: https://www.researchgate.net/profile/Mehmood-Ali-Noor/publication/335541725/figure/fig3/AS:798593015373825@1567410820912/Typical-nitrogen-response-curve_W640.jpg ).

Unfortunately you probably don't have sufficient observations to confidently intuit what the response curve might look like.

BTW, there are also nonparametric tests that treat the groups as ordered. See, e.g. Cuzick test and Jonckheere-Terpstra test.

Salvatore S. Mangiafico , if a different quadratic is used for group, you are out of df! I think the questioner has to answer the question I posed above about why this sample size was chosen before we should suggest more.

Hi, Daniel Wright .

For a quadratic contrast across groups, usually there would be a linear contrast which would take one degree of freedom, and then a quadratic contrast contrast that would take one degree of freedom. So I don't think there's a problem with burning up the degrees of freedom for such an approach... These are sometimes called "single degree-of-freedom contrasts". ...

I agree with the issue of sample size.... But, it's not unusual to have three experimental units per treatment in agricultural experiments. Granted, this is placing a lot of faith for consistent results from an inherently variable biological system. ... I guess we just inherent optimism from farmers. If we told them the probability of making a reasonable profit from farming year to year, 90% of them would probably drop out of the profession....

Salvatore S. Mangiafico , what are you assuming for the intercept? I was assuming to estimate it too.

On the sample size, the observed SD suggest it is a problem here. Here is the example Fisher gave with n=2.

\begin{quote}

The apparent paradox may be explained by pointing out that the difference of 10 tons between the experiments indicates the existence of uncontrolled circumstances so influential that in both cases the apparent benefit may be due to chance, whereas in the former case the relatively close agreement of the results suggests that the uncontrolled factors are not so very influential. Much of the advantage of further replication lies in the fact that with duplicates our estimate of the importance of the uncontrolled factors is so extremely hazardous. \\ \phantom{d} \hfill \citet[p.~112]{Fisher1925}

\end{quote}

Daniel Wright , I don't disagree with you about the sample size. I'm not sure what else to say about it. ... Though on the other hand, depending on the field of study and the purpose of the analysis, data on this scale aren't that unusual. And note that these particular data will produce the hallmarks of results we'd be happy with: p's < 0.01, r-square of about 0.8, non-overlapping confidence intervals of treatments. It's difficult to complain. But that's a whole other discussion.

The kind of polynomial contrasts I'm thinking of are detailed here: https://www.ndsu.edu/faculty/horsley/Polycnst.pdf .

The linear contrast uses one degree of freedom and the quadratic contrast uses one. But traditionally these aren't subtracted from the df for the error term, but there's acknowledgement that for four treatments, you shouldn't have more than three single-degree-of-freedom contrasts. ... But I'm not sure I agree that it's necessary to limit the number of post-hoc contrasts, especially if p-values are adjusted as in multiple comparisons.

So, the approach doesn't estimate the actual polynomial model that would fit the data. Though one could do that afterwards.

Below I have the R code and results for the data here analyzed with polynomial contrasts. For the R users, I prefer analyzing these kinds of contrasts with the emmeans or multcomp packages, but the aov approach is maybe more traditional in that it limits the df used in these contrasts.

________

library(ggplot2)

library(emmeans)

library(multcomp)

### Data, with log transformation

Treatment1 = c(34.39 , 13.75 , 14.75)

Treatment2 = c( 6.12 , 8.80 , 6.97)

Treatment3 = c( 6.08 , 3.48 , 2.96)

Treatment4 = c( 9.55 , 6.45 , 6.63)

Y = c(Treatment1, Treatment2, Treatment3, Treatment4)

Treatment = factor(rep(c("1","2","3","4"), each=3))

logY = log(Y)

Data = data.frame(Treatment, Y, logY)

### Plot

ggplot(Data, aes(x=Treatment, y=logY)) +

geom_dotplot(binaxis='y', stackdir='center')

### Define contrasts

linear = c(-3, -1, 1, 3)

quadratic = c( 1, -1, -1, 1)

Matriz = cbind(linear, quadratic)

contrasts(Data$Treatment) = Matriz

CList = list("Linear" = 1,

"Quadratic" = 2)

### Model, contrasts, and post-hoc

model = aov(logY ~ Treatment, data = Data)

summary(model, split=list(Treatment=CList))

marginal = emmeans(model, ~ Treatment)

# summary(marginal)

# pairs(marginal)

multcomp::cld(marginal, Letters=letters)

### Plot with means and confidence intervals

Sum = as.data.frame(summary(marginal))

ggplot(Sum, aes(x = Treatment, y = emmean)) +

geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),

width = 0.05, size = 0.5) +

geom_point(shape = 15, size = 4) +

theme_bw() +

theme(axis.title = element_text(face = "bold")) +

ylab("Marginal mean\n") +

xlab("\nTreatment")

### r-square

RSquare = summary(lm(logY ~ Treatment, data=Data))$r.squared

round(RSquare,3)

### Results

Df Sum Sq Mean Sq F value Pr(>F)

Treatment 3 3.794 1.2648 10.44 0.00386 **

Treatment: Linear 1 1.770 1.7696 14.60 0.00508 **

Treatment: Quadratic 1 1.918 1.9177 15.82 0.00407 **

Residuals 8 0.970 0.1212

---

Treatment emmean SE df lower.CL upper.CL .group

3 1.38 0.201 8 0.916 1.84 a

2 1.98 0.201 8 1.512 2.44 a

4 2.00 0.201 8 1.541 2.47 a

1 2.95 0.201 8 2.487 3.41 b

Confidence level used: 0.95

P value adjustment: tukey method for comparing a family of 4 estimates

significance level used: alpha = 0.05

NOTE: Compact letter displays can be misleading

because they show NON-findings rather than findings.

RSquare

0.796