Matrix to calculate effects using delta method
Aren't you looking for a simple matrix multiplication? Try:

x &lt;- as.matrix(dpcp[c("x1", "x2", "x3", "x4")])


Then

marg_effects &lt;- x %*% t(m4)


or equivalently (maybe a little faster):

marg_effects &lt;- tcrossprod(x, m4)