Commit a01da6f8 by Eric Coissac

temporary partial results on power

parent eadd4706
......@@ -356,44 +356,66 @@ n_rand <- 1000
<<estimate_power, cache=TRUE, message=FALSE, warning=FALSE, include=FALSE, dependson="estimate_power_setting">>=
if (! dir.exists("estimate_power_setting.cache")) {
dir.create("estimate_power_setting.cache")
}
compute <- FALSE
h1_sims = array(0,dim = c(n_sim,length(n_indivduals),length(p_qs),length(r2s),3))
for (k in seq_len(n_sim)) {
filename = sprintf("estimate_power_setting.cache/%d.rdata",k)
if (file.exists(filename)) {
h1_sims[k, , , , ] <- get(load(filename))
} else {
for (i in seq_along(n_indivduals))
for (j in seq_along(p_qs))
for (r in seq_along(r2s)) {
X <- simulate_matrix(n_indivduals[i],
p_qs[j],
equal.var = TRUE)
Y <- simulate_correlation(X,
p_qs[j],
r2 = r2s[r],
equal.var = TRUE)
h1_sims[k,i,j,r,1] <- attr(corls(X,Y,nrand = n_rand),"p.value")[1,2]
h1_sims[k,i,j,r,2] <- suppressWarnings(vegan::protest(X,Y,permutations = n_rand)$signif)
h1_sims[k,i,j,r,3] <- ade4::procuste.rtest(as.data.frame(X),
as.data.frame(Y),
nrepet = n_rand)$pvalue
}
res <- h1_sims[k, , , , ]
save(res,file = filename)
if (compute) {
if (!dir.exists("estimate_power_setting.cache")) {
dir.create("estimate_power_setting.cache")
}
}
h1_sims = array(0,dim = c(n_sim,length(n_indivduals),length(p_qs),length(r2s),3))
for (k in seq_len(n_sim)) {
filename = sprintf("estimate_power_setting.cache/%d.rdata",k)
if (file.exists(filename)) {
h1_sims[k, , , , ] <- get(load(filename))
} else {
for (i in seq_along(n_indivduals))
for (j in seq_along(p_qs))
for (r in seq_along(r2s)) {
X <- simulate_matrix(n_indivduals[i],
p_qs[j],
equal.var = TRUE)
Y <- simulate_correlation(X,
p_qs[j],
r2 = r2s[r],
equal.var = TRUE)
h1_sims[k,i,j,r,1] <- attr(corls(X,Y,nrand = n_rand),"p.value")[1,2]
h1_sims[k,i,j,r,2] <- suppressWarnings(vegan::protest(X,Y,permutations = n_rand)$signif)
h1_sims[k,i,j,r,3] <- ade4::procuste.rtest(as.data.frame(X),
as.data.frame(Y),
nrepet = n_rand)$pvalue
}
res <- h1_sims[k, , , , ]
save(res,file = filename)
}
unlink("estimate_power_setting.cache",
recursive = TRUE)
}
} else {
read_cache = function() {
h1_sims = array(0,dim = c(n_sim,
length(n_indivduals),
length(p_qs),
length(r2s),3))
for (k in seq_len(n_sim)) {
filename = sprintf("manuscript/estimate_power_setting.cache/%d.rdata",k)
if (file.exists(filename)) {
h1_sims[k, , , , ] <- get(load(filename))
}
}
h1_sims
}
h1_sims = read_cache()
}
#unlink("estimate_power_setting.cache",
# recursive = TRUE)
@
<<estimate_power_complementary, message=FALSE, warning=FALSE, include=FALSE>>
......
No preview for this file type
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment