Hi Ioana,

What you probably should use for that data is Linear Mixed Models (lmer) package.

Here is a link to an example that seems very similar to your data:

http://stats.stackexchange.com/questions/58745/using-lmer-for-repeated-measures-linear-mixed-effect-model.

And here is a link to a more general introduction to what linear mixed models are. They are also available in SPSS.

http://www.bodowinter.com/tutorial/bw_LME_tutorial2.pdf

I can highly recommend the book by Faraway that Dasapta mentioned above. The method described here is described in Faraway's second book:

http://www.amazon.com/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X

Another excellent resource is

"Applied Longitudinal Data Analysis: Modeling Change and Event Occurrence"

by Judith D. Singer and John B. Willett , which is written more for people who study psychology, so it is less mathematically heavy than the faraway books. The associated website has all kinds of nice info, including an excellent video of a seminar of the authors:

http://www.ats.ucla.edu/stat/examples/alda.htm

But  unless you have time to read two books, the two links in the beginning of my post offers a faster choice :)

kind regards

When there is problem of normality, you can try a data transformation to "normalized" them.

- log10 log2

- 1/x

- sqrt(x)

etc.

Data transformation are often used in different disciplines (microarray, data on population...) and may help you to keep the most powerful statistic test. http://udel.edu/~mcdonald/stattransform.html

Linear Mixed Models is a good idea.

You could also look at multi way nonparametric anova :

http://onlinelibrary.wiley.com/doi/10.1111/j.1442-9993.2001.01070.pp.x/full

http://stats.stackexchange.com/questions/3634/multi-way-nonparametric-anova

When all subjects are measured at all 4 time points then a mixel model is not required. Keep things simple!

What is your response variable? Is there a probably a meaningful transformation applicable? Or are these counts (->Poisson regression), proprotions (->beta-regression), concentrations (-> gamma model)?

Non-linear models will only be useful when the estimated parameters are non-linear in the model. Since you have a categorical predictor I have no idea how this should be possible (at least in a reasonable way)...

You may have a look in R at ?friedman.test maybe this solves your problem.

While it is true that as the data is balanced (we have same number of measurements for all time points) the main benefit of mixed models is lost, the lmm is more robust (when compared to repeated measures measures ANOVA) in other ways as well: It handles heteroscedasticity (different amount of variance between groups) as well as allow you to specify autoregressive correlation matrix for the residuals. It, however, does NOT handle the non-normality which i forgot to mention, so some kind transformation (e.g., logarithmic) would be necessary. So the analysis does get a bit complicated.

Friedman test seems like an excellent idea, though being a non-parametric test, it might be less powerful than the parametric test,

best regards

Hi

Thank you all for your answers. My outcome variables are physiological outcomes. One is the high frequency (HF) component of heart rate variability which was normalized (HFn), following standard procedure (divided by total power). The other one is heart rate. I am interested in the effects of Time (did these variables modified across subjects regardless of their experimental condition) and the Group by Time interaction. More precisely, the first 2 phases are 2 baseline phases (eyes close and watching a neutral film clip), then in phase 3 subjects got one of three types of emotion regulation instructions (this would be my categorical variable), and in phase 4 they were subjected to an emotional induction task. So I would be particularly interested in an interaction effect in phase 4.

My predictions are: (1) there was a significant effect of Time particularly from phases 1, 2 and 3 to Phase 4 (I expect no differences between phases 1, 2, and 3); (2) There is a significant interaction effect in Phase 4.

One problem I have regards outliers. I tried both the univariate approach around the MAD (median absolute deviation) recommended by https://dipot.ulb.ac.be/dspace/bitstream/2013/139499/1/Leys_MAD_final-libre.pdf and the Mahalanobis distance and there are a lot of values that would qualify as outliers. I don't think they should be just excluded, because maybe they are an important effect of the experimental design. Anyway this is the reason I was looking for a more robust form of ANOVA.

Best and thank you again for taking the time to help me

Ioana

Agree for Friedman test. Best solution for few time involved.

But it's a non parametric test (rank test).

