I am supporting a psychology experiment, and having problems analyzing some of the data.

By way of background, I’m a programmer at a research organization who’s taken several stat courses recently. A relatively junior researcher (Ph.D. candidate) asked me if I would help her analyze some data in preparation for an internal review.

The experiment involves measuring pupil diameter as a function of cognitive load.

In a nutshell, a person sits in front of a computer screen for about 20 minute reading. From time to time the subject is given a task, such as performing a mathematical calculation, that is meant to induce cognitive load. Sensors at the bottom of the computer monitor measure pupil dilation on the order of every few milliseconds. The idea is that under cognitive load, pupils have been shown to dilate.

The overall 20-minute time periods are divided into four sections, each of which has a different numbers of tasks. A representative 20-minute session looks like this:

Section         A       B       C       D
# Observations  13908   20416   19635   23269
# Tasks         0       11      28      6

At the end of a session, the data might look something like this:

Section         A           B       C       D
# Observations  13908       20416   19635   23269
# Tasks         0           11      28      6
Mean Dilation   4.073       4.195   4.216   4.118
StdDev          0.202       0.254   0.256   0.239

So the task is to determine if there are significant differences between the mean dilations across the sections, and I was asked to run one-way ANOVAs on the data.

In R I did the following on each data set:

df$section <- as.factor(df$section)

testANOVA <- oneway.test(dilation ~ section, data = df)


I got outputs like the following:

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

data:  df$dilation and df$section
F = 1436.55, num df = 3.00, denom df = 41130.44, p-value < 2.2e-16

I believe 2.2e-16 is the lowest p-value that R can give, and it did not seem right to me. I basically interpreted this as, “There is almost no possible way that the null hypothesis could be true.”

I'm relatively new to statistics and haven’t studied ANOVAs yet, but in just looking at the means and standard deviations of all the data sets it does not seem reasonable that the null hypothesis would be so strongly rejected across almost the entirety of the data sets.

The researcher told me early on not to worry about assumptions, but given those results, I started digging into the assumptions of a one-way ANOVA.

These are some of the issues I came up with regarding this data:

--Unequal variance. The variances are not equal across the blocks (although my understanding of the oneway.test in R is that it accounts for unequal variance).

--Lack of independence. The data from each segment comes from the same person during the same session.

--Lack of normalcy. Each block exhibits normal histograms when looked at in its entirety, but I started taking random samples of each block (without replacement), and once I got down to below 1,000 observations in a block, it became apparent that some of the blocks were pretty skewed, and in some cases the histogram almost looked like a uniform distribution.

Also of note, the segment sizes are different (greatly so in some cases).

I tried two different approaches to getting some insights: (1) I created simulated data using the real means and standard deviations, but imposing normal distributions; (2) I standardized the size of the segments of the real data by sampling (without replacement).

In both of the above scenarios, I noticed that with very large numbers of observations, I still invariably got very low p-values. If I started reducing the number of observations, the p-values started rising.

In digging deeper, I came across this post in which someone says that with a large number of observations, normality tests are essentially useless:

Is normality testing essentially useless?

Which made me wonder if the same holds true for one-way ANOVAs.

I’m obviously way out of my league at this point, and I’m going to tell the researcher that. But for my own edification I have a simple question:

Given the way I’ve described the data, (1) are there ways – via transformations and perhaps non-parametric approaches – to perform ANOVAs to glean some sorts of accurate insights into segment-level differences, or (2) with all of the assumption issues and whatnot, might there be a fundamental flaw in the design of the experiment such that it’s generally producing just noise?

