thr3ads.net - R help - [R] empty plots ! [May 2021]

If this information is useful, please help other people find it:
Share via:

varin sacha

2021-May-12 08:33 UTC

[R] 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"))
>> >
>> > #################
>> >
>> > ______________________________________________
>> > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more,
see
>> > https://stat.ethz.ch/mailman/listinfo/r-help
>> > PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> > and provide commented, minimal, self-contained, reproducible code.
>> >
>>
>> ______________________________________________
>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
> 
> ? ? ? ? [[alternative HTML version deleted]]
> 
> 
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
______________________________________________
R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Jim Lemon

2021-May-12 11:03 UTC

head link

[R] empty plots !

Hi varin,
Were you expecting image files? I don't see any plot device e.g. pdf()
in your code.

Jim

On Wed, May 12, 2021 at 6:34 PM varin sacha via R-help
<r-help at r-project.org> wrote:>
> 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"))
> >> >
> >> > #################
> >> >
> >> > ______________________________________________
> >> > R-help at r-project.org mailing list -- To UNSUBSCRIBE and
more, see
> >> > https://stat.ethz.ch/mailman/listinfo/r-help
> >> > PLEASE do read the posting guide
> >> http://www.R-project.org/posting-guide.html
> >> > and provide commented, minimal, self-contained, reproducible
code.
> >> >
> >>
> >> ______________________________________________
> >> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more,
see
> >> https://stat.ethz.ch/mailman/listinfo/r-help
> >> PLEASE do read the posting guide
> >> http://www.R-project.org/posting-guide.html
> >> and provide commented, minimal, self-contained, reproducible code.
> >>
> >
> >         [[alternative HTML version deleted]]
>
> >
> >
> > ______________________________________________
> > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

Rui Barradas

2021-May-12 12:03 UTC

head link

[R] empty plots !

Hello,

All power.* are filled with NA, please revise the code that produces 
those matrices, there's nothing wrong with the plotting code.

Also, suggested in an answer to your previous e-mail *with* code that 
you should save

old_par <- par(mfrow = c(4,2), cex = 0.45)
[...]
par(old_par)  # at the end, put the graphics device pars back.


Hope this helps,

Rui Barradas

?s 09:33 de 12/05/21, varin sacha via R-help escreveu:> 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"))
>>>>
>>>> #################
>>>>
>>>> ______________________________________________
>>>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and
more, see
>>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>>> and provide commented, minimal, self-contained, reproducible
code.
>>>>
>>>
>>> ______________________________________________
>>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more,
see
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>
>>  ? ? ? ? [[alternative HTML version deleted]]
> 
>>
>>
>> ______________________________________________
>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
> 
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

R help - May 2021 - empty plots !

[R] empty plots !

[R] empty plots !

[R] empty plots !