Assignment 1

You may discuss homework questions with other students, but whatever you end up submitting has to be your own work. Please combine all of your answers, the computer code and the figures into one PDF file, and submit everything via EMail at aleksey@binghamton.edu. Due date: 11:59 PM October 10, 2025 (Friday evening).

Question 1

This question is from the textbook: Chatterjee and Hadi (CH) page 23, Exercises 1.2. Give an example in any area of interest to you where regression analysis can be used as a data analytic tool to answer some questions of interest.

1. What is the question of interest.
2. Identify the response and the predictor variables.
3. Classify each of the variables as either quantitative or qualitative.
4. Which type of regression (see CH page 18, Table 1.15) can be used to analyze the data?
5. Give a possible form of the model and identify its parameters.

Provide your answers below this horizontal line.

1. Can we predict a pregnant woman’s risk level based on her age, body temperature, blood glucose levels, blood pressure, and heart rate?
2. Predictor variables (X): age, body temperature, blood glucose levels, blood pressure, and heart rate. Response variable (Y): pregnant woman’s risk level
3. Age: Quantitative Body temperature: Quantitative Blood glucose levels: Quantitative Blood pressure: Quantitative Heart rate: Quantitative Pregnant woman’s risk level: Qualitative
4. A logistic regression model could be used to analyze the data.
5. A possible form of the model would be yhat = β0 + Age(β1) + BodyTemp(β2) + GlucoseLevel(β3) + BloodPressure(β4) + HeartRate(β5) + Error Parameters: β0-5 and Error

---

Question 2

On Groundhog Day, February 2, a famous groundhog in Punxsutawney, PA is used to predict whether a winter will be long or not based on whether or not he sees his shadow. Jonathan Taylor collected data on whether Phil saw his shadow or not from here. Some of the data are saved in this table. Although Phil is on the East Coast, Jonathan (based in the bay area) wondered if the information says anything about whether or not he will experience a rainy winter out here in California. For this, Jonathan found rainfall data. The dataset is available for download at here. Here is how this was extracted using R packages.

library(rvest)
url = "http://cdec.water.ca.gov/cgi-progs/precip1/8STATIONHIST"
webpage = readLines(url)[476:593]
ind_of_rows_with_space = which(webpage == webpage[2])
webpage = webpage[-ind_of_rows_with_space]
# remove white space in each line
webpage = lapply(as.list(webpage), function(x){trimws(x)}) %>% unlist()
library(stringr)
temp = lapply(as.list(webpage), function(x){strsplit(x, " ", fixed = TRUE)}) %>% unlist()
temp = temp[-which(temp == temp[2])]
df = matrix(temp, ncol = 14, byrow = TRUE)

col.names.df = df[1, ]
df = df[-1, ] %>% data.frame()
colnames(df) = col.names.df
df = apply(df, 2, function(x){x = as.numeric(x)}) %>% data.frame()

write.csv(df, file = "rainfall.csv", row.names = FALSE)

