Skip to content

Commit

Permalink
lots o fiddling with plot options to look proper at small resolution
Browse files Browse the repository at this point in the history
  • Loading branch information
@cboettig
cboettig committed Apr 6, 2011
1 parent 6b253db commit e1c77f9
Show file tree
Hide file tree
Showing 5 changed files with 61 additions and 35 deletions.
15 changes: 8 additions & 7 deletions warningsignals/R/bootstrap_indicators.R
View file
Expand Up @@ -63,7 +63,7 @@ plot.bootstrap_tau <- function(taus, show_p = FALSE, show_error=TRUE, ...){
m <- length(taus[[1]]) ## number of indicators

## set up m x n plot-matrix, no margins on subplots, add outer margins
par(mfrow=c(m,n), oma=c(8,8,8,4), mar=c(0,0,0,0))
par(mfrow=c(m,n), oma=c(3,4,3,2), mar=c(0,0,0,0), ...)

for(j in 1:m){
for(i in 1:n){
Expand All @@ -72,15 +72,16 @@ plot.bootstrap_tau <- function(taus, show_p = FALSE, show_error=TRUE, ...){
if(i > 1){ yaxt <- "n"
} else { yaxt <- "s" }

plot(taus[[i]][[j]], show_p = show_p, show_error=show_error, xaxt = xaxt, yaxt=yaxt, ...)

if(j==1) mtext(data_names[i], NORTH<-3, cex=par()$cex.lab, line=2) ## data labels on top row
plot(taus[[i]][[j]], show_p = show_p, show_error=show_error,
xaxt = xaxt, yaxt=yaxt, ...)
if(j==1) mtext(data_names[i], NORTH<-3,
cex=par()$cex.lab, line=1) ## data labels on top
if(i==1){
mtext(taus[[i]][[j]]$signal, WEST<-2, line=4) ## statistic name on first column
mtext(taus[[i]][[j]]$signal, WEST<-2, line=3, cex=par()$cex.lab) ## statistic name on first column
mtext(expression(paste("Prob Density of ", tau)),
WEST<-2, line=2, cex=.7*par()$cex.lab) ## statistic name on first column
WEST<-2, line=2, cex=.6*par()$cex.lab) ## statistic name
}
if(j==m & i==2) mtext(expression(paste(tau)), SOUTH<-1, line=4, cex=par()$cex.lab) ## x-axis label
if(j==m & i==2) mtext(expression(paste(tau)), SOUTH<-1, line=3, cex=par()$cex.lab) ## x-axis label
}
}
}
Expand Down
45 changes: 33 additions & 12 deletions warningsignals/R/indicators.R
View file
Expand Up @@ -102,16 +102,37 @@ plot_indicator <- function(X, indicator=c("Autocorrelation", "Variance", "Skew",
}
Y <- compute_indicator(X, indicator, windowsize)
time_window <- time(X)[windowsize:length(X)]
plot(time_window, Y, xlim=c(start(X)[1], end(X)[1]), type="l", xlab="time", ylab=indicator, lwd=2, ...)
plot(time_window, Y, xlim=c(start(X)[1], end(X)[1]), type="l", xlab="time", ylab=indicator, ...)
abline(v=time_window[1], lty="dashed")

out <- cor.test(time(X)[windowsize:length(X)], Y, method=method)

w <- c(out$estimate, out$p.value)
text(xshift(xpos), yshift(ypos),
#substitute(paste(method, "coef = ", val, " (p ", pval, ")"), list(val=round(w[1],2),pval=format.pval(w[2]))),
paste(method, "coef=", round(w[1],2), "pval=", format.pval(w[2]) ), pos=4, cex=.9*par()$cex.lab
)
if (method=="kendall"){
## newline character doesn't work
# text(xshift(xpos), yshift(ypos),
# substitute(paste(tau == val, "\n (p== ", pval, ")"), list(val=round(w[1],2),pval=format.pval(w[2]))),
# pos=4, cex=par()$cex.lab)
text(xshift(0), yshift(ypos),
substitute(paste(tau == val),
list(val=round(w[1],2) )),
pos=4, cex=par()$cex.axis)
text(xshift(0), yshift(ypos-20),
substitute(paste("(",p == pval,")"),
list(pval=format.pval(w[2], digits=2))),
pos=4, cex=par()$cex.axis)
# legend("topleft", substitute(tau == val, list(val=round(w[1],2) )))
}
else if (method=="pearson") {
text(xshift(xpos), yshift(ypos),
substitute(paste(r == val, "\n (p== ", pval, ")"), list(val=round(w[1],2),pval=format.pval(w[2]))),
pos=4, cex=par()$cex.lab)
}
else if (method=="spearman") {
text(xshift(xpos), yshift(ypos),
substitute(paste(rho == val, " (p== ", pval, ")"), list(val=round(w[1],2),pval=format.pval(w[2]))),
pos=4, cex=par()$cex.lab)
}

}

