empty plots !

Dear Experts,

My R code was perfectly working since I decide to add a 5th correlation coefficient : hoeffdings' D.
fter a google search, I guess I need somewhere in my R code "unlist" but I don't know where !
Here below my R code with 1 error message. At the end I get my 8 plots but they are empty !
Many thanks for your precious help !

#################
set.seed(1)
library(energy)
library(independence)
library(TauStar)

# Here we define parameters which we use to simulate the data
# The number of null datasets we use to estimate our rejection reject #regions for an alternative with level 0.05
nsim=50

# Number of alternative datasets we use to estimate our power
nsim2=50

# The number of different noise levels used
num.noise <- 30

# A constant to determine the amount of noise
noise <- 3

# Number of data points per simulation

n=100

# Vectors holding the null "correlations" (for pearson, for spearman, for #kendall, for hoeffding and dcor respectively) for each of the nsim null datasets at a #given noise level
val.cor=val.cors=val.cork=val.dcor=val.hoe=rep(NA,nsim)

# Vectors holding the alternative "correlations" (for pearson, for #spearman, for kendall, for hoeffding and dcor respectively) for each of #the nsim2 #alternative datasets at a given noise level
val.cor2=val.cors2=val.cork2=val.dcor2=val.hoe2= rep(NA,nsim2)

# Arrays holding the estimated power for each of the 4 "correlation" types, #for each data type (linear, parabolic, etc...) with each noise level
power.cor=power.cors=power.cork=power.dcor=power.hoe= array(NA, c(8,num.noise))

## We loop through the noise level and functional form; each time we #estimate a null distribution based on the marginals of the data, and then #use that null distribution to estimate power
## We use a uniformly distributed x, because in the original paper the #authors used the same

for(l in 1:num.noise){

      for(typ in 1:8){

## This next loop simulates data under the null with the correct marginals #(x is uniform, and y is a function of a uniform with gaussian noise)

    for(ii in 1:nsim){
      x=runif(n)

#lin+noise
if(typ==1){
y=x+ noise *(l/num.noise)* rnorm(n)
}

#parabolic+noise
if(typ==2){
y=4*(x-.5)^2+  noise * (l/num.noise) * rnorm(n)
}

#cubic+noise
if(typ==3){
y=128*(x-1/3)^3-48*(x-1/3)^3-12*(x-1/3)+10* noise  * (l/num.noise) *rnorm(n)
}

#sin+noise
if(typ==4){
y=sin(4*pi*x) + 2*noise * (l/num.noise) *rnorm(n)
}

#their sine + noise
if(typ==5){
y=sin(16*pi*x) + noise * (l/num.noise) *rnorm(n)
}

#x^(1/4) + noise
if(typ==6){
y=x^(1/4) + noise * (l/num.noise) *rnorm(n)
}

#circle
if(typ==7){
y=(2*rbinom(n,1,0.5)-1) * (sqrt(1 - (2*x - 1)^2)) + noise/4*l/num.noise *rnorm(n)
}

#step function
if(typ==8){
y = (x > 0.5) + noise*5*l/num.noise *rnorm(n)
}

# We resimulate x so that we have the null scenario
x <- runif(n)

# Calculate the 5 correlations
val.cor[ii]=(cor(x,y))
val.cors[ii]=(cor(x,y,method=c("spearman")))
val.cork[ii]=(cor(x,y,method=c("kendal")))
val.dcor[ii]=dcor(x,y)
val.hoe[ii]=(hoeffding.D.test(x,y,na.rm=TRUE,collisions=TRUE))
}

## Next we calculate our 5 rejection cutoffs
cut.cor=quantile(val.cor,.95)
cut.cors=quantile(val.cors,.95)
cut.cork=quantile(val.cork,.95)
cut.dcor=quantile(val.dcor,.95)
cut.hoe=quantile(val.hoe,.95)

## Next we simulate the data again, this time under the alternative

    for(ii in 1:nsim2){
      x=runif(n)

#lin+noise
if(typ==1){
y=x+ noise *(l/num.noise)* rnorm(n)
}

#parabolic+noise
if(typ==2){
y=4*(x-.5)^2+  noise * (l/num.noise) * rnorm(n)
}

#cubic+noise
if(typ==3){
y=128*(x-1/3)^3-48*(x-1/3)^3-12*(x-1/3)+10* noise  * (l/num.noise) *rnorm(n)
}

#sin+noise
if(typ==4){
y=sin(4*pi*x) + 2*noise * (l/num.noise) *rnorm(n)
}

#their sine + noise
if(typ==5){
y=sin(16*pi*x) + noise * (l/num.noise) *rnorm(n)
}

#x^(1/4) + noise
if(typ==6){
y=x^(1/4) + noise * (l/num.noise) *rnorm(n)
}

#circle
if(typ==7){
y=(2*rbinom(n,1,0.5)-1) * (sqrt(1 - (2*x - 1)^2)) + noise/4*l/num.noise *rnorm(n)
}

#step function
if(typ==8){
y = (x > 0.5) + noise*5*l/num.noise *rnorm(n)
}

## We again calculate our 5 correlations
val.cor2[ii]=(cor(x,y))
val.cors2[ii]=(cor(x,y,method=c("spearman")))
val.cork2[ii]=(cor(x,y,method=c("kendal")))
val.dcor2[ii]=dcor(x,y)
val.hoe2[ii]=(hoeffding.D.test(x,y,na.rm=TRUE,collisions=TRUE))
}

## Now we estimate the power as the number of alternative statistics #exceeding our estimated cutoffs
power.cor[typ,l] <- sum(val.cor2 > cut.cor)/nsim2
power.cors[typ,l] <- sum(val.cors2 > cut.cor)/nsim2
power.cork[typ,l] <- sum(val.cork2 > cut.cor)/nsim2
power.dcor[typ,l] <- sum(val.dcor2 > cut.dcor)/nsim2
power.hoe[typ,l] <- sum(val.hoe2 > cut.hoe)/nsim2
}
}