So if you plan to use it, you can forget about the data tranformation (which could be tedious when you don't know where to go, even if Jochen Wilhelm gave you some leads for that... ;-) ).

By the way, I look at your table and I think you should reorganize your data by :

one column for one factor + the final response columns.

Could look like this :

Gr    ?Repetition?     HFn      HR     

3                 1             ,6524     77,73

3                 1             ,7846      83,80

...                ...               ...            ...

(If "Repetition" can be considered as a factor and if I understand well your data)

Be careful also with the name of your group, your group names have no reason to be numeric and have to be treated as character.

> is.numeric(3)

[1] TRUE

> is.numeric("3")

[1] FALSE

> is.character("3")

[1] TRUE

> is.character("emotion")

[1] TRUE

Good Luck !

Hi, 

if you want to model your data as a Linear Mixed Model, 

Package ‘robustlmm’

can be a viable option. It provides some support for the NON NORMAL distributed data.

best, 

http://cran.r-project.org/web/packages/robustlmm/index.html

Hej,

try it as simple as possible, the use of mixed models seems to be like shooting with a sledgehammer to crack especially as you have only 4 measurements (recommended  are something higher than 10/15 with more than 4/6 levels within one random factor and than Faraway is really the best choice to read and A. Zuur for more not mathematical interested people). The first thing I would do making plots of the changes you have detect/not detect, to get a feeling which model might fit best. Even from a plot you will get a feeling what kind of effect is available (linear /non linear) in general or for different subgroups.  Out-layers will tell you something about the natural variability  or about errors in you measurements/design. They are valuable, so before transforming your data or deleting some of them, take care!

Sven

Hi

Thanks, you are all very kind, sorry if I wasn't clear. My dependent variables are HR and HFn (they are both indici of cardiac activity). HFn1 means HFn at time 1 (first baseline), HFn2 means HFn at time 2 (second baseline), HFn3 means HFn at time 3 (each group gets an emotion regulation strategy), HFn4 means HFn at time 4.

Group is a categorical variable, it's with numbers because that's how SPSS does it. 1, 2, 3 each denote an emotion regulation strategy.

I am interested in the effect of Time ( I expect lower HFn and higher HR at time 4, across groups as compared to times 1, 2 and 3). I am also interested in the Time*Group interaction at time 4.

Best regards

Ioana

Ioana, I just tried to analyze the data you gave. The responses (HFn and HR) should be log-ed. The residuals look fine, no proplems, no outliers. There is a clear effect of "Time", but no relevant interaction of Gr x Time for neither of the responses.

Hi Jochen

That's incredibly kind of you, thank you very much! I also didn't get an interaction effect doing the analysis in SPSS (unless I used baseline values as covariates), but I was worried about outliers.

Just one question: why do you say they should be log-ed? I mean, how do you realize that? I ask this because HFn is considered normalized as it is, it's a standard procedure used for this kind of data in which power in the high frequency domain of heart rate variability (HF) is divided by total power, maybe not a proper normalization.

Best

Ioana

Ioana, normalization does not mean that the variable (or the residuals) will be normal distributed.

I realiezd it after lookiong at the residual plots. It was clearly visible that the variance depends on the predicted value and the normal-QQ plot showed a curved pattern. Both features indicate that the variable is roughly log-normal, and log-ing the variable will make it normal distributed. After using the log-ed values, the residual ploats look perfectly fine (as attached above).

Dear Ioana,

Dear Jochen,

I also look at the log transformation, for sure it improve the data normality but not enough to pass the Shapiro-Wilk test. Any other transformation did it.

Code for multiple Shapiro-Wilk test (results enclosed)

> tapply(AER2$HFn,interaction(AER2$Gr,AER2$Time,drop=TRUE),shapiro.test)

> tapply(log(AER2$HFn),interaction(AER2$Gr,AER2$Time,drop=TRUE),shapiro.test)

Your logged data (almost normalized) didn't also pass the Bartlett test.

All this thing told us that we should stick with non-parametric test. As there is repetition in your data, you can't use Friedman test...

So best solution, for me, is to rank your results and compute a standard anova with them.

>AER2$HFn bartlett.test(HFn~Time*Gr,data=AER2)

> bartlett.test(HFn~Gr*Time,data=AER2)

> anova(lm(HFn~Time*Gr,data=AER2))

Results (enclosed in figures) look like those obtained by Jochen.

So here, I think there is not one way to compute your data. You should choose your way and find your compromise between robustness and power.

Good Luck

Cyril

Dear Cyrill,

thanks for your critical review. I wonder why a Shapiro-test would be helpful to decide if a normal approximation is good enough. Each response (HR and HFn) has 408 values. The power to detect a statistically significant deviation from a normal distribution is very high. In other words: it is very likely to get "significant results" even for irrelevant deviations. The test does not (help to) judge what kind and amount of deviation is relevant. Therefore I do not see a reason why such a test should be performed at all.

Apart from this, the Shapiro-test on the residuals for the HR data is not significant (what I find quite unexpected, though):

> m = lm(log(values)~Time*Gr,xx,subset=Part=="HR")

> shapiro.test(resid(m))

Shapiro-Wilk normality test

data: resid(m)

W = 0.9961, p-value = 0.4068

Only for the HFn data is gabe a tiny p-value (

Dear Jochen,

Firsty I would like to thank you for your answer and for sharing your script. I work more with small population, that would explain some of my automatic reflexes in the tests I perform.

Anyway, I agree that the use of Shapiro test on big population would be likely to detect significant deviation. But I computed the Shapiro test on each 12 statistic groups, which are composed of 30 values and Shapiro test is quite efficient for this size of population (please find attach txt).

Why do I perform a Shapiro test ? : because in a first phase we don't know that 'Gr' is not significant and a group correspond with one batch of manipulation (which should be some how normal). And I would like to know if the plotted p-value is the "'right one" or if it is a kind a lower than the "real" p-value. But maybe I'm wrong here.

Why do I perform a non-parametric test ?: I said because of non-normality, but that was mainly because of the non homogeneity of variance (please see attach pictures). Both increase the Type I error rate, so it was a simple precautionary principle. (Even if variance are not "extremely" different and the number of group is small...)

Do I change really the question ?: I would say no, because Ioana finished her question by "I am mostly interested in the interaction effect between Time and Group" and didn't limit the effect detection on mean statistics. And finally, on that question, we get the same answer.

To conclude, I'm sure that you are more expert than me on the subject and Ioana should follow your instructions. I just explain my choices.

I would like to thank you, because I know more now on prerequisite evaluation.

Furthermore you seems to have George Box with you so...

“To make the preliminary test on variances is rather like putting to sea in a rowing boat to find out whether conditions are sufficiently calm for an ocean liner to leave port !” George Box, 1953

Have a nice day !

Cyril

Cyril, thank you for your coprehensive response. I would not consider myself being "are more expert than [you] on the subject". I just try to understand the logic behind the statistical procedures and share my insights (and sometimes my convictions). I find this formum a great opportunity to discuss such things and learn from each other.

I may share a comment on the difference that you analyszed the 12 groups separately and I analyzed the residuals of the whole model comprising (2x) 6 of these groups together:

I agree that the Shapiro-test is appropriate and is powerful for such "small samples" (and the power is still greater for larger samples, though). The "problem" to perform the test individually and independently on all subsets increases the FWER, so you are more likely to get one (or more) false-positives than you hoped to get. An adjustment would possibly be required then. This can be ignored when all data is analyzed together and only one single test is performed. Further, having a "statistically significant" deviation from the normal distribution in some groups but not in others is then a very unpleasant result. How would you reasonably interpret it? Why should the error distributions depend on the groups? If so, it would question if these groups are comparable at all (or if one would compare apples and peaches). The experimental design should assure that the errors in all groups are "of the same kind" (I hope you get the point). So either the experiment is flawed or the model is inappropriate. In both cases we would (should) not bother to analyze any data at all. Given the experiment is ok and the model is (roughly) appropriate, the combined analysis of all residuals together would be even more powerful to detect deviations from the normal distribution, and false-positive alarms among sub-groups would be avoided.

And finally, on that question, we get the same answer. - Not exactly ;) We both may find a non-significant interaction, but this is a minor part of the story. There is still the interesting scientific question of the strength (size?) of the interaction. It can still be that the uncertainty is so large to obscure a relevant effect size (if there was one). I can not judge what a relevant effect size would be here, this has to be judged by Ioana as the expert in the scientific field. The analysis should provide Ioana an estimate together with a "range of not-too-implausible values" for this estimate =confidence interval), given the data she has at hand. This can be sensibly interpreted. A p-value alone can not be interpreted at all. It may serve to reject a null hypothesis, but this does not gain insight (this takes a while to understand, especially when thousands of scientists would strongly disagree to this sentence! It's a rather long story through the logic of hypothesis tests and you will have to carefully read Neyman's and Fisher's original contributions [especially textbooks are often completely flawed in this topic!]).

To give an example (made-up, with many bolts and nuts, but it clarifies the principle): the survival times of cancer patients getting an anti-cancer drug A or B is investigated. The study gives a p-value for drug A of 0.001, for drug B of 0.23. Which drug would you take if you had cancer? Drug A, becuse it was "significant"? - I say: you can not decide. Too much information is missing. If I would tell you that drug A increased the survival by about 2 days, very precisely (and thus "significantly") but drug B increased the survival on average by 1 year, but for some patients by 10 or 20 years, whereas on others it had no effect, I'd guess that you would prefer drug B now. You can only sensibly interpret the results when you look at the effects, but not when you look (only) at p-values.

And this brings me to my skecpticism about the rank-analyses (tha do have their place, though). What is the effect? How do I interpret it?

And yes: the quote from Goerge Box is a deep truth - I surely agree here :)

