arsenic <- read.table("C:/Rexamples/arsenic_data.txt", header=T) View(arsenic) summary(arsenic) attach(arsenic) regression1 <- lm(conc ~ dist) summary(regression1) anova(regression1) attach(regression1) output1 <- cbind(dist,conc,fitted.values,residuals) View(output1) summary(output1) regr1_clim <- predict(regression1, se.fit=TRUE, interval="confidence") regr1_plim <- predict(regression1, se.fit=TRUE, interval="prediction") regr1_estimates <- cbind(dist, regr1_clim$fit, regr1_plim$fit[,2:3], conc, residuals) colnames(regr1_estimates) <- c("dist","fit","lclm","uclm","lclp","uclp","conc","resid") # fitted values conf limits for the mean and for a prediction # lclm,uclm confidence interval for mean response # lclp,uclp prediction interval for one observation regr1_estimates matplot(dist, cbind(conc,fitted.values), pch="*", type="pl") matplot(dist, cbind(residuals,0), pch="*", type="pl") matplot(dist, regr1_estimates[,2:6], lty = c(1,2,2,3,3), type = "l") # estimate the parameters # intercept out_contrast1 <- data.frame(intercept=c(1), dist=c(0)) estimate_intercept <- predict(regression1, out_contrast1, se.fit=TRUE, interval="confidence") estimate_intercept # mean out_contrast2 <- data.frame(intercept=c(1), dist=c(16.1)) estimate_mean <- predict(regression1, out_contrast2, se.fit=TRUE, interval="confidence") estimate_mean # mean response dist=20 out_contrast3 <- data.frame(intercept=c(1), dist=c(20)) estimate_mean <- predict(regression1, out_contrast3, se.fit=TRUE, interval="confidence") estimate_mean --------------------- NOTE ------------------------ The first line above which reads the data from the file "arsenic_data.txt" assumes that this file is in the directory "C:/Rexamples". Obey this convention or modify the command. -------------------------- results ---------------------------------- The remainder of this file is copied and pasted from the RStudio command window. ------------------ EXTRAS --------------------------- ---here are some extra commands which may be of interest --------- plot(regression1) --------------------------------------------- ------------- output --------------------- > arsenic <- read.table("C:/Rexamples/arsenic_data.txt", header=T) > View(arsenic) > summary(arsenic) conc dist Min. :0.3000 Min. : 2.0 1st Qu.:0.8475 1st Qu.: 8.5 Median :1.5800 Median :13.5 Mean :1.6280 Mean :16.1 3rd Qu.:2.0000 3rd Qu.:22.5 Max. :3.2600 Max. :36.0 > attach(arsenic) > regression1 <- lm(conc ~ dist) > summary(regression1) Call: lm(formula = conc ~ dist) Residuals: Min 1Q Median 3Q Max -1.0847 -0.2487 0.1066 0.3088 0.6864 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.88623 0.31072 9.289 1.47e-05 *** dist -0.07815 0.01611 -4.850 0.00127 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.5409 on 8 degrees of freedom Multiple R-squared: 0.7462, Adjusted R-squared: 0.7145 F-statistic: 23.53 on 1 and 8 DF, p-value: 0.001271 > anova(regression1) Analysis of Variance Table Response: conc Df Sum Sq Mean Sq F value Pr(>F) dist 1 6.8826 6.8826 23.526 0.001271 ** Residuals 8 2.3404 0.2925 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > attach(regression1) > output1 <- cbind(dist,conc,fitted.values,residuals) > View(output1) > summary(output1) dist conc fitted.values Min. : 2.0 Min. :0.3000 Min. :0.0728 1st Qu.: 8.5 1st Qu.:0.8475 1st Qu.:1.1278 Median :13.5 Median :1.5800 Median :1.8312 Mean :16.1 Mean :1.6280 Mean :1.6280 3rd Qu.:22.5 3rd Qu.:2.0000 3rd Qu.:2.2219 Max. :36.0 Max. :3.2600 Max. :2.7299 residuals Min. :-1.0847 1st Qu.:-0.2487 Median : 0.1066 Mean : 0.0000 3rd Qu.: 0.3088 Max. : 0.6864 > regr1_clim <- predict(regression1, se.fit=TRUE, interval="confidence") > regr1_plim <- predict(regression1, se.fit=TRUE, interval="prediction") Warning message: In predict.lm(regression1, se.fit = TRUE, interval = "prediction") : predictions on current data refer to _future_ responses > regr1_estimates <- cbind(dist, regr1_clim$fit, regr1_plim$fit[,2:3], conc, residuals) > colnames(regr1_estimates) <- c("dist","fit","lclm","uclm","lclp","uclp","conc","resid") > # fitted values conf limits for the mean and for a prediction > # lclm,uclm confidence interval for mean response > # lclp,uclp prediction interval for one observation > regr1_estimates dist fit lclm uclm lclp uclp 1 2 2.72992457 2.0741642 3.3856849 1.32077860 4.139071 2 4 2.57362321 1.9755562 3.1716902 1.19038251 3.956864 3 8 2.26102050 1.7648947 2.7571463 0.91870477 3.603336 4 10 2.10471914 1.6498181 2.5596202 0.77708753 3.432351 5 12 1.94841778 1.5256022 2.3712334 0.63143496 3.265401 6 15 1.71396575 1.3174340 2.1104975 0.40518458 3.022747 7 21 1.24506167 0.8106511 1.6794723 -0.07568935 2.565813 8 23 1.08876032 0.6183434 1.5591772 -0.24426735 2.421788 9 30 0.54170556 -0.1081336 1.1915447 -0.86469468 1.948106 10 36 0.07280149 -0.7652043 0.9108073 -1.42983765 1.575441 conc resid 1 3.19 0.46007543 2 3.26 0.68637679 3 1.82 -0.44102050 4 1.02 -1.08471914 5 1.85 -0.09841778 6 2.05 0.33603425 7 1.34 0.09493833 8 0.79 -0.29876032 9 0.66 0.11829444 10 0.30 0.22719851 > matplot(dist, cbind(conc,fitted.values), pch="*", type="pl") > matplot(dist, cbind(residuals,0), pch="*", type="pl") > matplot(dist, regr1_estimates[,2:6], lty = c(1,2,2,3,3), type = "l") > # estimate the parameters > # intercept > out_contrast1 <- data.frame(intercept=c(1), dist=c(0)) > estimate_intercept <- predict(regression1, out_contrast1, se.fit=TRUE, interval="confidence") > estimate_intercept $fit fit lwr upr 1 2.886226 2.169704 3.602747 $se.fit [1] 0.31072 $df [1] 8 $residual.scale [1] 0.5408773 > # mean > out_contrast2 <- data.frame(intercept=c(1), dist=c(16.1)) > estimate_mean <- predict(regression1, out_contrast2, se.fit=TRUE, interval="confidence") > estimate_mean $fit fit lwr upr 1 1.628 1.23358 2.02242 $se.fit [1] 0.1710404 $df [1] 8 $residual.scale [1] 0.5408773 > # mean response dist=20 > out_contrast3 <- data.frame(intercept=c(1), dist=c(20)) > estimate_mean <- predict(regression1, out_contrast3, se.fit=TRUE, interval="confidence") > estimate_mean $fit fit lwr upr 1 1.323212 0.9030168 1.743408 $se.fit [1] 0.182218 $df [1] 8 $residual.scale [1] 0.5408773 >