Untitled
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)