Have a nice day, too!

Dear Cyril

Dear Jochen

Thank you both for these very interesting comments and in particular for the scripts.  You are very kind!

I have mainly worked in SPSS before and non-parametric analysis is limited there, but it was my understanding that was particularly suitable for small data, which I didn't consider to be my case.  While this is by no means a big study, it is reasonably well-powered both in comparison with the literature and in regards to the type of effect sizes I was looking for.

I was indeed mostly interested in the interaction effect, but that was an exploratory hypothesis.

Best  regards

Ioana

@Ioana : Youre welcome, but I thank you also... that was great to discuss and learn with Jochen using your dataset. :-)

@Jochen : Because you finished with these questions (and because I agree with the rest) : "And this brings me to my skecpticism about the rank-analyses (tha do have their place, though). What is the effect? How do I interpret it?"

I would say : we should perform a rank analysis when we trust more the medians than the means. (see attach pictures - Same example as you with the two formula A and B, survival times in month. I tried to show : the efficiency on median test, the need of precision even with rank test - with a small population because it's my hobbyhorse ;-) ). Regarding median : A>B (significantly)... the survival of 45 months in B is singular and fortunate for this formulation.

Best regards

Cyril

Cyril, I think it is a misconception that the Wilcoxon test is about medians. Just to demonstrate this by example I attached another script.  Following the Popperian logic this example (although possibly arbitrary) falsifies the premise that the Wilcoxon-test (and derivatives) are about medians.

The test is about the entire distribution, with a particular sensitivity to location shifts. If and only if the shapes of the distributions are identical (except the "location"), then and only then the test is actually a test about the null hypothesis that the medians (or means or any other measure of the location!) are the same.

To specifically test the medians this test makes a rather strong assumption. It does not require the assumption of the (maximum-entropy!!) condition of normal distribution but instead it makes the assumption that the entire distributions are identical except the location (shift). This assumption is generally unreasonable in biological settings, especially when it is already obvious from the data that the variances are not identical. The test still gives you an answer to the question if you can expect one group giving higher values than the other, but this is not (neccesarily) related to the means or medians (or any other location measure) - as demonstrated by example.

The means are actually "expected values", they express an expectation. A mean difference is an expectation about a difference (between groups). I find it straight forward to ask a question like: "Is it reasonable (given the data) to expect a difference between the groups?" or, better phrased: "What amount of difference should one expect between the groups (given the data)?". These questions natually deal with means. Going to the log-scale asks for relative differences instead of absolute differences.

This is what I call "making sense of data", and this is what statistics shoul be for (rather than producing p-values and "significance stars" at any cost).

Best wishes!

Great answer, Jochen. 

Iona, you may also find the nice simulation study by Fagerland & Sandvik useful (link below).  It is concerned with the Wilcoxon-Mann-Whitney test for two independent groups (not repeated measures), but the general point is still relevant.

HTH.

Article The Wilcoxon-Mann-Whitney test under scrutiny

For your HFn values, you could also try a log transform. The HR values do look normally distributed.

Here is some R code for a repeated measures ANOVA of the HFn scores. If you want to exactly the same results as in SPSS you might have to use the Anova command from the car package (implementing Type III sums of squares). Either way, I am pretty sure that you will get a main effect of Time, but no interaction with group. I leave it as an exercise to adapt the code to the HR scores (which should not be log transformed).

library(reshape2)

library(ggplot2)

xx

Thank you Bruce and Jochen (and everyone)! I am trying all of these. As I said, I did get an interaction effect for HF in SPSS using ANCOVA (covariating for the first baseline), but I was skeptical about that.

Best

Ioana

Try with this package rptR (repeatability estimation for gaussian and non-gaussian data)