Expand Down Expand Up @@ -150,21 +171,21 @@ all_indicators <- function(X, indicators = c("Variance", "Autocorrelation", "Ske

## mar is inner margins, in order bottom, left, top, right.
## oma is outer margins, default to 0
par(mfrow=c(m+1,n), oma=c(8,8,8,4), mar=c(0,2,0,2),...)
par(mfrow=c(m+1,n), oma=c(3,3,2,.2), mar=c(0,1,0,1), ...)
for(i in 1:n){
plot(X[[i]], type="o", ylab="data", xaxt="n", ...)
mtext(data_names[i], NORTH<-3, cex=par()$cex.lab, line=2) ## data names on each col
if(i==1) mtext("data", WEST<-2, line=4, cex=par()$cex.lab) ## "data" y-axis label
plot(X[[i]], type="l", ylab="data", xaxt="n", ...)
mtext(data_names[i], NORTH<-3, cex=par()$cex.lab, line=1) ## data names on each col
if(i==1) mtext("data", WEST<-2, line=3, cex=par()$cex.lab, las=0) ## "data" y-axis label
}
## Starts on next row, and goes across datasets
for(j in 1:m){
if(j == m){ xaxt <- "s"
} else { xaxt <- "n"
}
for(i in 1:n){
plot_indicator(X[[i]], indicators[j], xaxt=xaxt, method=method, ...)
if(i==1) mtext(indicators[j], WEST<-2, line=4, cex=par()$cex.lab) ## stat name on each row
if(j==m) mtext("time", SOUTH<-1, line=4, cex=par()$cex.lab) ## x-axis label
plot_indicator(X[[i]], indicators[j], xaxt=xaxt, method=method, xpos=-15, ...)
if(i==1) mtext(indicators[j], WEST<-2, line=3, cex=.8*par()$cex.lab, las=0) ## stat name on each row
if(j==m) mtext("time", SOUTH<-1, line=2, cex=.8*par()$cex.lab) ## x-axis label
}
}
}
Expand Down
4 changes: 2 additions & 2 deletions warningsignals/R/montecarlotest.R
View file
Expand Up @@ -239,11 +239,11 @@ plot.pow <- function(pow, main="", legend=FALSE, type="density", test_dist=TRUE,
else if (shade_aic==TRUE & test_dist==TRUE){
legend("topright", c( paste("False Alarms (", aic_wrong*100, "%)", sep=""),
paste("Missed Events (", (1-aic_power)*100, "%)", sep="")),
pch=c(15,15), col=c(rgb(1,.5,1,.5), rgb(1,1,0,.5)))
pch=c(15,15), col=c(rgb(1,.5,0,.5), rgb(1,1,0,.5)))
}
else if (shade_aic==TRUE & test_dist==FALSE){
legend("topright", paste("False Alarms (", aic_wrong*100, "%)", sep=""),
pch=15, col=rgb(1,.5,1,.5))
pch=15, col=rgb(1,.5,0,.5))
}
}
}
Expand Down
6 changes: 4 additions & 2 deletions warningsignals/R/tau_dist_montecarlo.R
View file
Expand Up @@ -83,8 +83,10 @@ plot.tau_dist_montecarlo <- function(out, show_p=TRUE, show_error=FALSE, thresho
threshold_tail <- sort(null_dist)[ round(threshold*nboot) ]
power <- sum(test_dist > threshold_tail)/nboot
p <- 1-sum(null_dist < out$observed[1])/nboot
text(xshift(0), yshift(95), paste("Type I:", round(p,2)), pos=4 )
text(xshift(0), yshift(85), paste("Type II:", round(1-power,2)), pos=4 )
text(xshift(0), yshift(95), paste("Type I:", round(p,2)), pos=4,
cex=par()$cex.lab)
text(xshift(0), yshift(85), paste("Type II:", round(1-power,2)),
pos=4, cex=par()$cex.lab )
}
}

