Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Different results from JASP when I compare to ez package (repeated measures) #609

Closed
anovabr opened this issue Feb 10, 2020 · 9 comments
Closed
Assignees

Comments

@anovabr
Copy link

anovabr commented Feb 10, 2020

  • JASP version: 0.11.1
  • OS name and version: Windows 10
  • Analysis: Repeated measure (and ANCOVA)
  • Bug description: Different results from R (base--ez package) and JASP
  • Expected behavior: I imagine both results should be similar

Hello, I'm running two different analyses in JASP and R. The first analysis is an ANCOVA and the second one is a repeated-measures ANOVA. As always, people want to match the results across software and that drives my question. In the place where I work, R (ez package) is the default software and people will remain with its results.

When I compare the first analysis (ANCOVA), the results look great. They are similar between R (base) and JASP.

image

This makes me feel comfortable to ask the target-question.

I changed the format from my dataset to long and then I ran the RM ANOVA. Unfortunately, the output was not parallel. Just bolding some results, the main effects are different.
From ez package,
latino_group is 252
From JASP,
209
Condition from ez package is 144, but it is 172 from JASP.

Please check the image below:

image

Any comments are valuable.
Thank you,

The entire code is this one:
`#ancova
aov_outcome <- lm(M_4I2TOT ~ M_4I1TOT + factor(condition) * factor(latino_group), data = temp)
car::Anova(aov_outcome, type=3)
#results are ok

repeated measures (long format)

temp %>%
select(ID, latino_group,condition, M_4I1TOT,M_4I2TOT) %>%
pivot_longer(cols = M_4I1TOT:M_4I2TOT, #variables that have results
names_to = "time",
values_to = "results") -> temp_long

library(ez)
options(contrasts = c("contr.sum","contr.poly"))
ez_outcome <- ezANOVA(
data = temp_long,
dv = results,
wid = ID,
within = time,
between = .(latino_group, condition),
type = 3,
detailed = TRUE,
return_aov = TRUE)
summary(ez_outcome$aov)`

