Untitled

mail@pastecode.io avatar
unknown
plain_text
2 years ago
2.3 kB
12
Indexable
Never
#lab1.2

load("Data/Pb_Vasternorrland.rda")
summary(Pb_Vasternorrland)
head(Pb_Vasternorrland)

x = I(Pb_Vasternorrland$year - 1975)
Y <- Pb_Vasternorrland$Pb

plot(x,Y)

(Pb.model <- lm(Y ~ x, data = Pb_Vasternorrland))

## slope
# get summary
summary(Pb.model)
# confidens interval
confint(Pb.model, level = 0.95)
#coeficineter b0 b1 osv
coef(Pb.model)


## Predict ett år
# Predict the Pb value for the year 2011 with confidence intervals
newdata <- data.frame(x = 2011 - 1975)
prediction2011 <- predict(Pb.model, newdata = newdata, interval = "prediction")

# View the predicted value and confidence intervals
prediction2011

# Create a residual plot
plot(Pb.model, which = 1)

# Create a Q-Q plot
plot(Pb.model, which = 2)


## Part b

# Create a variable for the years from 1975
x = I(Pb_Vasternorrland$year - 1975)

# Create a variable for the natural logarithm of the Pb values
Y <- log(Pb_Vasternorrland$Pb)

# Fit a linear model with natural logarithm of Pb
Pb.model2 <- lm(Y ~ x, data = Pb_Vasternorrland)

## slope
# get summary
summary(Pb.model2)
# confidens interval
confint(Pb.model2, level = 0.95)
#coeficineter b0 b1 osv
coef(Pb.model2)


## Predict ett år
# Calculate the mean and standard deviation of the residuals from the log-linear model
residuals <- resid(Pb.model2)
residual_mean <- mean(residuals)
residual_sd <- sd(residuals)

# Calculate the 95% prediction interval for the logarithm of Pb in 2011
new_years <- data.frame(x = 2011 - 1975)
log_pred <- predict(Pb.model2, newdata = new_years, interval = "prediction", level = 0.95)

# Calculate the upper and lower limits of the prediction interval for the logarithm of Pb in 2011
log_pred_lower <- log_pred[1,2]
log_pred_upper <- log_pred[1,3]

# Calculate the upper and lower limits of the interval for the observed log-lead concentrations in 2011
observed_lower <- log_pred_lower - 2*residual_sd
observed_upper <- log_pred_upper + 2*residual_sd

# Create a residual plot
plot(Pb.model2, which = 1)

# Create a Q-Q plot
plot(Pb.model2, which = 2)


## Part C
summary(Pb.model2)
# confidens interval
confint(Pb.model2, level = 0.95)
#coeficineter b0 b1 osv
coef(Pb.model2)
a <- exp(3.13)
CI_low <- exp(2.93)
CI_high <- exp(3.33)

est <- exp(-0.081168)

95_low_slope <- exp(-0.08903472)

95_high_slope <- exp(-0.07330116)