First, load the Daggity library and define a DAG. We use the following variables: * X: the treatment of AR vs. images. * M: the mediator of choices. * Y: the outcome of product returns. * U: the unobserved confounders.

library(dagitty)

g <- dagitty('dag {
    M [pos="2,-0.5"]
    U [pos="3,-1"]
    X [pos="1,0"]
    Y [pos="4,0"]
    
    X -> M -> Y
    M <- U -> Y
    X -> Y
}')

plot(g)

From the plot, we cannot identify the mediation path due to confounding effect, but the total effect of \(X \rightarrow Y\) is still identifiable.

adjustmentSets(g, exposure = "X", outcome = "Y",
               effect = "total", type = "minimal")
##  {}

The adjustment set is an empty set, so all non-causal paths are blocked. To further verify this, I ran a simple simulation based on the DAG with a linear model:

myfun <- function(b_xm, b_my, b_xy) {
  
  n <- 5000
  x <- rnorm(n)
  u <- rnorm(n)
  m <- b_xm*x + 0.5*u + rnorm(n)
  y <- b_xy*x + b_my*m + u + rnorm(n)
  
  data <- data.frame(x,y,m)
  lm_xy <- lm(y~0+x,data)
  
  out <- c(b_xy+b_xm*b_my, lm_xy$coefficients[1])
  names(out) <- c("True Total Effects", "Estimated Total Effects")
  
  return(out)
  
}

You can set the path values to any values, run the function and check:

mypara <- rnorm(3)
myfun(mypara[1],mypara[2],mypara[3])
##      True Total Effects Estimated Total Effects 
##                2.365445                2.322361