Matrix to calculate effects using delta method I have written some code to calculate effects using the Delta method.

I have a dataframe dpcp with variables x1, x2, x3, x4 and a matrix of 1000 draws from a multivariate normal, m4[1000,4]. 

This code works to calculate the effects, but it runs very slowly - it currently takes upwards of 5 hours for only 2000 observations.  How can I improve the runtime?

n=nrow(dpcp)
for (i in 1:n)

 {

     for (j in 1:1000)
      {        
     marg_effects[i,j]=(m4[j,1]*dpcp[i,]$x1)+(m4[j,2]*dpcp[i,]$x2)+(m4[j,3]*dpcp[i,]$x3)+(m4[j,4]*dpcp[i,]$x4)          
 } 
   }
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)