Question

2. R programming

2·The data set prostate in the faraway package is froma study on 97 men with prostate cancer who were due to receive a radical prostatectomy We are interest is in predicting lpsa (log prostate specific antigen) with Icavol (log cancer volume). (a) Draw a scatterplot -does a simple linear regression model seem reasonable? (b) Without using the R function Im), compute the values T,Y, Sxx, Syy and Sxy. Com- pute the ordinary least squares estimates of the intercept and slope for the simple linear regression model, and draw the fitted line on your plot from part (a). (c) Obtain the estimate of σ2 and find the estimated standard errors of A, and A. d) Find the estimated covariance between Ao and β1. (e) Carry out t-tests for the two null hypotheses β)-0 and A-0, reporting the value of the test statistic and a p-value in each case. (f) Use the R function Im() to fit the regressions of Ipsa on Icavol

Answer 1

Solution:

R Programming:

a) R code

#install the faraway package if it is not already installed
install.packages('faraway')

library(faraway)
names(prostate)
#a) Draw a scatter plot
plot(prostate$lcavol,prostate$lpsa,xlab="lcavol",ylab="lpsa",main="lpsa vs lcavol")

#get this plot

o O 0 2 4 Icavol

We can see that there is an overall positive linear relationship between lspa and lcavol. The log of prostate specific antigen (lspa) seems to increase with the increase in log cancer vol (lcavol).

A simple linear regression model seems reasonable.

b) The regression line that we want to fit is

$y=\beta_0+\beta_1x+\epsilon$

where y = lspa

$\beta_0$ is the intercept of the regression line

$\beta_1$ is the slope coefficient corresponding to x=lcavol

$\epsilon \stackrel{iid}\sim \mathcal{N}(0,\sigma^2)$ is a random error

We calculate the following

$\begin{align*} \bar{x}&=\frac{\sum x}{n}\\ \bar{y}&=\frac{\sum y}{n}\\ S_x&=\sum(x_i-\bar{x})^2\\ S_y&=\sum(y_i-\bar{y})^2\\ S_{xy}&=\sum(x_i-\bar{x})(y_i-\bar{y})\\ \end{align*}$

and the estimates of slope and intercept using

$\begin{align*} \hat{\beta}_1&=\frac{S_{xy}}{S_x}\\ \hat{\beta}_0&=\bar{y}-\hat{\beta}_1\bar{x} \end{align*}$

The fitted value of y is

$\begin{align*} \hat{y}=\hat{\beta}_0+\hat{\beta}_1x \end{align*}$

The following R code does all these

#part b)
y<-prostate$lpsa
x<-prostate$lcavol
#sample means
xbar<-mean(x)
ybar<-mean(y)
#sum of sqaures
Sx<-sum((x-xbar)^2)
Sy<-sum((y-ybar)^2)
Sxy<-sum((x-xbar)*(y-ybar))
#estimate the value of slope
beta1hat<-Sxy/Sx
#Estimate the value of intercept
beta0hat<-ybar-beta1hat*xbar
sprintf('The estimated value of the intercept is %.4f',beta0hat)
sprintf('The estimated value of the slope is %.4f',beta1hat)
sprintf('The estimated regression line is %.4f+%.4fx',beta0hat,beta1hat)
#calculate the fitted values
yhat<-beta0hat+beta1hat*x
#Draw the fitted line on to the plot from part a)
lines(sort(x),yhat[order(x)],col="red")

# get these outputs

get this plot

c&d) An estimate of $\begin{align*} \sigma^2 \end{align*}$ is

$\begin{align*} SSE=S_y-\hat{\beta}_1S_{xy}\\ \hat{\sigma}^2=MSE=\frac{SSE}{n-2} \end{align*}$

The standard errors of coefficients are

$\begin{align*} s.e(\hat{\beta}_1)&=\sqrt{\frac{\hat{\sigma^2}}{S_x}}\\ s.e(\hat{\beta}_0)&=\sqrt{\frac{\hat{\sigma^2}\sum x^2}{nS_x}}\\ cov(\hat{\beta}_0,\hat{\beta}_1)&=\frac{-\hat{\sigma^2}\bar{x}}{S_x}\\ \end{align*}$

R code

#part c)
#get the number of observations
n<-length(x)
# get the sum of square error
sse<-Sy-beta1hat*Sxy
#get mean square error, which is the estimate of sigma^2
mse<-sse/(n-2)
#estimates of stamdard errors
sb1<-sqrt(mse/Sx)
sb0<-sqrt(mse*sum(x^2)/(n*Sx))
sprintf('The estimated value of sigma^2 %.4f',mse)
sprintf('The standard error of beta1 %.4f',sb1)
sprintf('The standard error of beta0 %.4f',sb0)

#part d)
cov<--mse*xbar/Sx
sprintf('The estimated covariance between beta0&beta 1 %.4f',cov)

#get the following outputs

> sprint f("The estimated value of sigma ^2 %.4f",mse) > sprintf("The standard error of betal %.4f', sbl) [1] "The standard error of betal 0.0682" > sprintf("The standard error of beta0 %.4f', sbO) [1] "The standard error of beta0 0.1219 5 > #part d) > sprintf("The estimated covariance between beta0&beta 1 %.4f', cov) [1] "The estimated covariance between beta0&beta 1 -0.0063"

e) We want to test the following hypotheses for $\begin{align*} \beta_i=0 \end{align*}$ where i=0,1

$\begin{align*} &H_0:\beta_i=0\leftarrow\text{null hypothesis}\\ &H_a:\beta_i\ne 0\leftarrow\text{alternative hypothesis}\\ &\alpha=0.05\leftarrow\text{level of significance to test the hypotheses} \end{align*}$

The test statistics is

$\begin{align*} t=\frac{\hat{\beta}_i-\beta_{iH_0}}{s.e(\hat{\beta}_i)}=\frac{\hat{\beta}_i-0}{s.e(\hat{\beta}_i)}=\frac{\hat{\beta}_i}{s.e(\hat{\beta}_i)} \end{align*}$

this is a 2 tailed test (the alternative hypothesis has "not equal to")

The p-value is

$\begin{align*} \text{p-value}=P(T>t)&plus;P(T<-t) \end{align*}$

the degrees of freedom for t statistics is n-2

Following is the R code

#part e)
#test statistics for beta 0
tb0<-beta0hat/sb0
#p-value of beta0 = P(T>tb0)+P(T<-tb0)
pb0<-pt(abs(tb0),df=n-2,lower.tail=FALSE)+ pt(-abs(tb0),df=n-2,lower.tail=TRUE)
sprintf('The test statistics to test beta0=0 is %.4f, the p-value is %.4f',tb0,pb0)

#test statistics for beta 1
tb1<-beta1hat/sb1
#p-value of beta1 = P(T>tb1)+P(T<-tb1)
pb1<-pt(abs(tb1),df=n-2,lower.tail=FALSE)+ pt(-abs(tb1),df=n-2,lower.tail=TRUE)
sprintf('The test statistics to test beta1=0 is %.4f, the p-value is %.4f',tb1,pb1)

# get these

We will reject the null hypothesis if the p-value is less than the significance level of alpha=0.05

Here for both $\begin{align*} \beta_0,\beta_1 \end{align*}$ the p-values are less than 0.05.

Hence we reject the null hypothesis.

We conclude that there is sufficient evidence to support the claim that the coefficients are significant.

f) Use lm()

R code

#part f) use lm()
m<-lm(lpsa~lcavol,data=prostate)
summary(m)

# get these

we can see that what we have calculated in part a to e), match with this output