R O N D A ‎‏‏‎ ‏‏‎ ‎‎‎‎F. ‏‏‎‏‏‎ ‎ ‎L O
  • Home
  • About Me
  • CV
  • Publications
  • R and Statistics
  • Contact info

R AND STATISTICS

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
Picture
Picture
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.
Powered by Create your own unique website with customizable templates.
  • Home
  • About Me
  • CV
  • Publications
  • R and Statistics
  • Contact info