If it's the first, I'd be inclined to keep digging around, just to educate myself. But if it's the latter, I don't really want to spend a whole lot of time on it.

  • $\begingroup$ Are you running separate tests on each person, or are you running one test on everyone? $\endgroup$ – gung - Reinstate Monica Feb 10 '15 at 14:55
  • $\begingroup$ @gung. There are two different configurations of the test, and within each configuration the order of the segments is sometimes shuffled. $\endgroup$ – rwjones Feb 10 '15 at 15:54
  • $\begingroup$ No, I mean are you running separate ANOVAs on each person, or are you running one ANOVA on everyone? The latter would be normal, but from your description I wonder if you are doing the former. $\endgroup$ – gung - Reinstate Monica Feb 10 '15 at 16:17
  • 1
    $\begingroup$ Those are some large-ish sample sizes. Note that frequentist hypothesis tests are biased towards rejecting the null as sample size gets big. $\endgroup$ – Alexis Feb 10 '15 at 19:20
  • 1
    $\begingroup$ Agree with @Glen_b that temporal autocorrelation is likely to be very important. I would start by binning the data by some reasonable temporal frequency (1 second? 5 seconds?), then compute and plot acf() for each segment. $\endgroup$ – Ben Bolker Feb 11 '15 at 2:27

The assumption of normality is on the residuals, not on the original data. No worries there. If you are worried about homogeneity of variance try a transformation. I'm a little confused if these observations are all from one person, or from many people? You may have a nested design.

You have a very large n (observations), so it's not entirely surprising that you have a significant p-value. It is the job of the researcher to ask the next question - is statistical significance in this case physiologically significant?

In any case, I'd recommend a linear regression, rather than ANOVA, for example, lm(dilation ~ tasks, mydata). With an ANOVA you lose some "predictive" power. If you do have a nested design, I'd look into the lme4 package and mixed linear models using lmer().

Cheers and good luck!

  • $\begingroup$ As I mentioned to gung, the observations are from 20 different people. There are two different configurations of the test, and within each configuration the segments might get shuffled. $\endgroup$ – rwjones Feb 10 '15 at 15:56
  • $\begingroup$ In that case, these are nested random effects. I would consider lmer() from lme4. Something similar to: lmer(dilation ~ tasks + (1|configuration:tasks) + (1|people:configuration), mydata). $\endgroup$ – danno Feb 10 '15 at 16:19
  • $\begingroup$ I was wondering about your comment that the assumption of normality applies to residuals, not data. On this site, for example, they say it's about the underlying population GraphPad Statistics. $\endgroup$ – rwjones Feb 10 '15 at 16:49
  • $\begingroup$ It's a common misconception. See <theanalysisfactor.com/the-distribution-of-dependent-variables>. You should use a Q-Q plot to investigate normality of residuals. Do not bother doing a test of normality on the residuals - instead make a judgement from the Q-Q plot. ANOVAs have some wiggle room here. If you are very concerned, use a non-parametric ANOVA. However, I still suggest you look into mixed linear models, as you have nested random effects that will be difficult to incorporate into other statistical tests. $\endgroup$ – danno Feb 10 '15 at 17:14
  • $\begingroup$ I have done a lot of reading since first posting this question. First, yes, I do believe that I need to look at the lme4 package and mixed linear models. Second, since each subject in the experiment has four treatments sequentially, I don't believe it's a nested design (at least as I understand it). Third, I continue to be confused about your normality comments. The link you provided does not work, and literally every article or book I've seen says that normality of the independent variables (not the residuals) is the thing. Anyway, overall your guidance was spot on. Thanks. $\endgroup$ – rwjones Mar 4 '15 at 20:11

The variances are so nearly equal that any impact of unequal variances on inference will be small.

Normality in such large samples may not be a big issue either. With such big samples you should have good level-robustness, and if you have so much data that the potential adverse impact on power may be a nonissue. If this was the only problem, you could always look at permutation tests for this.

The main issue as I see it will be that you have sequences of measurements of dilation over time for a given subject -- effectively repeated measures/time series data.

These will be dependent (two measurements taken a few milliseconds apart are much more likely to be very similar than two measurements taken minutes apart). You can't avoid this issue and must tackle it head on, I think by explicitly including the tendency to be correlated in any modelling and any inference.

Does the mean of dilation change over time (what does a time series plot look like, for example)? What does the ACF and PACF of each time series look like (or of the residuals for some suitable model for the mean)?

This can certainly impact significance levels, perhaps a great deal.

  • 1
    $\begingroup$ @danno's answer is reasonable, but I think autocorrelation is the crux. +1. $\endgroup$ – Ben Bolker Feb 11 '15 at 2:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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