I am passionate about teaching statistics and programming skills for use in psychology. My weapons of choice:
R for data management and statistical analyses
Python, PsychoPy, and Pavlovia for building and deploying in-lab/online experiments
R for data management and statistical analyses
Python, PsychoPy, and Pavlovia for building and deploying in-lab/online experiments
|
|
See my OSF page for R scripts from my projects.
Custom R Toolbox
I have written a few custom R functions for personal use. It's not the most efficient code, nor is it dependency-free. Feel free to use if useful!
Eta-squared partial for aov_ez models from library(afex)
# eta squared partial function
# this is because library(afex) only reports generalized eta square, which is technically the better effect size measure for repeated-measures ANOVA designs
# alas, it is a necessary evil to also know how to compute effect sizes other researchers typically use and report
# requires library(apaTables)
# x = model object name from library(afex)
# n = which # line from anova table you want
# ci = how wide is the confidence interval (i.e., 0.90, 0.95)
etasqp <- function(x, n, ci) {
require(apaTables)
sseffect <- summary(x)$univariate.tests[n,1]
sserror <- summary(x)$univariate.tests[n,3]
es <- sseffect/(sseffect+sserror)
f <-summary(x)$univariate.tests[n,5]
ndf <- summary(x)$univariate.tests[n,2]
ddf <- summary(x)$univariate.tests[n,4]
cis <- get.ci.partial.eta.squared(F.value=f, df1=ndf, df2=ddf, conf.level=ci)
stuff <- c("partial eta sq",es,"confidence intervals",cis)
return(stuff)
}
Emmeans but with up to 2 decimal places
# emmeans() rounds marginal means to one decimal place, and I wrote this function so that it would Not Do That
# requires library(emmeans)
# emm_object = emmeans(lm_object, specs = pairwise ~factor1:factor2)
better_emmeans <- function(emm_object) {
x <- data.frame(group = summary(emm_object)$emmeans[[1]],
marg_mean = round(summary(emm_object)$emmeans[[2]],2),
se = round(summary(emm_object)$emmeans[[3]],2))
return(x)
}
Post-hoc t-test with Cohen's d and 95% CI for emmeans object
# allows you to test any 2 levels of a factor in a model (lm/aov), from an emmeans object
# and it spits out a handy dandy cohen's d and associated 95% CIs
# requires library(MBESS)
# emm_object = emmeans(lm_object, specs = pairwise ~factor1:factor2)
# (flexible for main effect or interaction, just remove `:factor2` above)
# group 1 = what row # is the first group in emm_object$emmeans
# group 2 = what row # is the other group in emm_object$emmeans
# contrast_row = what row # is the $contrast of interest
phtt <- function (emm_object, group1, group2, contrast_row) {
require(MBESS)
tstat <- summary(emm_object)$contrasts[[5]][contrast_row]
n1 <- emm_object[["emmeans"]]@grid[[".wgt."]][group1]
n2 <- emm_object[["emmeans"]]@grid[[".wgt."]][group2]
results <- ci.smd(ncp = tstat,
n.1 = n1,
n.2 = n2,
conf.level=0.95)
x <- data.frame(contrast = summary(emm_object)$contrasts[[1]][contrast_row],
t = round(tstat, 2),
df = summary(emm_object)$contrasts[[4]][contrast_row],
p = round(summary(emm_object)$contrasts[[6]][contrast_row], 4),
cohend = round(results$smd, 2),
lowerCI = round(results$Lower.Conf.Limit.smd, 2),
higherCI = round(results$Upper.Conf.Limit.smd, 2))
return(x)
# reproducible example. requires library(emmeans), and my custom function phtt() in your environment
mod1 <- lm(weight ~ feed, data = chickwts) # you should default have chickwts in R environment
emm_object = emmeans(mod1, specs = pairwise ~ feed)
emm_object # check which contrast you want
# I'm choosing soybean ($emmeans row 5) vs. sunflower ($ emmeans row 6), which is row 15 in $contrasts
phtt(emm_object, group1 = 5, group = 6, contrast_row = 15)
# you should get these results:
# contrast: soybean - sunflower
# t = -3.82
# df = 65
# p = 0.0039
# cohen d = -1.5
# lowerCI = -2.37
# higherCI = 0.61
Photo on banner was taken at Whistler, in British Columbia, Canada.