1. Using the ggplot2 package, make a box plot of the mean monthly rainfall (total annual rainfall divided by 12 months) rainfall in Northern California comparing the years Phil sees his shadow versus the years he does not. [Read Phil’s data from this link. Download and read rainfall data from here. [Use dplyr::filter to select years of rainfall data that Phil sees his shadow or not. Use dplyr::left_join to merge Phil’s shadow data and rainfall data.]

Below, I am providing my answer for you. Study what each line does

library(dplyr)

url_groundhog = 'http://people.math.binghamton.edu/qiao/data501/data/groundhog.table'
url_rainfall = 'http://people.math.binghamton.edu/qiao/data501/data/rainfall.csv'
groundhog <- read.csv(url(url_groundhog))
rainfall <- read.csv(url(url_rainfall))
names(rainfall)[1] <- "year" # change the name of the first column so that left_join below
# knows how to merge the two tables.
merged <- left_join(groundhog,rainfall)

print(head(merged))

merged$mean_rainfall = merged$Total / 12
## group_by, then summarise. This is not needed by this question. Just a showcase of how it is used.
## Note that this is similar to the groupby method in pandas
merged %>% group_by(shadow) %>% summarise(mean(mean_rainfall))

library(ggplot2)
ggplot(data = merged, aes(y=mean_rainfall, x =shadow, color=shadow)) + geom_boxplot()

2. Construct a 95% confidence interval for the difference between the mean monthly rainfall (total annual rainfall divided by 12 months) in years Phil sees his shadow and years he does not. This is a test comparing the means of two population. Feel free to use existing built-in function in R for this part, or use formula provided in the class. In any case, assume that the two groups have the same variance.

t.test(mean_rainfall ~ shadow, data = merged, var.equal = TRUE)

    Two Sample t-test

data:  mean_rainfall by shadow
t = 0.54731, df = 19, p-value = 0.5905
alternative hypothesis: true difference in means between group N and group Y is not equal to 0
95 percent confidence interval:
 -1.139317  1.946150
sample estimates:
mean in group N mean in group Y 
       4.701333        4.297917 

3. Interpret the procedure used to construct in part 2. What do we really know about confidence intervals? What do confidence intervals promise to achieve?

I used a two sample t-test to test the difference in rain fall between the years that Phil sees his shadow and the ones that he doesn’t see his shadow. When Phil did not see his shadow, there was 4.70 inches of rainfall per month. When he did see his shadow, there was 4.30 inches of rainfall per month. The confidence interval is (-1.14, 1.95).

Confidence intervals provide us with a range of values likely to contain a population parameter that is unknown. For this problem, the confidence intervals of 95% shows us a range of values likely to contain the difference in mean rainfall between shadow and no-shadow years. The 95% confidence interval means that if we were to repeat this study over and over again, 95% of the resulting intervals would contain the true parameter value.

4. At level, α = 0.05 would you reject the null hypothesis that the average rainfall in Northern California during the month of February was the same in years Phil sees his shadow versus years he does not?

I would fail to reject the null hypothesis at α = 0.05 because the p-value is greater than 0.05, so it is not statistically significant.

5. What assumptions are you making in forming your confidence interval and in your hypothesis test?

6. The data should be collected using a random sampling method.

7. Each observation in the data set should be independent of every other observation.

8. The two groups have the same population variance of monthly rainfall.

---

Question 3

The data set walleye in the package alr4 (remember you may have to run install.packages(“alr4”)) of data measured on walleye fish in Wisconsin.

1. Create a box plot (use ggplot2) of length, for age in 5:8 separately. Use filter in dplyr package to filter rows from the data frame.

library(dplyr)
library(magrittr)
library(ggplot2)
library(alr4)

head(walleye)

str(walleye)

'data.frame':   3198 obs. of  3 variables:
 $ age   : int  1 1 1 1 1 1 1 1 1 1 ...
 $ length: num  215 193 203 201 232 ...
 $ period: int  1 1 1 1 1 1 1 1 1 1 ...

walleye_sub = walleye %>% filter(age >= 5 & age <= 8)
walleye_sub$age = factor(walleye_sub$age)
ggplot(data = walleye_sub, aes(y = length , x = age)) + geom_boxplot(fill = "orange")

2. Compute the sample mean and the sample standard deviation length in the four age groups (5:8) respectively (use dplyr). Your answer should make use of the group_by and summarise functions and take up only one line.

walleye_sub %>% group_by(age) %>% summarise(mean_length = mean(length, na.rm = TRUE), sd_length = sd(length, na.rm = TRUE))

3. Create a histogram (use ggplot2) of length within age of 5:8 putting the plots in a 2x2 grid in one file.

ggplot(walleye_sub, aes(x = length)) + geom_histogram(binwidth=3, color="white") + facet_wrap(~ age)

4. Compute a 95% confidence interval for the difference in the mean length in years 5 and 7. What assumptions are you making?

age5 = walleye %>% filter(age == 5) %>% .$length
age7 = walleye %>% filter(age == 7) %>% .$length
t.test(age5, age7, alternative = "two.sided", conf.level = 0.95)$conf.int

[1] -44.37574 -33.28132
attr(,"conf.level")
[1] 0.95

The confidence interval is (-44.38, -33.28). We assume independent observations, approximately normal length distributions within each age group, and no equal-variance requirement.

5. At level α = 10%, test the null hypothesis that the average length in the group age==5 is the same as the in the group age==7. What assumptions are you making? What can you conclude?

t.test(age5, age7)$p.value

[1] 3.896149e-34

The p-value is 3.896149e-34, which is less than 0.05. So we will reject the null hypothesis at α = 10%. We assume that there are independent observations in each group and between groups and approximate normality of lengths within each age. We can conclude that there is average length of age 5 is significantly different from average length of age 7.

6. Repeat the test in 5. using the function lm. Hint: use the subset with age 5 and age 7 only, fit a regression model of length on age (in which age is a factor/categorical variable). Then check the significance of the slope.

age57 = walleye %>% filter(age == c("5", "7")) %>% .[,c("age","length")]
age57$age = factor(age57$age)
summary(lm(length ~ age - 1, data = age57))

Call:
lm(formula = length ~ age - 1, data = age57)

Residuals:
    Min      1Q  Median      3Q     Max 
-60.776 -18.275  -2.989  17.538  78.658 

Coefficients:
     Estimate Std. Error t value Pr(>|t|)    
age5  364.762      2.203   165.6   <2e-16 ***
age7  401.020      3.255   123.2   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 28.38 on 240 degrees of freedom
Multiple R-squared:  0.9944,    Adjusted R-squared:  0.9944 
F-statistic: 2.13e+04 on 2 and 240 DF,  p-value: < 2.2e-16

The p-value < 2e-16 for the slope of the above regression suggesting that we reject the null hypothesis at α = 10% and conclude that the mean length differs between ages 5 and 7.

---

Question 4

The following questions are based on the anscombe data in R.

1. Plot the 4 data sets (x1, y1), (x2, y2), (x3, y3), (x4, y4) on a 2-by-2 grid of plots using ggplot2 and gridExtra package. The results should look like Figure 2.3 in the textbook on page 30. Add the two variable names to each plot as the main title on each plot. Add the axis-labels using bquote to each plot. For example,

library(ggplot2)
library(magrittr)
library(dplyr)
head(anscombe)

p1 <- ggplot(data = anscombe) +
geom_point(aes(x = x1, y = y1)) +
xlab(bquote(x[1])) +
ylab(bquote(y[1])) +
ggtitle(bquote(x[1]~"versus"~y[1])) +
theme(plot.title = element_text(hjust = 0.5))

p1

Hint: use the gridExtra package to display multiple plots in the same figure.

p1 <- ggplot(data = anscombe) + 
  geom_point(aes(x = x1, y = y1)) +
  xlab(bquote(x[1])) +
  ylab(bquote(y[1])) +
  ggtitle(bquote(x[1]~"versus"~y[1])) +
  theme(plot.title = element_text(hjust = 0.5))

p2 <- ggplot(data = anscombe) + 
  geom_point(aes(x = x2, y = y2)) +
  xlab(bquote(x[2])) +
  ylab(bquote(y[2])) +
  ggtitle(bquote(x[2]~"versus"~y[2])) +
  theme(plot.title = element_text(hjust = 0.5))

p3 <- ggplot(data = anscombe) + 
  geom_point(aes(x = x3, y = y3)) +
  xlab(bquote(x[3])) +
  ylab(bquote(y[3])) +
  ggtitle(bquote(x[3]~"versus"~y[3])) +
  theme(plot.title = element_text(hjust = 0.5))

p4 <- ggplot(data = anscombe) + 
  geom_point(aes(x = x4, y = y4)) +
  xlab(bquote(x[4])) +
  ylab(bquote(y[4])) +
  ggtitle(bquote(x[4]~"versus"~y[4])) +
  theme(plot.title = element_text(hjust = 0.5))


library(gridExtra)

grid.arrange(p1, p2, p3, p4, nrow = 2)

2. Fit a regression model to the data sets:

1. y1 ∼ x1

fit1 = lm(y1 ~ x1, data = anscombe)
summary(fit1)

Call:
lm(formula = y1 ~ x1, data = anscombe)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.92127 -0.45577 -0.04136  0.70941  1.83882 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   3.0001     1.1247   2.667  0.02573 * 
x1            0.5001     0.1179   4.241  0.00217 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.237 on 9 degrees of freedom
Multiple R-squared:  0.6665,    Adjusted R-squared:  0.6295 
F-statistic: 17.99 on 1 and 9 DF,  p-value: 0.00217

2. y2 ∼ x2

fit2 = lm(y2 ~ x2, data = anscombe)
summary(fit2)

Call:
lm(formula = y2 ~ x2, data = anscombe)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.9009 -0.7609  0.1291  0.9491  1.2691 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)    3.001      1.125   2.667  0.02576 * 
x2             0.500      0.118   4.239  0.00218 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.237 on 9 degrees of freedom
Multiple R-squared:  0.6662,    Adjusted R-squared:  0.6292 
F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002179

