Calculate omega-squared in R

Calculate omega-squared in R.R

# Compute omega-squared and partial omega-squared
# By Arnoud Plantinga
# Based on http://stats.stackexchange.com/a/126520

# Functions ---------------------------------------------------------------

# Omega-squared
Omegas <- function(mod){
  aovMod <- mod
  if(!any(class(aovMod) %in% 'aov')) aovMod <- aov(mod)
  sumAov     <- summary(aovMod)[[1]]
  residRow   <- nrow(sumAov)
  msError    <- sumAov[residRow,3]
  dfEffects  <- sumAov[1:{residRow-1},1]
  ssEffects  <- sumAov[1:{residRow-1},2]
  msEffects  <- sumAov[1:{residRow-1},3]
  ssTotal    <- rep(sum(sumAov[1:residRow, 2]), 3)
  Omegas <- abs((ssEffects - dfEffects*msError)/(ssTotal + msError))
  names(Omegas) <- rownames(sumAov)[1:{residRow-1}]
  Omegas
}

# Partial omega-squared
partialOmegas <- function(mod){
  aovMod <- mod
  if(!any(class(aovMod) %in% 'aov')) aovMod <- aov(mod)
  sumAov     <- summary(aovMod)[[1]]
  residRow   <- nrow(sumAov)
  dfError    <- sumAov[residRow,1]
  msError    <- sumAov[residRow,3]
  nTotal     <- nrow(model.frame(aovMod))
  dfEffects  <- sumAov[1:{residRow-1},1]
  ssEffects  <- sumAov[1:{residRow-1},2]
  msEffects  <- sumAov[1:{residRow-1},3]
  partOmegas <- abs((dfEffects*(msEffects-msError)) /
                      (ssEffects + (nTotal -dfEffects)*msError))
  names(partOmegas) <- rownames(sumAov)[1:{residRow-1}]
  partOmegas
}

# Example -----------------------------------------------------------------

data(ToothGrowth)

lm1 <- lm(len ~ supp * dose, data=ToothGrowth)
summary(lm1)

# Eta-squared
require(lsr)
etaSquared(lm1)

# Omega-squared
Omegas(lm1)
partialOmegas(lm1)

This is amazing! Thanks for posting this.

If you want vertical output in the Omegas function, replace:

  names(Omegas) <- rownames(sumAov)[1:{residRow-1}]
  Omegas

with:

  result <- data.frame(Omegas)
  rownames(result) <- rownames(sumAov)[1:{residRow-1}]
  colnames(result) <- "Omega-squared"
  result

It does not affect the results, but you should change 3 to residRow-1 for ssTotal definition (https://gist.github.com/arnoud999/e677516ed45e9a11817e#file-calculate-omega-squared-in-r-r-L17).

On the other hand, you do not need to replicate the vector, because you divide vector by the scalar in the next line.

This should work with type I type II and type III sum of squares. I am not sure if everything is correct in there but I think it is an improvement.

compute_omega<-function(model,ss="I") {
  n_total=nrow(model.frame(model))
  ss1<-data.frame(summary(model)[[1]],check.names=FALSE)
  ss2<-data.frame(car::Anova(model,type="II"),check.names=FALSE)
  ss3<-data.frame(car::Anova(model,type="III"),check.names=FALSE)
  ss2$`Mean Sq`<-ss2$`Sum Sq`/ss2$Df
  ss3$`Mean Sq`<-ss3$`Sum Sq`/ss3$Df
  ss2<-ss2[,c("Df","Sum Sq","Mean Sq","F value","Pr(>F)")]
  ss3<-ss3[,c("Df","Sum Sq","Mean Sq","F value","Pr(>F)")]
  
  if(ss=="I")
    summary_aov<-ss1
  if(ss=="II")
    summary_aov<-ss2
  if(ss=="III"){
    summary_aov<-ss3[2:nrow(ss3),]
    intercept<-ss3[1,]
  }
  
  residual_row<-nrow(summary_aov)
  ms_error<-summary_aov[residual_row,3]
  df_effects<-summary_aov[1:(residual_row-1),1]
  ss_effects<-summary_aov[1:(residual_row-1),2]
  ms_effects<-summary_aov[1:(residual_row-1),3]
  ss_total<-rep(sum(summary_aov[1:residual_row,2]),3)
  omega<-abs((ss_effects-df_effects*ms_error)/(ss_total+ms_error))
  names(omega)<-trimws(rownames(summary_aov)[1:(residual_row-1)])
  df_error<-summary_aov[residual_row,1]
  partial_omega<-abs((df_effects*(ms_effects-ms_error))/(ss_effects+(n_total-df_effects)*ms_error))
  names(partial_omega)<-trimws(rownames(summary_aov)[1:(residual_row-1)])
  
  if(ss=="III")
    summary_aov<-rbind_all(intercept,summary_aov)
  
  summary_aov<-data.frame(ss=ss,
                          comparisons=trimws(row.names(summary_aov)),
                          summary_aov,
                          check.names=FALSE)
  result<-data.frame(ss=ss,
                     comparisons=names(omega),
                     omega=omega,
                     partial_omega=partial_omega)
  
  result<-merge(summary_aov,result,all=TRUE,sort=FALSE)
  return(result)
}