nitric.lm1 <- lm(y ~ x1+x2+x3, data=nitric)

drop1(nitric.lm1)

anova(nitric.lm1)

nitric.lm2 <- update(nitric.lm1, .~ . -x3)

plot(nitric.lm2)

anova(nitric.lm2)

nitric.lm3 <- lm(y~x1+x2, data=nitric, subset=1:20)

predict(nitric.lm2, se.fit=T, ci.fit=T)

predict(nitric.lm2, pi.fit=T)

predict(nitric.lm2, pred.new, se.fit=T, ci.fit=T)

predict(nitric.lm2, pi.fit=T)



sd <- function(x) {  
# This calculates the variance of the argument
  
# first and then takes square root of it. 

# The square root is returned

s <- sqrt(var(x))

return(s) 
}




ci <- function(mod,a=0.95) { 

se <- sqrt(diag(summary(mod)$cov.unscaled))

se <- se * summary(mod)$sigma

aa <- (1+a)/2

upper <- coef(mod) + (qt(aa, mod$df.residual) * se)

lower <- coef(mod) - (qt(aa, mod$df.residual) * se)

return(upper, lower)
}


butterfly <- function(color = 8)
{
	theta <- seq(from=0.0, to=24 * pi, len = 2000)
	radius <- exp(cos(theta)) - 2 * cos(4 * theta)
	radius <- radius + sin(theta/12)^5
	x <- radius * sin(theta)
	y <-  - radius * cos(theta)
	plot(x, y, type = "l", axes = F, xlab = "", ylab = "", col = color)
}