3. y3 ∼ x3

fit3 = lm(y3 ~ x3, data = anscombe)
summary(fit3)

Call:
lm(formula = y3 ~ x3, data = anscombe)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1586 -0.6146 -0.2303  0.1540  3.2411 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   3.0025     1.1245   2.670  0.02562 * 
x3            0.4997     0.1179   4.239  0.00218 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.236 on 9 degrees of freedom
Multiple R-squared:  0.6663,    Adjusted R-squared:  0.6292 
F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002176

4. y4 ∼ x4

fit4 = lm(y4 ~ x4, data = anscombe)
summary(fit4)

Call:
lm(formula = y4 ~ x4, data = anscombe)

Residuals:
   Min     1Q Median     3Q    Max 
-1.751 -0.831  0.000  0.809  1.839 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   3.0017     1.1239   2.671  0.02559 * 
x4            0.4999     0.1178   4.243  0.00216 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.236 on 9 degrees of freedom
Multiple R-squared:  0.6667,    Adjusted R-squared:  0.6297 
F-statistic:    18 on 1 and 9 DF,  p-value: 0.002165

using the command lm. Verify that all the fitted models have the exact same coefficients (up to numerical tolerance).

3. Using the command cor, compute the sample correlation for each data set.

attach(anscombe)
cor(x1,y1)

