### Helper function related to linear models
# Calculate number of sigma for each observation, using linear model built with all other observations
CalcLmSigma <- function(x, y, plot=FALSE, xlab=NA, ylab=NA) {
# x, y Vectors to fit linear model lm(y~x). Must have the same length
# plot Whether to make the regression plot
# xlab, ylab Axis labels for the plot
c <- cbind(x, y);
sigma <- sapply(1:nrow(c), function(i) {
x<-c[-i, ];
mdl<-lm(x[, 2]~x[, 1]);
coe<-as.vector(coefficients(mdl));
prd<-c[i, 1]*coe[2]+coe[1];
sig<-summary(mdl)$sigma;
(c[i, 2]-prd)/sig;
});
if (plot) {
par(mar=c(5,5,2,2));
col <- rep('#88888888', length(x));
col[abs(sigma) >= 1.5] <- '#FF888888';
col[abs(sigma) >= 3.0] <- '#FF000088';
plot(x, y, col=col, pch=19, cex=2, xlab=xlab, ylab=ylab, cex.lab=1.5);
points(x, y, col='black', pch=19, cex=.2);
abline(lm(y~x), col='blue', lty=2);
if (cor(c[, 1], c[, 2], use='pair') > 0) loc <- 'topleft' else loc <- 'topright';
legend(loc, legend=c('Num_Sigma < 1.5', 'Num_Sigma >= 1.5', 'Num_Sigma >=3.0'), cex=1, bty='n', pch=19,
col=c('#88888888', '#FF888888', '#FF000088'));
abline(lm(y~x), col='blue', lty=2);
}
names(sigma) <- names(y);
sigma;
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.