# Pendul Example -- source code file pendul.r

# Read raw dataset.
raw = read.table("pendul.txt", header=T)
print(raw)

# Create transformed variables.
len = raw$Length * 0.0254
per = raw$Time / 30.0
sqrtlen = sqrt(len)
lenlen = len * len
loglen = log(len)
logper = log(per)
pendul = data.frame(len=len,
	                per=per,
	                lenlen=lenlen,
	                sqrtlen=sqrtlen,
	                loglen=loglen,
                    logper=logper)
print(pendul)	                  
	           
pdf("pendul.pdf")

# Plot Period vs. Length
plot(len, per, main="Original Data",
	         xlab="Length (Meters)",
               ylab="Period (Seconds)")

# Model 1
ms = "Model 1: Simple Linear Regression"
cat("####################################\n")
cat(ms,  "#\n")
cat("####################################\n")
model1 = lm(per ~ len, data=pendul)
print(summary(model1))
r = residuals(model1)
p = fitted(model1)

plot(p, r, main=ms,
	       xlab="Predicted Values",
	       ylab="Residuals")
abline(h=0, lty="dashed")
qqnorm(r, main=ms,
	      ylab="Residuals")

# Model 2
ms = "Model 2: Quadratic Regression"
cat("####################################\n")
cat(ms,      "#\n")
cat("####################################\n")
model2 = lm(per ~ len +lenlen, data=pendul)
print(summary(model2))
r = residuals(model2)
p = fitted(model2)

plot(p, r, main=ms,
	       xlab="Predicted Values",
	       ylab="Residuals")
abline(h=0, lty="dashed")
qqnorm(r, main=ms,
	      ylab="Residuals")

# Model 3
ms = "Model 3: Using Sqrt Transform"
cat("####################################\n")
cat(ms, "#\n")
cat("####################################\n")
model3 = lm(per ~ sqrtlen, data=pendul)
print(summary(model3))
r = residuals(model3)
p = fitted(model3)

plot(p, r, main=ms,
	       xlab="Predicted Values",
	       ylab="Residuals")
abline(h=0, lty="dashed")
qqnorm(r, main=ms,
	      ylab="Residuals")

# Model 4
ms = "Model 4: Using Log-log transform"
cat("####################################\n")
cat(ms,   "#\n")
cat("####################################\n")
model4 = lm(logper ~ loglen, data=pendul)
print(summary(model4))
r = residuals(model4)
p = fitted(model4)

plot(p, r, main=ms,
	       xlab="Predicted Values",
	       ylab="Residuals")
abline(h=0, lty="dashed")
qqnorm(r, main=ms,
	       ylab="Residuals")	 
	         
dev.off( )