Calculation of the simple effects is not working. This seems to be an attempt to calculate relative variance, as per partial eta, but the numbers are not matching up in the model I have. My model, as pasted below (this is moderation1) cannot calculate the etasq because (5,3) does not exist.
Value df SS F p
Unequal 0.2484115 1 1.216679 0.3168459 0.5743497
Equal -0.4111375 1 3.328362 0.8667674 0.3533478
Residuals . NA 150 575.995752 NA NA

create interaction model of interest

modinter <- na.omit(model1)

examine simple effects

(moderation1 <- testInteractions(modinter, fixed = c("iv2", "iv3"), across="iv1", adjustment="none"))

etasqiv1.l_l <- (moderation1[1,3]) / ((moderation1[1,3]) + (moderation1[5,3])) # eta-squared iv2-l, iv3-l, across iv3
etasqiv1.h_l <- (moderation1[2,3]) / ((moderation1[2,3]) + (moderation1[5,3])) # eta-squared iv2-h, iv3-l, across iv3
etasqiv1.l_h <- (moderation1[3,3]) / ((moderation1[3,3]) + (moderation1[5,3])) # eta-squared iv2-l, iv3-h, across iv3
etasqiv1.h_h <- (moderation1[4,3]) / ((moderation1[4,3]) + (moderation1[5,3])) # eta-squared iv2-h, iv3-h, across iv3

In the section with the correlation matrix instead of having iv1, you have iv1_num. I did not convert my IVs to numerics and they still ran. If you want people to convert their variables so the function works, please explicitly add instructions for how people can convert their variables.

dat3 = na.omit(dat %>%
dplyr::select(iv1_num, iv2, dv) %>%
rename(NEW_NAME_IV1 = iv1_num, # relabel whatever you want your variables to be named in the manuscript, cannot contain spaces though
NEW_NAME_IV2 = iv2,
NEW_NAME_DV = dv))

The following code works, as long as the DV actually had the parameters 1-9. As I was testing it with data between -1 and 1, the graph didn't load any data until I changed the y-limits/y-breakdown. Consider having the user input these values into an object and putting the object in the code.

(figure1 = ggplot(dat1, aes(x = iv1, y = mean, group = iv2)) +
geom_bar(stat = "identity", position = "dodge", aes(fill = iv2)) +
geom_errorbar(limits, position=dodge, width=0.25) +
apatheme +
coord_cartesian(ylim=c(-1,1)) + # this is the range of the y-axis
scale_y_continuous(breaks=seq(-1,1,.5)) + # this is the increment ticks on the y-axis
ylab('dv') +
xlab('iv1') +
scale_fill_grey())

Same issue as #7, script is missing objects required to run the figure.

This is supposed to test a 3-way model, but on line 64 the function only requires two IVs instead of three. See below taken from the 3way script

aov_models <- dat %>%
dplyr::select(starts_with("avg_")) %>% # this line tells the map() only use your DVs (all start "avg_" in my datasets)
map(~Anova(lm(. ~ iv1 * iv2, data = dat), type = 3))

aov_models

It feels like I'm looking at part of a script. The figure code requires objects such as y_axis_low, y_axis_high, and y_increment, but I do not see these objects being defined in the script.

Also, why is calling figure_1 and the code for saving the figure commented out?

On line 75 you have a mis-type, only need 1 DV:

linear model to dive into any significant 2-ways

model1 <- lm(dv ~ dv ~ iv1 * iv2 * iv3, data=dat)