cor.drop <-
function(x, y=NULL, cases=NULL, delete.cases=NULL, na.action=c("listwise.deletion","mean.replace"),method=c("pearson","kendall","spearman"), alpha=0.05)
{
# Verifica di na.action
na.action <- tolower(na.action[1])
na.action <- match.arg(na.action)
# Verifica di method
method <- tolower(method[1])
mathod <- match.arg(method)
# Riorganizzazione dati
if(is.vector(x)) {
if(is.null(y))
stop("You must specify the vector \'y\'")
else
df <- data.frame(x,y)
if(is.null(cases))
df$cases <- 1:nrow(df)
else
df$cases <- cases
df$cases <- as.character(df$cases)
} else {
df <- x
if(is.matrix(df))
df <- data.frame(df)
if(ncol(df) == 2)
df$cases <- as.character(1:nrow(df))
colnames(df) <- c("x","y","cases")
}
if(na.action == "listwise.deletion" | na.action == "pairwise.deletion")
df <- na.omit(df)
else {
na.check <- list(x = is.na(df$x), y = is.na(df$y))
if(sum(na.check$x) > 0)
df$x[which(na.check$x)] <- mean(df$x[-which(na.check$x)])
if(sum(na.check$y) > 0)
df$y[which(na.check$y)] <- mean(df$y[-which(na.check$y)])
}
na.check <- is.na(df$cases)
if(sum(na.check) > 0)
df$cases[which(na.check)] <- paste("NA",1:sum(na.check),sep=".")
if(!is.null(delete.cases))
df <- df[-which(df$cases %in% as.character(delete.cases)),]
# Calcolo delle correlazioni
r.test <- cor.test(df[,1],df[,2],method=method)
n <- nrow(df)
rdel <- numeric(n)
for(i in 1:n)
rdel[i] <- cor(df[-i,1],df[-i,2],method=method)
r.diff <- rdel-r.test$estimate
m <- mean(r.diff)
z <- qnorm(1-alpha/2)
se <- sd(r.diff)/sqrt(length(r.diff))
output <- list(
data = df,
overall.test = r.test,
method = method,
correlations = data.frame(cases=df$cases,cor=rdel),
diff.ci = c(lower=m-z*se,upper=m+z*se)+r.test$estimate,
na.action = na.action
)
class(output) <- "cor.drop"
return(output)
}
print.cor.drop <-
function(x,...)
{
sig <- ifelse(x$overall.test$p.value < 0.001, "(p < 0.001)",
ifelse(x$overall.test$p.value < 0.01, "(p < 0.01)",
ifelse(x$overall.test$p.value < 0.05, "(p < 0.05)","(p = n.s.)")))
method <- paste0(
toupper(substr(x$method,1,1)),
substr(x$method,2,nchar(x$method)),"\'s"
)
cat("\nOverall",method,"correlation: r =",round(x$overall.test$estimate,4),sig,"\n\n")
cat("Correlations calculated deleting single cases:\n")
x$correlations$cor <- round(x$correlations$cor,4)
print(x$correlations)
}
plot.cor.drop <-
function(x, vertical=FALSE, main=NULL, axis.lim=NULL, mar.bottom=5.1, mar.left=4.1, mar.top=4.1, mar.right=2.1,...)
{
y <- x$correlations$cor
s <- x$correlations$cases
r <- x$overall.test$estimate
n <- nrow(x$data)
# Definizione dettagli grafici:
cor.lab <- paste0(
toupper(substr(x$method,1,1)),
substr(x$method,2,nchar(x$method)),
"\'s correlation index"
)
mar <- c(mar.bottom,mar.left,mar.top,mar.right)
cex <- 1.7
vertical <- !vertical
# Calcolo dei limiti dell'asse orizzontale
lim <- range(y)
fcround <- function(lim) {
rnd <- round(lim,1)
lim[1] <- ifelse(lim[1]-rnd[1] >=0, rnd[1], rnd[1]-0.1)
lim[2] <- ifelse(lim[2]-rnd[2] <=0, rnd[2], rnd[2]+0.1)
return(lim)
}
if(is.null(axis.lim))
axis.lim <- lim <- fcround(lim=range(y))
else
lim <- fcround(lim=axis.lim)
# Visualizzazione grafico:
par(mar=mar)
if(!vertical) {
stripchart(y ~ s, vertical=vertical, las=2,
pch=21, cex=cex, bg="grey",
xlim=axis.lim, xlab="", ylab="")
rect(x$diff.ci[1],0,x$diff.ci[2],n+1,col="gray90",border=NA)
for(i in 1:n) segments(-1,i,1,i,col="gray70")
for(i in seq(lim[1],lim[2],by=0.1))
segments(i,0,i,n+1,lty=2,col=ifelse(round(i,1)!=0,"gray70","red3"))
segments(r,0,r,n+1,lty=1,lwd=2,col="blue")
par(mar=mar,new=TRUE)
stripchart(y ~ s, vertical=vertical, las=2, pch=21, cex=cex, bg="grey",
main=main, xlim=axis.lim, xlab=cor.lab, ylab="")
} else {
stripchart(y ~ s, vertical=vertical, las=2, pch=21, cex=cex, bg="grey",
ylim=axis.lim, xlab="", ylab="")
rect(0,x$diff.ci[1],n+1,x$diff.ci[2],col="gray90",border=NA)
for(i in 1:n) segments(i,-1,i,1,col="gray70")
for(i in seq(lim[1],lim[2],by=0.1))
segments(0,i,n+1,i,lty=2,col=ifelse(round(i,1)!=0,"gray70","red3"))
segments(0,r,n+1,r,lty=1,lwd=2,col="blue")
par(mar=mar,new=TRUE)
stripchart(y ~ s, vertical=vertical, las=2, pch=21, cex=cex, bg="grey",
main=main, ylim=axis.lim, xlab="", ylab=cor.lab)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.