[1] 0.8164205

cor(x2,y2)

[1] 0.8162365

cor(x3,y3)

[1] 0.8162867

cor(x4,y4)

[1] 0.8165214

4. Fit the same models in 2. but with the x and y reversed. Using the command summary, does anything about the results stay the same when you reverse x and y?

When you reverse x and y, t-value and p-value for the predictor, the R^2 and adjusted R^2, the F-statistic, and the significance conclusion stays the same. However, the coefficients differ drastically, along with the residual standard error and standard error of the coefficients.

1. x1 ∼ y1

reversefit1 = lm(x1 ~ y1, data = anscombe)
summary(reversefit1)

Call:
lm(formula = x1 ~ y1, data = anscombe)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.6522 -1.5117 -0.2657  1.2341  3.8946 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -0.9975     2.4344  -0.410  0.69156   
y1            1.3328     0.3142   4.241  0.00217 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.019 on 9 degrees of freedom
Multiple R-squared:  0.6665,    Adjusted R-squared:  0.6295 
F-statistic: 17.99 on 1 and 9 DF,  p-value: 0.00217

2. x2 ∼ y2

reversefit2 = lm(x2 ~ y2, data = anscombe)
summary(reversefit2)

Call:
lm(formula = x2 ~ y2, data = anscombe)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.8516 -1.4315 -0.3440  0.8467  4.2017 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -0.9948     2.4354  -0.408  0.69246   
y2            1.3325     0.3144   4.239  0.00218 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.02 on 9 degrees of freedom
Multiple R-squared:  0.6662,    Adjusted R-squared:  0.6292 
F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002179

