#### R tips
#########################################################
## R tips: how to use R helps?
# want to know how to use function "runif"
?runif
# want to know how to generate random numbers from normal distribution
help.search("normal distribution")
?Normal
#########################################################

#########################################################
## R tips: how to print or save R graphs?
# if use menu operations: [1] click the R graph window; 
#                         [2] (menu)File --> print or File --> Save as
# if use R command lines: 
x <- rnorm(1000, mean=1, sd=2)
postscript("normal.ps",paper="special",width=8,height=7,horizontal=F)
par(mfrow=c(1,1))
hist(x, freq=F)
curve(dnorm(x, mean=1, sd=2), lwd=2, lty=2, add=T)
dev.off()
#########################################################


### Example: Table 8.1-1 on page 440
## read data from a file, if data file is in the working folder
data = read.table("table811.dat", head=T)
x = data$Weight
y = data$GPM

## plot Figure 8.1-2
par(mfrow=c(1,1))
plot(x, y, xlab="Weight (in 1000 lb)", ylab="GPM (gallongs/100 miles)", 
     main="Figure 8.1-2: Scatter plot of gpm against weight")

# print graph into ps file
postscript("figure812.ps",paper="special",width=8,height=7,horizontal=F)
par(mfrow=c(1,1))
plot(x, y, xlab="Weight (in 1000 lb)", ylab="GPM (gallongs/100 miles)", 
     main="Figure 8.1-2: Scatter plot of gpm against weight")
dev.off()

# print graph into eps file
postscript("figure812.eps",paper="special",width=8,height=7,horizontal=T)
par(mfrow=c(1,1))
plot(x, y, xlab="Weight (in 1000 lb)", ylab="GPM (gallongs/100 miles)",  
     xlim=c(0,5), ylim=c(0,8), 
     main="Figure 8.1-2: Scatter plot of gpm against weight")
dev.off()

## estimate coefficients (beta0, beta1)
xbar = mean(x)
ybar = mean(y)
xybar = mean(x*y)
x2bar = mean(x*x)
beta1 = (xybar - xbar*ybar)/(x2bar - xbar^2)
beta0 = ybar - beta1*xbar
beta1
#[1] 1.638996
beta0
#[1] -0.3630888

## use R function "lm" directly
lm(y~x)

# see more output of "lm"
fit1 = lm(y~x)
summary(fit1)

# plot output of "lm"
par(mfrow=c(2,2))
plot(fit1)

## plot residual plot directly
par(mfrow=c(1,1))
plot(fit1$fitted, fit1$residual, xlab="Fitted values", ylab="Residuals", main="Residuals vs. Fitted")
abline(0, 0, lty=2)