The dput of the dataset is below:
dput(temp) structure(list(ID = structure(c(10002, 10006, 10009, 10010, 10017, 10018, 10023, 10025, 10026, 10031, 10035, 10036, 10042, 10043, 10047, 10048, 10049, 10061, 10065, 10072, 10077, 10081, 10082, 10083, 10085, 10086, 10087, 10089, 10090, 10094, 10100, 10104, 10105, 10106, 10110, 10111, 10112, 10114, 10116, 10117, 10121, 10122, 10144, 10147, 10148, 10150, 10153, 10155, 10164, 10198, 10201, 10204, 10206, 10223, 10224, 10226, 10234, 10239, 10241, 10242, 10243, 10248, 20002, 20003, 20004, 20012, 20021, 20026, 20028, 20047, 20048, 20049, 20081, 20093, 20106, 20107, 20114, 20127, 20129, 20133, 20160, 20165, 20169, 20197, 20202, 20205, 20206, 20216, 20217, 20219, 20220, 20221, 20224, 20225, 20229, 20230, 20231, 20236, 20243, 20284, 20285, 20292, 20296), label = "ID", format.spss = "F5.0"), latino_group = c("white", "white", "white", "white", "white", "white", "latino", "latino", "white", "latino", "white", "latino", "white", "white", "latino", "white", "latino", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "latino", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "latino", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "white", "latino", "white", "latino", "latino", "latino", "white", "white", "latino", "latino", "latino", "latino", "white", "latino", "latino", "latino", "latino", "latino", "latino", "latino", "latino", "latino", "latino", "latino", "latino", "white", "white", "latino", "latino", "white", "latino", "white"), condition = c("attention-control", "attention-control", "experimental", "attention-control", "experimental", "experimental", "attention-control", "attention-control", "attention-control", "attention-control", "experimental", "experimental", "experimental", "experimental", "attention-control", "attention-control", "experimental", "experimental", "attention-control", "experimental", "attention-control", "attention-control", "attention-control", "experimental", "attention-control", "attention-control", "experimental", "experimental", "experimental", "experimental", "attention-control", "experimental", "attention-control", "attention-control", "experimental", "attention-control", "experimental", "experimental", "experimental", "attention-control", "attention-control", "attention-control", "experimental", "experimental", "attention-control", "attention-control", "experimental", "attention-control", "experimental", "attention-control", "experimental", "attention-control", "experimental", "experimental", "attention-control", "attention-control", "experimental", "attention-control", "experimental", "experimental", "attention-control", "experimental", "experimental", "attention-control", "attention-control", "attention-control", "experimental", "experimental", "experimental", "experimental", "experimental", "attention-control", "experimental", "experimental", "attention-control", "experimental", "experimental", "attention-control", "experimental", "attention-control", "experimental", "attention-control", "attention-control", "experimental", "attention-control", "attention-control", "attention-control", "experimental", "attention-control", "experimental", "experimental", "attention-control", "experimental", "attention-control", "experimental", "attention-control", "experimental", "experimental", "experimental", "attention-control", "attention-control", "experimental", "experimental"), M_4I1TOT = c(7, 3, 8, 12, 7, 9, 3, 4, 8, 3, 5, 1, 12, 8, 2, 6, 10, 7, 7, 8, 9, 5, 12, 7, 9, 8, 8, 7, 13, 10, 7, 5, 4, 11, 10, 4, 4, 7, 7, 7, 3, 8, 10, 3, 7, 9, 8, 11, 7, 3, 4, 12, 8, 8, 7, 13, 10, 7, 13, 7, 8, 5, 10, 12, 5, 5, 2, 5, 9, 3, 6, 8, 7, 7, 1, 4, 7, 5, 8, 4, 3, 4, 4, 5, 6, 7, 2, 7, 6, 9, 2, 5, 5, 3, 6, 1, 10, 5, 10, 7, 3, 3, 5), M_4I2TOT = c(9, 2, 13, 8, 7, 12, 3, 1, 9, 6, 9, 7, 8, 8.2078195323884, 3, 5, 12, 10, 6, 6, 9, 7, 12, 8, 6, 5, 9, 6, 13, 13, 11, 13, 5, 9.17142301161172, 12, 8, 4, 13, 13, 6, 11, 9, 13, 3, 6, 10, 13, 8, 11, 5, 4, 10.2470110225081, 12, 10, 9, 10, 13, 13, 13, 12, 5, 9, 11, 11, 8, 4, 6.76864538709267, 13, 13, 11, 11, 5, 11, 8.06448852972179, 5, 9, 7, 6, 12, 4, 10, 5, 5, 2, 4, 8, 2, 10, 7, 11, 4, 8, 10, 2, 9, 2, 8, 7, 10, 7, 5, 12, 10)), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, -103L))

The attached file is the CSV file of my dataset.

temp.zip

@JorisGoosen
Copy link
Contributor

Hey @JohnnyDoorn could you have a look at this?
If it is a bug somehow it would be good to get it out before 0.12.

@JohnnyDoorn
Copy link

Hi @anovabr,

Thanks for taking the effort to discuss this. I just looked at your code in R, and recreated JASP's results to see where the difference lies. JASP uses the afex package, and you can reprodce JASP's results with the following R code:
library(afex)
result <- afex::aov_car(results ~ conditiontimelatino_group + Error(ID/time), data=datLong, type= 3, include_aov = TRUE)
summary(result)

I tried to reproduce EZ's output with afex, and it seems that if I specify the sum of squares type to 2, it is almost the same (except for the latino_group variable). I think that due to the unbalanced design, the sums of square type matters a lot. It is weird that ez produces these different results, as you specify the type to 3. Upon closer inspection, the ezANOVA function returns exactly the same output for specifying type to 1, 2, or 3, which seems dubious.
I am not sure what is going on there, but I have more faith in afex, as this package actually gives different results based on the different sums of squares (which it should). What do you think?

Kind regards,
Johnny

@anovabr
Copy link
Author

anovabr commented Feb 11, 2020

I don't know either what's going on with ez package, I'll try to reach out to the developers and discover this issue. I always used ez package and now I'm afraid of its previous results...
Thanks for your support. I hope I had contributed to the improvement of the JASP. Keep up the good work you're doing with this software.

@mike-lawrence
Copy link

If I recall correctly, afex simply calls ez beneath the hood, so I'm surprised by this discrepancy. I'll dive in now.

@mike-lawrence
Copy link

mike-lawrence commented Feb 11, 2020

So the main results returned by ezANOVA match afex/jasp:

#load packages 
library(tidyverse)
library(afex)
library(ez)

#read data & reshape to long
readr::read_csv('temp.csv') %>% 
	select(ID, latino_group,condition,  M_4I1TOT,M_4I2TOT) %>%
	pivot_longer(cols = M_4I1TOT:M_4I2TOT, #variables that have results
				 names_to = "time",
				 values_to = "results") -> temp_long

#make sure variables are factors and have their contrasts (if any) are explicitly 
# set (so we're not relying on defaults)
temp_long$ID = factor(temp_long$ID)

temp_long$time = factor(temp_long$time)
contrasts(temp_long$time) = 'contr.sum'

temp_long$condition = factor(temp_long$condition)
contrasts(temp_long$condition) = 'contr.sum'

temp_long$latino_group = factor(temp_long$latino_group)
contrasts(temp_long$latino_group) = 'contr.sum'

#run afex
afex_outcome <- aov_car(results ~ condition * time * latino_group + Error(ID/time),
						data=temp_long,
						type= 3,
						include_aov = TRUE)
summary(afex_outcome)

#run ez
ez_outcome <- ezANOVA(
	data = temp_long,
	dv = results,
	wid = ID,
	within = time,
	between = .(condition,latino_group),
	type = 3,
	detailed = TRUE,
	return_aov = TRUE
	)

ez_outcome$ANOVA

#check that all p-values are precisely equal
all(ez_outcome$ANOVA$p==summary(afex_outcome)$univariate.tests[ez_outcome$ANOVA$Effect,6])

But as you've observed, the $aov element provided by ezANOVA when return_aov is TRUE isn't matching up (even with ez's own summary). Digging into that now.

@mike-lawrence
Copy link

mike-lawrence commented Feb 11, 2020

Ah, I think I see what's going on now. The return_aov argument for ez was added in 2011 as a niche feature requested by a solitary user that wanted the equivalent stats::aov object for some post-hoc contrast calculations. It was simple to add but I didn't even think to address the unbalanced-designs case. Normally this wouldn't be a huge deal as all the ez documentation shows that the $ANOVA result from ezANOVA is the thing you want to look at for results, but I'm guessing that after afex came out people have started getting confused between it's "include_aov" and ez's "return_aov", expecting they are the same thing and/or designed for the same purpose.

I'll add a warning about return_aov being broken for unbalanced designs (I'm unlikely to fix since afex has taken over popularity & feature-wise; I'm also now philosophically against the whole concept of making a blunt/no-thinking tool like ANOVA even easier to access, but that's a whole separate discussion).

@JohnnyDoorn
Copy link

Hi @mike-lawrence,

Thank you for your quick response and diving into the issue, I am glad it is sorted out now! I think ez is a great tool (just showed it to my R students 2 weeks ago), but it's indeed an interesting point about the no-thinkingness, and something we also run into with JASP...

I am closing this issue now, @anovabr please reopen if you still have any questions!

Kind regards,
Johnny

@singmann
Copy link

Ah, I think I see what's going on now. The return_aov argument for ez was added in 2011 as a niche feature requested by a solitary user that wanted the equivalent stats::aov object for some post-hoc contrast calculations.

I believe this user was me, on stats.stackexchange. See my comment here:
https://stats.stackexchange.com/a/6212/442

Let me also address one more thing:

If I recall correctly, afex simply calls ez beneath the hood, so I'm surprised by this discrepancy. I'll dive in now.

If this is not clear from the discussion above, that is not correct. afex does not call ez at all. Both call car::Anova for calculating ANOVAs. The only code that afex uses from ez is for calculation of generalized eta-squared (this code is copied with appropriate attribution). All other afex code is written independent of ez. [In the sense of without looking at the ez code. The design of afex is of course highly influenced by ez, hence the function name aov_ez.]

@mike-lawrence
Copy link

mike-lawrence commented Feb 12, 2020

Let me also address one more thing:

If I recall correctly, afex simply calls ez beneath the hood, so I'm surprised by this discrepancy. I'll dive in now.

If this is not clear from the discussion above, that is not correct.

Yes, my apologies for this self-aggrandizing mis-memory. I knew we were connected somehow (we both call car::Anova as you say), but should have investigated my vague recollection before making any comment.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

5 participants