3. x3 ∼ y3

reversefit3 = lm(x3 ~ y3, data = anscombe)
summary(reversefit3)

Call:
lm(formula = x3 ~ y3, data = anscombe)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.9869 -1.3733 -0.0266  1.3200  3.2133 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -1.0003     2.4362  -0.411  0.69097   
y3            1.3334     0.3145   4.239  0.00218 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.019 on 9 degrees of freedom
Multiple R-squared:  0.6663,    Adjusted R-squared:  0.6292 
F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002176

4. x4 ∼ y4

reversefit4 = lm(x4 ~ y4, data = anscombe)
summary(reversefit4)

Call:
lm(formula = x4 ~ y4, data = anscombe)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.7859 -1.4122 -0.1853  1.4551  3.3329 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -1.0036     2.4349  -0.412  0.68985   
y4            1.3337     0.3143   4.243  0.00216 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.018 on 9 degrees of freedom
Multiple R-squared:  0.6667,    Adjusted R-squared:  0.6297 
F-statistic:    18 on 1 and 9 DF,  p-value: 0.002165

5. Compute the SSE, SST and R2 value for each data set. Use the commands mean, sum, predict and / or resid. (Use the original models, i.e. yi ∼ xi so only 4 SSE values)

lm1 = lm(y1 ~ x1, data = anscombe)
SSE = sum((resid(lm1))^2)
SST = sum((anscombe$y1 - mean(anscombe$y1))^2)
SSR = SST - SSE
R2 = 1 - SSE / SST
c(SSE,SSR,SST,R2)

[1] 13.7626900 27.5100009 41.2726909  0.6665425

lm2 = lm(y2 ~ x2, data = anscombe)
SSE = sum((resid(lm2))^2)
SST = sum((anscombe$y2 - mean(anscombe$y2))^2)
SSR = SST - SSE
R2 = 1 - SSE / SST
c(SSE,SSR,SST,R2)

[1] 13.776291 27.500000 41.276291  0.666242

lm3 = lm(y3 ~ x3, data = anscombe)
SSE = sum((resid(lm3))^2)
SST = sum((anscombe$y3 - mean(anscombe$y3))^2)
SSR = SST - SSE
R2 = 1 - SSE / SST
c(SSE,SSR,SST,R2)

[1] 13.756192 27.470008 41.226200  0.666324

lm4 = lm(y4 ~ x4, data = anscombe)
SSE = sum((resid(lm4))^2)
SST = sum((anscombe$y4 - mean(anscombe$y4))^2)
SSR = SST - SSE
R2 = 1 - SSE / SST
c(SSE,SSR,SST,R2)

[1] 13.7424900 27.4900009 41.2324909  0.6667073

6. Using the ggplot2 package, re-plot the data, adding the regression line to each plot. (Use the original models, i.e. yi ∼ xi so only 4 plots)

p1 <- ggplot(data = anscombe) + 
  geom_point(aes(x = x1, y = y1)) +
  xlab(bquote(x[1])) +
  ylab(bquote(y[1])) +
  ggtitle(bquote(x[1]~"versus"~y[1])) +
  theme(plot.title = element_text(hjust = 0.5)) +
  geom_abline(intercept = coef(lm1)[["(Intercept)"]], slope = coef(lm1)[["x1"]])

