Monte Carlo p-value for linear regression I'm trying to understand a simple monte carlo application so I can create empirical p-values for a linear regression. I'm not sure if I am following the right procedure - I am comparing the slope of my explanatory variable between the real and randomized dependent variable. 
I understand that lmPerm package does what I wish to do in a faster and more complicated way, but I wish to recreate it myself to understand the procedure as I need to scale it up for data with a complicated structure. 
library(lmPerm)
library(dplyr)

states <- as.data.frame(state.x77)
summary(lmp(Murder~ Income+Frost, data=states, perm="Prob", Ca=0, maxIter=100000))
summary(lm(Murder~ Income+Frost, data=states))

iterations = 10000
set.seed(1)
(realSlope <- coef(lm(Murder~ Income+Frost, data=states))["Income"])

randomDataSlopes <- data.frame(row.names = 1:iterations)
for(iteration in 1:iterations) {
  if (iteration %% 1000 == 0) print(iteration)
  #shuffle muder variable
  states <- mutate(states, randomMurder = base::sample(Murder, replace = FALSE))
  model <- lm(randomMurder ~ Income+Frost, states)
  randomDataSlope <- coef(model)["Income"]
  randomDataSlopes[iteration, "slope"] <- randomDataSlope
}

randomDataSlopes <- abs(randomDataSlopes$slope) #convert to absolute values for two sided test (or should it be centered first?)

emerpicalPValue <- (sum(randomDataSlopes >= abs(realSlope)) + 1) / (iterations+1)
paste0("My emperical p-value:", emerpicalPValue) #0.4458
paste0("lmp emperical p-value:", summary(lmp(Murder~ Income+Frost, data=states, perm="Prob", Ca=0, maxIter=100000))$coefficients["Income", "Pr(Prob)"])  #0.3640
paste0("Regular lm/t-stat:", summary(lm(Murder~ Income+Frost, data=states))$coefficients["Income", "Pr(>|t|)"])   #0.366705

Because my p-values seem to be higher than both lmp and lm, I'm wondering if I am doing it right. So I'm wondering if the method I am following is correct, or if I should be shuffling the data differently. 
 A: You just need to run say $B=10,000$ regressions, where during each $b$th iteration $(b=1,2,\ldots,B)$ you randomly permute (only) the $y$-variable and don't do anything to the $x$-variables. (permuting is simply shuffling, or randomly rearranging values in a vector). The p-value for each $j$th coefficient will be the number of times the coefficient's $Z^{(b)}_j=\beta_j/s.e.(\beta_j)$ from all regressions exceeds the same coefficient's single Z-value, $Z_j$, from the unaltered single run of the data, divided by 10,000.
The p-value for the $j$th regression coefficient is:
\begin{equation}
p_j=\frac{\#\{ b: Z^{(b)}_j > Z_j\} }{B}
\end{equation}
So what you are doing is first running the simple model, and note each value of $Z_j$ for each $j$th coefficient.  Then, (a) permute the $y$-values, (b) run the same model, (c) obtaining the $Z^{(b)}_j$ scores for each coefficient from each run -- repeat steps (a)-->(c) 10,000 times.  Next, count the number of times  $Z^{(b)}_j$ exceeds $Z_j$ for each variable, and then divide the count by 10,000. 
You could also shuffle each $x$-variable singly (leaving $y$ alone) and determine its p-value using the same approach, and repeat for each $x$-variable.  This will preserve the correlation between the other $x$-variables while permuting the single $x$-variable of interest during all 10,000 iterations.  
