CODEWITHSUNDEEP
r
Fuv8
Fuv8

Reputation: 905

Cycles of permutation tests

I have to do the following job: generate a random matrix, apply a linear model, then shuffle the matrix, apply the linear model, then shuffle again the matrix and apply a linear model,...., for 20 times and for each time I have to save the p-value. This job has to be done for 1000 times.

I wrote the following piece of code, but I'm not able to run a for cycle in a for cycle. Here what I wrote: 

B=1000   
n=100
b =20   
my.seed=1    
my.intercept<-0     
my.slope<-1      
res <- data.frame(matrix(ncol = 4, nrow = B))     
colnames(res) <- c("Estimate", "St_Err", "t_val", "P_val")     

for (i in 1:B)     
    for (j in 1:b){
        x1=rbinom(n, 1, 0.5)     
        e=rnorm(n, 1, 1)       
        my_model=lm(y~x1)      
        y <- true.intercept + true.slope*x1+e     
        res[i,] <- data.frame(summary(lm(y ~x1))$coefficients)
    }
}

I don't know how to save the results of the for cycle on j and then save the p values of the 20 permutations on the total 1000 permutations in order at the end to have finally a data.frame with 20 rows (because of the initial matrix is shuffled 20 times) and 1000 columns (because the permutations are 1000)

Upvotes: 2

Views: 104

Answers (1)

Carl Witthoft
Carl Witthoft

Reputation: 21502

So you want an output matrix p_vals[i,j] ? Then, instead of res[i,]<- stuff, you can steal some code from the print method of summary.lm :

[print " p-value:",]  format.pval(pf(x$fstatistic[1L], 
                x$fstatistic[2L], x$fstatistic[3L], lower.tail = FALSE)

where pf is grabbing data from summary(lm(etc))$fstatistic

so you would do 

fstats<- summary(lm(y ~x1))$fstatistic
pvals[i,j]<- pf(fstats[1L], fstats[2L], fstats[3L], lower.tail = FALSE)

Upvotes: 1

Related Questions