p2 <- ggplot(data = anscombe) + 
  geom_point(aes(x = x2, y = y2)) +
  xlab(bquote(x[2])) +
  ylab(bquote(y[2])) +
  ggtitle(bquote(x[2]~"versus"~y[2])) +
  theme(plot.title = element_text(hjust = 0.5)) +
  geom_abline(intercept = coef(lm2)[["(Intercept)"]], slope = coef(lm2)[["x2"]])

p3 <- ggplot(data = anscombe) + 
  geom_point(aes(x = x3, y = y3)) +
  xlab(bquote(x[3])) +
  ylab(bquote(y[3])) +
  ggtitle(bquote(x[3]~"versus"~y[3])) +
  theme(plot.title = element_text(hjust = 0.5)) +
  geom_abline(intercept = coef(lm3)[["(Intercept)"]], slope = coef(lm3)[["x3"]])

p4 <- ggplot(data = anscombe) + 
  geom_point(aes(x = x4, y = y4)) +
  xlab(bquote(x[4])) +
  ylab(bquote(y[4])) +
  ggtitle(bquote(x[4]~"versus"~y[4])) +
  theme(plot.title = element_text(hjust = 0.5)) +
  geom_abline(intercept = coef(lm4)[["(Intercept)"]], slope = coef(lm4)[["x4"]])


library(gridExtra)

grid.arrange(p1, p2, p3, p4, nrow = 2)

---

Question 5

In a recent, exciting, but also controversial Science article, Tomasetti and Vogelstein attempt to explain why cancer incidence varies drastically across tissues (e.g. why one is much more likely to develop lung cancer rather than pelvic bone cancer). The authors show that a higher average lifetime risk for a cancer in a given tissue correlates with the rate of replication of stem cells in that tissue. The main inferential tool for their statistical analysis was a simple linear regression, which we will replicate here.

You can download the dataset as follows:

tomasetti = read.csv("http://people.math.binghamton.edu/qiao/data501/data/Tomasetti.csv")
head(tomasetti)

The dataset contains information about 31 tumor types. The Lscd (Lifetime stem cell divisions) column refers to the total number of stem cell divisions during the average lifetime, while Risk refers to the lifetime risk for cancer of that tissue type.

1. Fit a simple linear regression model to the data with log(Risk) as the response variable and log(Lscd) as the predictor variable.

tomasettifit = lm(log(Risk) ~ log(Lscd), data = tomasetti)
summary(tomasettifit)

Call:
lm(formula = log(Risk) ~ log(Lscd), data = tomasetti)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8013 -1.0721  0.1434  0.9945  2.7873 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -17.52379    1.66436  -10.53 2.02e-11 ***
log(Lscd)     0.53260    0.07316    7.28 5.12e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.725 on 29 degrees of freedom
Multiple R-squared:  0.6463,    Adjusted R-squared:  0.6341 
F-statistic:    53 on 1 and 29 DF,  p-value: 5.117e-08

2. Plot the estimated regression line and the data.

tomasettiplot <- ggplot(data = tomasetti) + 
  geom_point(aes(x = log(Lscd), y = log(Risk))) +
  xlab(bquote(log(Lscd))) +
  ylab(bquote(log(Risk))) +
  ggtitle(bquote(log(Lscd)~"versus"~log(Risk))) +
  theme(plot.title = element_text(hjust = 0.5)) +
  geom_abline(intercept = coef(tomasettifit)[["(Intercept)"]], slope = coef(tomasettifit)[["log(Lscd)"]])

tomasettiplot

3. Add upper and lower 95% prediction bands to predict the response given a range of covariates on the plot, using predict(). That is, produce one line for the upper limit of each interval over a sequence of densities, and one line for the lower limits of the intervals.

newdata = data.frame(Lscd = seq(min(tomasetti$Lscd), max(tomasetti$Lscd), length.out = 200))
newdata[,2:4] = predict(tomasettifit, newdata, interval = "prediction")
colnames(newdata)[2:4] = c("fit","pred.lwr","pred.upr")

