z <- indirect/sqrt(b^2*sa^2 + a^2*sb^2 + sa^2*sb^2)
pval <- 2*pnorm(-abs(z));

return(c(indirect = indirect, z = z, pval = pval))

indirect <- a*b
indirect

indirect <- a*b
indirect

## resample data n times
resample <- data[i,]

## refit regression model on resampled data
model.M.resample <- lm(formula(model.M), data = resample)
model.Y.resample <- lm(formula(model.Y), data = resample)

a <- coef(model.M.resample)[apath] ## a path
b <- coef(model.Y.resample)[bpath] ## b path
c <- coef(model.Y.resample)[apath] ## c path

## effect decomposition
## indirect effect
indirect.boot <- a*b
direct.boot <- c
return(c(direct.boot = direct.boot, indirect.boot = indirect.boot))

## resample data based on index i
resample.data <- data[i, ]

## re-estimate the regression models
model.m1.resample <- lm(formula(model.m1), data = resample.data)
model.m2.resample <- lm(formula(model.m2), data = resample.data)
model.Y.resample <- lm(formula(model.Y), data = resample.data)

## extract coefficients and compute indirect effects
a1 <- coef(model.m1.resample)["cond"]
a2 <- coef(model.m2.resample)["cond"]
b1 <- coef(model.Y.resample)["pmi"]
b2 <- coef(model.Y.resample)["import"]
c <- coef(model.Y.resample)["cond"]

indirect_pmi <- a1*b1
indirect_import <- a2*b2
direct <- c

## output vector, name - values
out <- c(indirect_pmi = indirect_pmi,
         indirect_import = indirect_import,
         total_indirect = indirect_pmi + indirect_import,
         direct = direct)
return(out)

## resample data
resample.data <- data[i, ]

## re-estimate model of M and Y based on resampled data
multicat.M.rspl <- lm(formula(multicat.model.M), data = resample.data)
multicat.Y.rspl <- lm(formula(multicat.model.Y), data = resample.data)

## recover coefficients
a1 <- coef(multicat.M.rspl)["d1"]
a2 <- coef(multicat.M.rspl)["d2"]
b <- coef(multicat.Y.rspl)["inter"]
cp1 <- coef(multicat.Y.rspl)["d1"]
cp2 <- coef(multicat.Y.rspl)["d2"]

## effect decomposition
relative.ind.d1 <- a1*b
relative.ind.d2 <- a2*b
relative.direct.d1 <- cp1
relative.direct.d2 <- cp2

out <- c(relative.ind.d1 = relative.ind.d1,
         relative.ind.d2 = relative.ind.d2,
         relative.direct.d1 = relative.direct.d1,
         relative.direct.d2 = relative.direct.d2)

return(out)

## changes in R-sq (as a whole)
summary(multicat.model)$r.squared - summary(multicat.model.null)$r.squared
## anova (whether R-sq change is significant?)
anova(multicat.model, multicat.model.null)

## changes in R-sq (as a whole)
summary(multicat.model.Y)$r.squared - summary(multicat.model.null2)$r.squared
## anova (whether R-sq change is significant?)
anova(multicat.model.null2, multicat.model.Y)

## resample data
resample.data <- data[i, ]

## re-estimate model of M and Y based on resampled data
multicat.M.null.rspl <- lm(inter ~ 1, data = resample.data)
multicat.M.rspl <- lm(formula(multicat.model.M), data = resample.data)
multicat.Y.rspl <- lm(formula(multicat.model.Y), data = resample.data)

## recover quantity of interest
rsq.change <- summary(multicat.M.rspl)$r.squared - summary(multicat.M.null.rspl)$r.squared
b <- coef(multicat.Y.rspl)["inter"]

omnibus.indirect.multicat <- rsq.change*b

return(omnibus.indirect.multicat)