#########################################################
# BodyBrain Example -- R source code file body-brain.r #
#########################################################

# Read and print data frame.
bodybrain = read.table("body-brain.txt", header=T)
cat("bodybrain data frame:\n")
print(bodybrain)

# Create and print data frame for log-log transform.
loglog = data.frame(logbody=log(bodybrain$BodyWeight),
                    logbrain=log(bodybrain$BrainWeight))
cat("Transformed loglog dataframe:\n")
print(loglog)

# Obtain regression Model 1.
model1 = lm(BrainWeight ~ BodyWeight, data=bodybrain)
cat("Regression Analysis for Model1:\n")
print(summary(model1))

# Obtain regression Model 2.
model2 = lm(logbrain ~ logbody, data=loglog)
cat("Regression Analysis for Model2:\n")
print(summary(model2))

# Create graphs for Model 1.
pdf("body-brain.pdf")
attach(bodybrain)
p = fitted(model1)
r = residuals(model1)
plot(BodyWeight, BrainWeight, 
   main="Scatterplot of BrainWeight vs. BodyWeight",
   xlab="Body Weight (kilos)",
   ylab="Brain Weight (grams)")
abline(model1)
plot(p, r, main="Model 1: Residual Plot",
	       xlab="Predicted Values",
           ylab="Residual Values")
abline(h=0, lty="dashed")
qqnorm(r, main="Model 1: Normal Plot of Residuals",
          ylab="Residuals")

# Create graphs for Model 2.
attach(loglog)
p = fitted(model2)
r = residuals(model2)
plot(logbody, logbrain, 
	 main="Scatterplot of Log(BrainWeight) vs. Log(BodyWeight)",
   xlab="Log of Body Weight",
   ylab="Log of Brain Weight")
abline(model2)
plot(p, r, main="Model 2: Residual Plot",
           xlab="Predicted Values",
           ylab="Residual Values")
abline(h=0, lty="dashed")
qqnorm(r, main="Model 2: Normal Plot of Residuals",
          ylab="Residuals")

dev.off( )