library(ggplot2)

predictplot = ggplot(data = tomasetti) + 
  geom_point(mapping = aes(x = log(Lscd), y = log(Risk)), size = 2) +
  geom_line(data = newdata, aes(x = log(Lscd), y = pred.lwr), color = "red", linetype = "dotted", linewidth = 1) +
  geom_line(data = newdata, aes(x = log(Lscd), y = pred.upr), color = "red", linetype = "dotted", linewidth = 1) +
  theme_bw()

predictplot

4. Interpret the above bands at a Lscd = 10^10.

predict(tomasettifit, newdata = data.frame(Lscd = 1e10), interval="predict")

        fit       lwr       upr
1 -5.260253 -8.845771 -1.674736

5. Add upper and lower 95% confidence bands for the conditional mean response on the plot, using predict(). That is, produce one line for the upper limit of each interval over a sequence of densities, and one line for the lower limits of the intervals.

newdata[,5:6] = predict(tomasettifit, newdata, interval = "confidence")[,2:3]
colnames(newdata)[5:6] = c("conf.lwr", "conf.upr")

library(ggplot2)

plot = ggplot(data = tomasetti) + 
  geom_point(mapping = aes(x = log(Lscd), y = log(Risk)), size = 2) +
  geom_line(data = newdata, aes(x = log(Lscd), y = conf.lwr), color = "blue", linetype = "dashed", linewidth = .75) +
  geom_line(data = newdata, aes(x = log(Lscd), y = conf.upr), color = "blue", linetype = "dashed", linewidth = .75) +
  theme_bw() +
  geom_line(data = newdata, aes(x = log(Lscd), y = pred.lwr), color = "red", linetype = "dotted", linewidth = .75) +
  geom_line(data = newdata, aes(x = log(Lscd), y = pred.upr), color = "red", linetype = "dotted", linewidth = .75)

plot

6. Interpret the above bands at a Lscd = 10^10.

predict(tomasettifit, newdata = data.frame(Lscd = 1e10), interval="confidence")

        fit       lwr       upr
1 -5.260253 -5.901806 -4.618701

7. Test whether the slope in this regression is equal to 0 at level α = 0.05. State the null hypothesis, the alternative hypothesis, the conclusion, and the p-value.

summary(tomasettifit)

Call:
lm(formula = log(Risk) ~ log(Lscd), data = tomasetti)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8013 -1.0721  0.1434  0.9945  2.7873 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -17.52379    1.66436  -10.53 2.02e-11 ***
log(Lscd)     0.53260    0.07316    7.28 5.12e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.725 on 29 degrees of freedom
Multiple R-squared:  0.6463,    Adjusted R-squared:  0.6341 
F-statistic:    53 on 1 and 29 DF,  p-value: 5.117e-08

H0: β1 = 0 H1: β1 ≠ 0

The p-value is 5.12e-08, which is less than .05, so we would reject the null hypothesis. Therefore, there is strong evidence that the slope is not zero.

8. What are assumptions you made for part (7) above.

I assumed linearity, that observations are independent, and that error variance is constant across values of log(Lscd).

9. Give a 95% confidence interval for the slope of the regression line.

confint(tomasettifit, level = 0.95)

                  2.5 %      97.5 %
(Intercept) -20.9277953 -14.1197921
log(Lscd)     0.3829708   0.6822268

10. Interpret your interval in (9). What does this interval promise to achieve?

Since the interval does not include 0, the relationship is positive and statistically significant.The interval promises that if we repeat the study over and over again, about 95% of those intervals would contain the true slope β1.

11. Report the R2 of the model.

summary(tomasettifit)$r.squared

[1] 0.6463325

12. Report the adjusted R2 of the model.

summary(tomasettifit)$adj.r.squared

[1] 0.6341371

