climate <- read.csv("https://www.isid.ac.in/~deepayan/RT2018/notes/data/annual.csv")

climate <-
    transform(climate,
              VolcanoF = cut(Volcano, breaks = quantile(Volcano,
                                                        probs = c(0, 1/3, 2/3, 1)),
                             labels=c("low", "medium", "high"),
                             include.lowest = TRUE))



## Fit the model Temp ~ CO2. Report the p-value for CO2 and the value of R2 for the model.

s3 <- summary(lm(Temp ~ CO2, climate))
round(coefficients(s3)["CO2", "Pr(>|t|)"], 4)
round(s3$r.squared, 4)

## Fit the model Temp ~ CH4. Report the p-value for CH4 and the value of R2 for the model.

s4 <- summary(lm(Temp ~ CH4, climate))
round(coefficients(s4)["CH4", "Pr(>|t|)"], 4)
round(s4$r.squared, 4)

## Fit the model Temp ~ VolcanoF. Is VolcanoF significant at level
## 0.05? Report the relevant p-value.

s5 <- anova(lm(Temp ~ VolcanoF, climate))
round(s5["VolcanoF", "Pr(>F)"], 4)

## Fit a model with Temp as response and main effects for CO2 and
## VolcanoF as predictors.

fm6 <- lm(Temp ~ CO2 + VolcanoF, climate)
summary(fm6)
anova(fm6)

## - Is VolcanoF significant in the presence of CO2 at level 0.05?
## - Report the relevant p-value.

round(anova(fm6)["VolcanoF", "Pr(>F)"], 4)

## - Is CO2 significant in the presence of VolcanoF at level 0.05?
## - Report the relevant p-value.

round(coefficients(summary(fm6))["CO2", "Pr(>|t|)"], 4)

## Fit a model with Temp as response and CO2, NH4, and NO2 as predictors.

fm7 <- lm(Temp ~ CO2 + CH4 + NO2, climate)
cs7 <- coefficients(summary(fm7))

## Is CH4 significant at level 0.05 in the presence of the other predictors?

cs7["CH4", "Pr(>|t|)"]

## Is NO2 significant at level 0.05 in the presence of the other predictors?

cs7["NO2", "Pr(>|t|)"]

## Is CO2 significant at level 0.05 in the presence of the other predictors?

cs7["CO2", "Pr(>|t|)"]

## What is the estimated correlation between the estimated coefficients of NO2 and CO2?

cov2cor(vcov(fm7))
round(cov2cor(vcov(fm7))["NO2", "CO2"], 4)

## Report the coefficient of determination R2 for this model.

summary(fm7)$r.squared
round(summary(fm7)$r.squared, 4)


## Compute and report the predictive R2 for this model.

predR2 <- function(formula)
{
    n <- nrow(climate)
    e0 <- numeric(n)
    e <- numeric(n)
    y <- climate$Temp
    for (i in seq_len(n))
    {
        e0[i] <- y[i] - mean(climate$Temp[-i])
        fm <- lm(formula, climate[-i,])
        e[i] <- y[i] - predict(fm, climate[i,])
    }
    T2 <- sum(e0^2)
    S2 <- sum(e^2)
    (T2 - S2) / T2
}

round(predR2(Temp ~ CO2 + CH4 + NO2), 4)