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.


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
(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.


1 Answer 1


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.

  • $\begingroup$ Thank you, I think I see - I should be running sample on 'Income' (y-variable of interest) and not Murder (x-variable). When I do that, the p-values agree. Do you know how I would treat interactions of say two y-variables? shuffle both in the same way? $\endgroup$
    – leonfrench
    Mar 9, 2017 at 15:59
  • $\begingroup$ I've worked out how lmp does the interactions - the two variables are shuffled/sampled in the same way, then included as an interaction, for example: model <- lm(Murder ~ Frost + Illiteracy + randomFrost * randomIlliteracy, states) $\endgroup$
    – leonfrench
    Mar 9, 2017 at 18:13
  • $\begingroup$ For interactions between x-variables (y-variables don't have interactions), you should probably only shuffle one of the two variables, because you need the main effect of each variable in the model as well the interaction term. Shuffling one of the two interacting variables will break up the interaction, at the same time forcing the null hypothesis on the shuffled variable by itself (main effect). $\endgroup$
    – user32398
    Mar 9, 2017 at 20:41
  • $\begingroup$ I see your point, but doing the above - shuffling both variables in the same way/order and including the unshuffled versions as main effects gave me p-values that matched the lmp results. Which is pretty close to your suggestion. $\endgroup$
    – leonfrench
    Mar 10, 2017 at 17:44
  • $\begingroup$ I meant shuffle one of the main effects before you calculate the interaction term, and then this single variable is shuffled but will affect both the main effect and its counterpart in the interaction. You could just shuffle the interaction term alone -- and it sounds like that's what you did, which is understandable, but the interaction term is affected by each variables main effect, so you'll need to take that into consideration. Recall, when you're looking at interaction terms and their significance, the main effect does not need to be significant it just needs to be in the model. $\endgroup$
    – user32398
    Mar 10, 2017 at 17:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.