13. Report an estimate of the variance of the errors in the model.

summary(tomasettifit)$sigma^2

[1] 2.975007

14. Provide an interpretation of the R2 you calculated above, ideally to your neighbor who does not know much about statistics.

About 65% of the differences in log(lifetime cancer risk) across these 31 tissues can be explained by differences in log(lifetime stem-cell divisions). The remaining 35% of the variation is likely due to other factors and noise not in the model.

15. According to a Reuters article “Plain old bad luck plays a major role in determining who gets cancer and who does not, according to researchers who found that two-thirds of cancer incidence of various types can be blamed on random mutations and not heredity or risky habits like smoking.” Is this a correct interpretation of R2?

No it’s no the correct interpretation of R2.

---

Question 6

This question is from our textbook CH Exercises 2.1, Page 53.

In order to investigate the feasibility of starting a Sunday edition for a large metropolitan newspaper, information was obtained from a sample of 34 newspapers concerning their daily and Sunday circulations (in thousands) (Source: Gale Directory of Publications, 1994). The data can be read from the book’s Website: http://www1.aucegypt.edu/faculty/hadi/RABE5/Data5/P054.txt.

1. Read the data using read.table (separator of column is tab and the data frame has variable names).

newspaper = read.table("https://www1.aucegypt.edu/faculty/hadi/RABE5/Data5/P054.txt", sep = "\t", header = TRUE)
head(newspaper)

2. Construct a scatter plot of Sunday circulation versus daily circulation.

library(ggplot2)

newsplot <- ggplot(data = newspaper) + 
  geom_point(aes(x = Daily, y = Sunday)) +
  xlab("Daily Circulation") +
  ylab("Sunday Circulation") +
  ggtitle(bquote("Daily Circulation versus Sunday Circulation")) +
  theme(plot.title = element_text(hjust = 0.5))

newsplot

3. Does the plot suggest a linear relationship between daily and Sunday circulation?

Yes, it does suggest a linear relationship between the two variables.

4. Fit a regression line predicting Sunday circulation from daily circulation (Use lm()).

newsfit = lm(Sunday ~ Daily, data = newspaper)

new_newsplot <- ggplot(data = newspaper) + 
  geom_point(aes(x = Daily, y = Sunday)) +
  xlab("Daily Circulation") +
  ylab("Sunday Circulation") +
  ggtitle(bquote("Daily Circulation versus Sunday Circulation")) +
  theme(plot.title = element_text(hjust = 0.5)) +
  geom_abline(intercept = coef(newsfit)[["(Intercept)"]], slope = coef(newsfit)[["Daily"]])

new_newsplot

5. Is there a significant relationship between Sunday circulation and daily circulation? Justify your answer by a statistical test (Use F test in anova()).

There is a significant relationship between Sunday circulation and daily circulation because the F critical value is 4.15 and the F computed value is 358.53, and F computed value is greater than F critical value. Also, the p-value is less than .05 which shows a significant relationship.

anova(newsfit)

Analysis of Variance Table

Response: Sunday
          Df  Sum Sq Mean Sq F value    Pr(>F)    
Daily      1 4292653 4292653  358.53 < 2.2e-16 ***
Residuals 32  383136   11973                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

qf(0.95, 1, 32)

[1] 4.149097

6. Indicate what hypothesis your are testing and your conclusion for the test in part (5).

H0: There is no relationship between Sunday circulation and Daily circulation. (β1 = 0) H1: There is a relationship between Sunday circulation and Daily circulation. (β1 ≠ 0)

We reject the null hypothesis and accept the alternative.

7. Using the anova table produced in part (5), compute the proportion of the variability in Sunday circulation is accounted for by daily circulation.

SSR <- anova(newsfit)["Daily", "Sum Sq"]
SSE <- anova(newsfit)["Residuals", "Sum Sq"]
R2  <- SSR / (SSR + SSE)
R2

[1] 0.9180597