# Import data

rocket<-read.csv(file="c:/RData/rocket.csv", header=T)
rocket

names(rocket)  # list all variable names in rocket data frame
dim(rocket)  # there are 3 variables and 20 observations in rocket data frame

# Enter data and create a data fram
y<-c(2158.70, 1678.15, 2316.00, 2061.00, 2207.50, 1708.30, 
     1784.70, 2575.00, 2357.90, 2256.70, 2165.20, 2399.55, 1779.80, 
     2336.75, 1765.30, 2053.50, 2414.40, 2200.50, 2654.20, 1753.70)
x<-c(15.50, 23.75, 8.00, 17.00, 5.50, 19.00, 24.00, 2.50, 7.50, 
     11.00, 13.00, 3.75, 25.00, 9.75, 22.00, 18.00, 6.00, 12.50, 
     2.00, 21.50)
rocket1<-data.frame(y,x)
rocket1
names(rocket1)<-c("strength", "age")
head(rocket1)  
rocket1<-edit(rocket1) # can enter and/or edit data in an edit windows 

# Scatter plot

plot(rocket$age, rocket$strength) # make a scatter plot

attach(rocket)
par(mfrow=c(2,2))
plot(age, strength)  # points are circles
plot(age, strength, pch=16)  # points are solid circles
plot(age, strength, pch=16, cex=2)  # cex control size
plot(age, strength, pch=16, cex=2, col=2) # col control colour

par(mfrow=c(1,1))
plot(age, strength, pch=16, cex=1, col=2, xlab="age of propellant", 
     ylab="shear strength", main="Scatter plot", cex.main=2)  # add title

# Plot a regression line on scatter plot
rocket.lm<-lm(strength~age, data=rocket)  

# fit a simple linear regression model
plot(strength ~ age, data=rocket, pch=16, main="Scatter plot")
abline(reg=rocket.lm)  # plot the regression line on scatter plot

# Fit a simple linear regression model
rocket.lm<-lm(strength~age, data=rocket)
rocket.lm
summary(rocket.lm) # estimated coefficients, tests
names(rocket.lm)
coef(rocket.lm)
fitted.values(rocket.lm)
resid(rocket.lm)
fitted<-fitted.values(rocket.lm)
residual<-resid(rocket.lm)
cbind(strength, fitted, residual)
round(cbind(strength, fitted, residual),2)

# Compute ANOVA table
rocket.lm<-lm(strength~age, data=rocket)
anova(rocket.lm)

# Compute confidence intervals by formula
n<-length(strength)
n
beta1<-coef(rocket.lm)[2]
beta0<-coef(rocket.lm)[1]
SE.beta1<-2.889
SE.beta0<-44.184
qt(.975, n-2)

# 95% C.I. for beta1
#
c(beta1-qt(.975, n-2)*SE.beta1, beta1+qt(.975, n-2)*SE.beta1) 

# 95% C.I. for beta0
#
c(beta0-qt(.975, n-2)*SE.beta0, beta0+qt(.975, n-2)*SE.beta0) 

# 95% C.I. for sigma^2
#
MSE<-9236
qchisq(.975, n-2)
qchisq(.025, n-2)
c((n-2)*MSE/qchisq(.975, n-2), (n-2)*MSE/qchisq(.025, n-2)) 

# Correlation analysis

delivery<-read.csv(file="c:/RData/Delivery.csv", header=T)
delivery
attach(delivery)

r<-cor(Time, Cases) # compute the correlation coefficeint
r

# Test for pupulation correlation coefficient

n<-length(Time)
t<-r/sqrt((1-r^2)/(n-2))  # test for rao=0
t
qt(.975, n-2)
pvalue<-2*(1-pt(abs(t),n-2))
pvalue