# Rubin's rule, applied to absolute effect size or relative effect size (Cohen's d)?

Cohen's d is a way to describe the effect size relative to the standard deviation of the data.

For instance in the case of the difference between the means of two populations

$$\begin{array}{} \text{absolute effect size} &=& \bar{x_1} - \bar{x_2} \\ \text{relative effect size} &=& \frac{\bar{x_1} - \bar{x_2}}{\hat\sigma} &=& \text{Cohen's d} \\ \end{array}$$

If we want to apply Rubin's Rules to pool the results of multiple types of imputation of the same data, then should we apply the rules to the absolute effect size or to cohen's d?

Say, we have the following two approaches. We have some data and different imputations of it. We could do one of the following two:

• Compute the different $$d$$ values and their standard error for the different imputations and apply the Rubin's Rules to it to get a pooled $$d$$ and it's variance.
• Compute the different absolute effects $$\bar{x_1} - \bar{x_2}$$ and the different population variance estimates $$\hat{\sigma}$$ for the different imputations and apply the Rubin's Rules to each seperately. From those two results compute a pooled $$d$$ and it's variance.

Can the second approach be done (or maybe some other way of applying the rules to the absolute effects) and could it be more accurate?

Below are results for a comparison made by 1000 simulations (R-code below)

The steps taken are:

• generate two vectors We use two samples of sizes 20 and normal distributions $$X \sim N(0,1)$$ and $$Y \sim N(0.2,1)$$ For this case Cohen's d would be equal to 0.2

• randomly delete data We selected randomly 10 values out of the 40 values.

• impute the data For this I used the mice r-package. The standard method that it uses is predictive mean matching.

• compute Cohen's d This is done in multiple ways.

1. In one case I apply Rubin's rule to the $$d$$ computed for each imputation.
2. In another case, I apply Rubin's rule to the difference of the means and to the estimate of the standard deviation. Finally, to compute $$d$$ I take the ratio of those two.
3. A third case. In this simplified case imputation is not at all needed. We can compute Cohen's d for two samples of unequal size and there is no need to fill the missing data.

Imputing the data is usefull when the method can not deal with missing data (for instance a linear regression where some of the regressor values $$X$$ are missing). Also possible might be that one wishes to create more balanced data and use imputation to generate equal representation of classes in both variables (e.g. demographic data), but in this case there might be potentially alternatives by applying different weights based on the classes or by using a linear mixed effects model.

The results are

Method           Mean of 1000 simulations  Mean of squared error
1 relative       0.2124757                 0.1502107
2 absolute       0.2118989                 0.1493271
3 no imputation  0.2033886                 0.1302398


Code:

library(mice)

nd <- 20
nimp <- 10
effect <- 0.2
m <- 5

### function to create data with imputations from mice
set <- function(imp, n) {
data_imp <- imp$$data data_imp[imp$$where[,1],1] <- imp$$imp[[1]][,n] data_imp[imp$$where[,2],2] <- imp\$imp[[2]][,n]
data_imp
}

### function to do the simulation
sim <- function() {
### generate data
data <- data.frame(x = rnorm(nd,0,1), y = rnorm(nd,effect,1))

### sample nimp points to remove (make equal to NA)
imp <- sample(0:(nd*2-1),nimp)
for (i in imp) {
k <- floor(i/30)+1
l <- (i%%30)+1
data[l,k] <- NA
}

### perform 5 imputations with mice
imp2 <- mice(data[], m = m, printFlag = F)

### vectors to contain results from for-loop
ve_d <- rep(0,5)
sd_d <- rep(0,5)

ve_m <- rep(0,5)
sd_m <- rep(0,5)

ve_v <- rep(0,5)
sd_v <- rep(0,5)

for (i in 1:m) {  ### repeatedly compute statistics for different imputations
m1 <- set(imp2,i)
par1 <- mean(m1[,2])-mean(m1[,1])            ### difference of means
par2 <- sqrt(0.5*(var(m1[,1])+var(m1[,2])))  ### pooled variance estimate
par3 <- par1/par2                            ### Cohen's d
nu <- 2*nd-2

### standard deviations of par1, par2, par3
sd1 <- par2*sqrt(2/nd)
sd2 <- par2 * sqrt(gamma(nu/2)/gamma((nu+1)/2) * nu/2- 1)
sd3 <- sqrt((nd+nd)/(nd*nd) * (nu)/(nu-2)+ par3^2 * (nu/(nu-2)-(1-3/(4*nu-1))^-2))

ve_m[i] <- par1
sd_m[i] <- sd1
ve_v[i] <- par2
sd_v[i] <- sd2
ve_d[i] <- par3
sd_d[i] <- sd3
}

### Rubin's Rule applied to Cohen's d
est_d <- mean(ve_d)
var_d <- mean(sd_d^2) + var(ve_d) * (1+1/m)

### Rubin's Rule applied to difference
est_m <- mean(ve_m)
var_m <- mean(sd_m^2) + var(ve_m) * (1+1/m)

### Rubin's Rule applied to variance
est_v <- mean(ve_v)
var_v <- mean(sd_v^2) + var(ve_v) * (1+1/m)

### Straightforward Cohen's d
x <- data$$x[!is.na(data$$x)]
y <- data$$y[!is.na(data$$y)]
nx <- length(x)
ny <- length(y)
cohend <- (mean(y)-mean(x))/sqrt( ((nx-1)*var(x)+ (ny-1)*var(y))/(nx+ny-2))

### return the results
r <- list(est_d = est_d, var_d = var_d,
est_m = est_m, var_m = var_m,
est_v = est_v, var_v = var_v,
cohend = cohend)
return(r)
}

set.seed(1)
results <- replicate(10^3,sim(), simplify = TRUE)

# mean
effect
mean(as.numeric(results[1,]))
mean(as.numeric(results[3,])/as.numeric(results[5,]))
mean(as.numeric(results[7,]))
# variance
mean(as.numeric(results[2,]))
var(as.numeric(results[1,]))

# error
mean((as.numeric(results[1,])-effect)^2)
mean((as.numeric(results[3,])/as.numeric(results[5,])-effect)^2)
mean((as.numeric(results[7,])-effect)^2)