Question

2. The data set prostate in the faraway package is from a 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 lcavol (log cancer volume). (a) Draw a scatterplot - does a simple linear regression model seem reasonable? (b) Without using the R function Im(0, compute the values , 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 AO and A (d) Find the estimated covariance between Bo and B (e) Carry out t-tests for the two null hypotheses Bo 0 and B0, reporting the value of the test statistic and a p-value in each case. (f) Use the R function lm() to fit the regressions of Ipsa on lcavol.

Answer 1

R code with explanations (all statements starting with # are comments)

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

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

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