Expand Down
26 changes: 14 additions & 12 deletions warningsignals/demos/paper_plots.R
View file
@@ -1,8 +1,10 @@
equire(warningsignals)
require(warningsignals)

PNG=TRUE
PNG=FALSE

JPEG=TRUE

#load output of figure1.R
load("5554763401.Rdat")
source("../R/indicators.R")
Expand All @@ -11,10 +13,10 @@ indicators <- c("Variance", "Autocorrelation", "Skew", "Kurtosis")
#indicators <- c("Variance", "Autocorrelation")


if(PNG){ png(file="figure1.png", height=length(indicators)*240, width=3*480)
if(JPEG){ jpeg(file="Boettiger_fig1.jpg", height=length(indicators)*89/4, width=89, units="mm", quality=100, res=150)
} else { cairo_pdf(file="figure1.pdf",height=length(indicators)*7/2, width=3*7) }
all_indicators( list(Deteriorating=deteriorating, Constant=constant, Empirical=deut3),
indicators=indicators, cex.axis=2, cex.lab=2.3)
indicators=indicators, cex.axis=.5, cex.lab=.6, lwd=.5)
dev.off()


Expand All @@ -25,10 +27,10 @@ load("5554848679.Rdat")
#deterior_taus <- deterior_taus[1:2]; constant_taus <- constant_taus[1:2]; deut3_taus <- deut3_taus[1:2]


if(PNG) { png(file="figure2.png", width=3*480/3, height=480*length(constant_taus)/3)
if(JPEG) { jpeg(file="Boettiger_fig2.jpg", height=length(indicators)*89/3, width=89, units="mm", quality=100, res=150)
} else { cairo_pdf(file="figure2.pdf", width=3*7/3, height=7*length(constant_taus)/3) }
plot.bootstrap_tau(list(Deteriorating=deterior_taus, Constant=constant_taus, Empirical=deut3_taus),
cex.axis=1, show_p=FALSE, ylim=c(0,2.8))
cex.axis=.6, cex.lab=.8, show_p=FALSE, ylim=c(0,2.8))
dev.off()


Expand All @@ -37,20 +39,20 @@ dev.off()
load("35563325713.Rdat")
data_names <- c("Deteriorating", "Constant", "Empirical")

if(PNG) { png(file="figure3.png", width=480*1.5, height=480/1.5)
if(JPEG) { jpeg(file="Boettiger_fig3.jpg", height=89/2, width=89, units="mm", quality=100, res=150)
} else { cairo_pdf(file="figure3.pdf", width=7*3/3, height=3.5) }

par(mfrow=c(1,3), oma=c(8,8,8,4), mar=c(0,0,0,0))
par(mfrow=c(1,3), oma=c(3,3,3,1), mar=c(0,0,0,0), cex.lab=.8, cex.axis=.8)
plot(deterior_mc,show_text = c("p","power"), xlab="", main="", cex.lab=1, ylim=c(0,.4), xlim=c(0,90))
mtext(data_names[1], NORTH<-3, cex=par()$cex.lab, line=2) ## data labels on top row
mtext("Probability density", WEST<-2, line=4) ## statistic name on first column
mtext(data_names[1], NORTH<-3, cex=par()$cex.lab, line=1) ## data labels on top row
mtext("Probability density", WEST<-2, line=2, cex=par()$cex.lab) ## statistic name on first column

plot(constant_mc,show_text = c("p","power"), xlab="", main="", cex.lab=1, ylim=c(0,.4), yaxt="n", ylab="", xlim=c(0,90))
mtext(data_names[2], NORTH<-3, cex=par()$cex.lab, line=2) ## data labels on top row
mtext("Likelihood Ratio", SOUTH<-1, line=4) ## x-axis label
mtext(data_names[2], NORTH<-3, cex=par()$cex.lab, line=1) ## data labels on top row
mtext("Likelihood Ratio", SOUTH<-1, line=2, cex=par()$cex.lab) ## x-axis label

plot(deut3_mc,show_text = c("p","power"), xlab="", main="", cex.lab=1, ylim=c(0,.4), yaxt="n", ylab="", xlim=c(0,90) )
mtext(data_names[3], NORTH<-3, cex=par()$cex.lab, line=2) ## data labels on top row
mtext(data_names[3], NORTH<-3, cex=par()$cex.lab, line=1) ## data labels on top row

dev.off()

Expand Down

0 comments on commit e1c77f9

Please sign in to comment.