## The rest of the code is for plotting the image
par(mfrow = c(4,2), cex = 0.45)
plot((1:30)/10, power.cor[1,], ylim = c(0,1), main = "Linear", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b')
points((1:30)/10, power.cors[1,], pch = 2, col = "green", type = 'b')
points((1:30)/10, power.cork[1,], pch = 3, col = "blue", type = 'b')
points((1:30)/10, power.dcor[1,], pch = 4, col = "red", type = 'b')
points((1:30)/10, power.hoe[1,], pch = 5, col = "purple", type = 'b')
legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor","hoe" ), pch = c(1,2,3,4,5), col = c("black","green","blue","red", "purple"))

plot((1:30)/10, power.cor[2,], ylim = c(0,1), main = "Quadratic", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b')
points((1:30)/10, power.cors[2,], pch = 2, col = "green", type = 'b')
points((1:30)/10, power.cork[2,], pch = 3, col = "blue", type = 'b')
points((1:30)/10, power.dcor[2,], pch = 4, col = "red", type = 'b')
points((1:30)/10, power.hoe[2,], pch = 5, col = "purple", type = 'b')
legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor","hoe" ), pch = c(1,2,3,4,5), col = c("black","green","blue","red", "purple"))

plot((1:30)/10, power.cor[3,], ylim = c(0,1), main = "Cubic", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b')
points((1:30)/10, power.cors[3,], pch = 2, col = "green", type = 'b')
points((1:30)/10, power.cork[3,], pch = 3, col = "blue", type = 'b')
points((1:30)/10, power.dcor[3,], pch = 4, col = "red", type = 'b')
points((1:30)/10, power.hoe[3,], pch = 5, col = "purple", type = 'b')
legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor","hoe" ), pch = c(1,2,3,4,5), col = c("black","green","blue","red", "purple"))

plot((1:30)/10, power.cor[5,], ylim = c(0,1), main = "Sine: period 1/8", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b')
points((1:30)/10, power.cors[5,], pch = 2, col = "green", type = 'b')
points((1:30)/10, power.cork[5,], pch = 3, col = "blue", type = 'b')
points((1:30)/10, power.dcor[5,], pch = 4, col = "red", type = 'b')
points((1:30)/10, power.hoe[5,], pch = 5, col = "purple", type = 'b')
legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor","hoe" ), pch = c(1,2,3,4,5), col = c("black","green","blue","red", "purple"))

plot((1:30)/10, power.cor[4,], ylim = c(0,1), main = "Sine: period 1/2", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b')
points((1:30)/10, power.cors[4,], pch = 2, col = "green", type = 'b')
points((1:30)/10, power.cork[4,], pch = 3, col = "blue", type = 'b')
points((1:30)/10, power.dcor[4,], pch = 4, col = "red", type = 'b')
points((1:30)/10, power.hoe[4,], pch = 5, col = "purple", type = 'b')
legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor","hoe" ), pch = c(1,2,3,4,5), col = c("black","green","blue","red", "purple"))

plot((1:30)/10, power.cor[6,], ylim = c(0,1), main = "X^(1/4)", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b')
points((1:30)/10, power.cors[6,], pch = 2, col = "green", type = 'b')
points((1:30)/10, power.cork[6,], pch = 3, col = "blue", type = 'b')
points((1:30)/10, power.dcor[6,], pch = 4, col = "red", type = 'b')
points((1:30)/10, power.hoe[6,], pch = 5, col = "purple", type = 'b')
legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor","hoe" ), pch = c(1,2,3,4,5), col = c("black","green","blue","red", "purple"))

plot((1:30)/10, power.cor[7,], ylim = c(0,1), main = "Circle", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b')
points((1:30)/10, power.cors[7,], pch = 2, col = "green", type = 'b')
points((1:30)/10, power.cork[7,], pch = 3, col = "blue", type = 'b')
points((1:30)/10, power.dcor[7,], pch = 4, col = "red", type = 'b')
points((1:30)/10, power.hoe[7,], pch = 5, col = "purple", type = 'b')
legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor","hoe" ), pch = c(1,2,3,4,5), col = c("black","green","blue","red", "purple"))

plot((1:30)/10, power.cor[8,], ylim = c(0,1), main = "Step function", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b')
points((1:30)/10, power.cors[8,], pch = 2, col = "green", type = 'b')
points((1:30)/10, power.cork[8,], pch = 3, col = "blue", type = 'b')
points((1:30)/10, power.dcor[8,], pch = 4, col = "red", type = 'b')
points((1:30)/10, power.hoe[8,], pch = 5, col = "purple", type = 'b')
legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor","hoe" ), pch = c(1,2,3,4,5), col = c("black","green","blue","red", "purple"))
#################









Le mardi 11 mai 2021 ? 20:00:49 UTC+2, varin sacha via R-help <r-help at r-project.org> a ?crit :





Dear all,

Many thanks for your responses.

Best
S.







Le lundi 10 mai 2021 ? 17:18:59 UTC+2, Bill Dunlap <williamwdunlap at gmail.com> a ?crit :





Also, normalizePath("power.pdf").

On Sun, May 9, 2021 at 5:13 PM Bert Gunter <bgunter.4567 at gmail.com> wrote:

?getwd

Bert Gunter

"The trouble with having an open mind is that people keep coming along and
sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )


On Sun, May 9, 2021 at 2:59 PM varin sacha via R-help <r-help at r-project.org>
wrote:

Rui,

The created pdf.file is off-screen device. Indeed after dev.off() I should
view the pdf file on my computer. But I don't find it. Where do I find the
pdf.file ?

Regards,



Le dimanche 9 mai 2021 ? 22:44:22 UTC+2, Rui Barradas <
ruipbarradas at sapo.pt> a ?crit :





Hello,

You are not closing the pdf device.
The only changes I have made to your code are right at the beginning of
the plotting instructions and at the end of the code.


## The rest of the code is for plotting the image
pdf(file = "power.pdf")
op <- par(mfrow = c(4,2), cex = 0.45)

[...]

par(op)
dev.off()
#################

The comments only line is your last code line.
The result is attached.

Hope this helps,

Rui Barradas

?s 19:39 de 09/05/21, varin sacha via R-help escreveu:

Dear R-experts,

I am trying to get the 8 graphs like the ones in this paper :
https://statweb.stanford.edu/~tibs/reshef/comment.pdf
My R code does not show any error message neither warnings but I d'on't

get what I would like to get (I mean the 8 graphs), so I am missing
something. What's it ? Many thanks for your precious help.

#################
set.seed(1)
library(energy)

# Here we define parameters which we use to simulate the data
# The number of null datasets we use to estimate our rejection reject

#regions for an alternative with level 0.05

nsim=50

# Number of alternative datasets we use to estimate our power
nsim2=50

# The number of different noise levels used
num.noise <- 30

# A constant to determine the amount of noise
noise <- 3

# Number of data points per simulation
n=100

# Vectors holding the null "correlations" (for pearson, for spearman,

for kendall and dcor respectively) for each # of the nsim null datasets at
a #given noise level

val.cor=val.cors=val.cork=val.dcor=rep(NA,nsim)

# Vectors holding the alternative "correlations" (for pearson, for

#spearman, for kendall and dcor respectively) #for each of the nsim2
alternative datasets at a given noise level

val.cor2=val.cors2=val.cork2=val.dcor2= rep(NA,nsim2)


# Arrays holding the estimated power for each of the 4 "correlation"

types, for each data type (linear, #parabolic, etc...) with each noise level

power.cor=power.cors=power.cork=power.dcor= array(NA, c(8,num.noise))

## We loop through the noise level and functional form; each time we

#estimate a null distribution based on #the marginals of the data, and then
#use that null distribution to estimate power

## We use a uniformly distributed x, because in the original paper the

#authors used the same

for(l in 1:num.noise) {

       for(typ in 1:8) {

## This next loop simulates data under the null with the correct

marginals (x is uniform, and y is a function of a #uniform with gaussian
noise)

     for(ii in 1:nsim) {
       x=runif(n)

#lin+noise
if(typ==1) {
y=x+ noise *(l/num.noise)* rnorm(n)
}

#parabolic+noise
if(typ==2) {
y=4*(x-.5)^2+  noise * (l/num.noise) * rnorm(n)
}

#cubic+noise
if(typ==3) {
y=128*(x-1/3)^3-48*(x-1/3)^3-12*(x-1/3)+10* noise  * (l/num.noise)

*rnorm(n)

}

#sin+noise
if(typ==4) {
y=sin(4*pi*x) + 2*noise * (l/num.noise) *rnorm(n)
}

#their sine + noise
if(typ==5) {
y=sin(16*pi*x) + noise * (l/num.noise) *rnorm(n)
}

#x^(1/4) + noise
if(typ==6) {
y=x^(1/4) + noise * (l/num.noise) *rnorm(n)
}

#circle
if(typ==7) {
y=(2*rbinom(n,1,0.5)-1) * (sqrt(1 - (2*x - 1)^2)) + noise/4*l/num.noise

*rnorm(n)

}

#step function
if(typ==8) {
y = (x > 0.5) + noise*5*l/num.noise *rnorm(n)
}

# We resimulate x so that we have the null scenario
x <- runif(n)

# Calculate the 4 correlations
val.cor[ii]=(cor(x,y))
val.cors[ii]=(cor(x,y,method=c("spearman")))
val.cork[ii]=(cor(x,y,method=c("kendal")))
val.dcor[ii]=dcor(x,y)
}

## Next we calculate our 4 rejection cutoffs
cut.cor=quantile(val.cor,.95)
cut.cors=quantile(val.cors,.95)
cut.cork=quantile(val.cork,.95)
cut.dcor=quantile(val.dcor,.95)

## Next we simulate the data again, this time under the alternative

     for(ii in 1:nsim2) {
       x=runif(n)

#lin+noise
if(typ==1) {
y=x+ noise *(l/num.noise)* rnorm(n)
}

#parabolic+noise
if(typ==2) {
y=4*(x-.5)^2+  noise * (l/num.noise) * rnorm(n)
}

#cubic+noise
if(typ==3) {
y=128*(x-1/3)^3-48*(x-1/3)^3-12*(x-1/3)+10* noise  * (l/num.noise)

*rnorm(n)

}

#sin+noise
if(typ==4) {
y=sin(4*pi*x) + 2*noise * (l/num.noise) *rnorm(n)
}

#their sine + noise
if(typ==5) {
y=sin(16*pi*x) + noise * (l/num.noise) *rnorm(n)
}

#x^(1/4) + noise
if(typ==6) {
y=x^(1/4) + noise * (l/num.noise) *rnorm(n)
}

#circle
if(typ==7) {
y=(2*rbinom(n,1,0.5)-1) * (sqrt(1 - (2*x - 1)^2)) + noise/4*l/num.noise

*rnorm(n)

}

#step function
if(typ==8) {
y = (x > 0.5) + noise*5*l/num.noise *rnorm(n)
}

## We again calculate our 4 "correlations"
val.cor2[ii]=(cor(x,y))
val.cors2[ii]=(cor(x,y,method=c("spearman")))
val.cork2[ii]=(cor(x,y,method=c("kendal")))
val.dcor2[ii]=dcor(x,y)
}

## Now we estimate the power as the number of alternative statistics

#exceeding our estimated cutoffs

power.cor[typ,l] <- sum(val.cor2 > cut.cor)/nsim2
power.cors[typ,l] <- sum(val.cors2 > cut.cor)/nsim2
power.cork[typ,l] <- sum(val.cork2 > cut.cor)/nsim2
power.dcor[typ,l] <- sum(val.dcor2 > cut.dcor)/nsim2
}
}

save.image()

## The rest of the code is for plotting the image
pdf("power.pdf")
par(mfrow = c(4,2), cex = 0.45)
plot((1:30)/10, power.cor[1,], ylim = c(0,1), main = "Linear", xlab =

"Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b')

points((1:30)/10, power.cors[1,], pch = 2, col = "green", type = 'b')
points((1:30)/10, power.cork[1,], pch = 3, col = "blue", type = 'b')
points((1:30)/10, power.dcor[1,], pch = 4, col = "red", type = 'b')
legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"),

pch = c(1,2,3), col = c("black","green","blue","red"))

plot((1:30)/10, power.cor[2,], ylim = c(0,1), main = "Quadratic", xlab =

"Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b')

points((1:30)/10, power.cors[2,], pch = 2, col = "green", type = 'b')
points((1:30)/10, power.cork[2,], pch = 3, col = "blue", type = 'b')
points((1:30)/10, power.dcor[2,], pch = 4, col = "red", type = 'b')
legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"),

pch = c(1,2,3), col = c("black","green","blue","red"))

plot((1:30)/10, power.cor[3,], ylim = c(0,1), main = "Cubic", xlab =

"Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b')

points((1:30)/10, power.cors[3,], pch = 2, col = "green", type = 'b')
points((1:30)/10, power.cork[3,], pch = 3, col = "blue", type = 'b')
points((1:30)/10, power.dcor[3,], pch = 4, col = "red", type = 'b')
legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"),

pch = c(1,2,3), col = c("black","green","blue","red"))

plot((1:30)/10, power.cor[5,], ylim = c(0,1), main = "Sine: period 1/8",

xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b')

points((1:30)/10, power.cors[5,], pch = 2, col = "green", type = 'b')
points((1:30)/10, power.cork[5,], pch = 3, col = "blue", type = 'b')
points((1:30)/10, power.dcor[5,], pch = 4, col = "red", type = 'b')
legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"),

pch = c(1,2,3), col = c("black","green","blue","red"))

plot((1:30)/10, power.cor[4,], ylim = c(0,1), main = "Sine: period 1/2",

xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b')

points((1:30)/10, power.cors[4,], pch = 2, col = "green", type = 'b')
points((1:30)/10, power.cork[4,], pch = 3, col = "blue", type = 'b')
points((1:30)/10, power.dcor[4,], pch = 4, col = "red", type = 'b')
legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"),

pch = c(1,2,3), col = c("black","green","blue","red"))

plot((1:30)/10, power.cor[6,], ylim = c(0,1), main = "X^(1/4)", xlab =

"Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b')

points((1:30)/10, power.cors[6,], pch = 2, col = "green", type = 'b')
points((1:30)/10, power.cork[6,], pch = 3, col = "blue", type = 'b')
points((1:30)/10, power.dcor[6,], pch = 4, col = "red", type = 'b')
legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"),

pch = c(1,2,3), col = c("black","green","blue","red"))

plot((1:30)/10, power.cor[7,], ylim = c(0,1), main = "Circle", xlab =

"Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b')

points((1:30)/10, power.cors[7,], pch = 2, col = "green", type = 'b')
points((1:30)/10, power.cork[7,], pch = 3, col = "blue", type = 'b')
points((1:30)/10, power.dcor[7,], pch = 4, col = "red", type = 'b')
legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"),

pch = c(1,2,3), col = c("black","green","blue","red"))

plot((1:30)/10, power.cor[8,], ylim = c(0,1), main = "Step function",

xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b')

points((1:30)/10, power.cors[8,], pch = 2, col = "green", type = 'b')
points((1:30)/10, power.cork[8,], pch = 3, col = "blue", type = 'b')
points((1:30)/10, power.dcor[8,], pch = 4, col = "red", type = 'b')
legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"),

pch = c(1,2,3), col = c("black","green","blue","red"))

#################

        [[alternative HTML version deleted]]