wheatear0 <- read.table("C:/Rexamples/wheateardata.txt", header=T) # sort wheatear0 so that x=tcell is increasing (for plotting convenience) xorder <- order(wheatear0$tcell) wheatear <- wheatear0[xorder,] View(wheatear) summary(wheatear) attach(wheatear) regression1 <- lm(mass ~ tcell) summary(regression1) anova(regression1) attach(regression1) output1 <- cbind(tcell,mass,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(tcell, regr1_clim$fit, regr1_plim$fit[,2:3], mass, residuals) colnames(regr1_estimates) <- c("tcell","fit","lclm","uclm","lclp","uclp","mass","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(tcell, cbind(mass,fitted.values), pch="*", type="pl") matplot(tcell, cbind(residuals,0), pch="*", type="pl") matplot(tcell, regr1_estimates[,2:6], lty = c(1,2,2,3,3), type = "l") # estimate the parameters # intercept out_contrast1 <- data.frame(intercept=c(1), tcell=c(0)) estimate_intercept <- predict(regression1, out_contrast1, se.fit=TRUE, interval="confidence") estimate_intercept # mean out_contrast2 <- data.frame(intercept=c(1), tcell=c(.32395)) estimate_mean <- predict(regression1, out_contrast2, se.fit=TRUE, interval="confidence") estimate_mean --------------------- NOTE ------------------------ The first line above which reads the data from the file "wheateardata.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. -------------------------------------------------------------------- > wheatear0 <- read.table("C:/Rexamples/wheateardata.txt", header=T) > # sort wheatear0 so that x=tcell is increasing (for plotting convenience) > xorder <- order(wheatear0$tcell) > wheatear <- wheatear0[xorder,] > View(wheatear) > summary(wheatear) mass tcell Min. :3.330 Min. :0.183 1st Qu.:6.290 1st Qu.:0.251 Median :6.810 Median :0.312 Mean :7.204 Mean :0.324 3rd Qu.:8.180 3rd Qu.:0.411 Max. :9.950 Max. :0.508 > attach(wheatear) > regression1 <- lm(mass ~ tcell) > summary(regression1) Call: lm(formula = mass ~ tcell) Residuals: Min 1Q Median 3Q Max -3.1429 -0.7327 0.3448 0.7472 3.2736 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.911 1.112 3.517 0.00230 ** tcell 10.165 3.296 3.084 0.00611 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.426 on 19 degrees of freedom Multiple R-squared: 0.3336, Adjusted R-squared: 0.2986 F-statistic: 9.513 on 1 and 19 DF, p-value: 0.006105 > anova(regression1) Analysis of Variance Table Response: mass Df Sum Sq Mean Sq F value Pr(>F) tcell 1 19.339 19.340 9.5129 0.006105 ** Residuals 19 38.627 2.033 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > attach(regression1) > output1 <- cbind(tcell,mass,fitted.values,residuals) > View(output1) > summary(output1) tcell mass fitted.values Min. :0.183 Min. :3.330 Min. :5.771 1st Qu.:0.251 1st Qu.:6.290 1st Qu.:6.463 Median :0.312 Median :6.810 Median :7.083 Mean :0.324 Mean :7.204 Mean :7.204 3rd Qu.:0.411 3rd Qu.:8.180 3rd Qu.:8.089 Max. :0.508 Max. :9.950 Max. :9.075 residuals Min. :-3.1429 1st Qu.:-0.7327 Median : 0.3448 Mean : 0.0000 3rd Qu.: 0.7472 Max. : 3.2736 > 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(tcell, regr1_clim$fit, regr1_plim$fit[,2:3], mass, residuals) > colnames(regr1_estimates) <- c("tcell","fit","lclm","uclm","lclp","uclp","mass","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 tcell fit lclm uclm lclp uclp mass 1 0.183 5.771488 4.601242 6.941733 2.565953 8.977022 6.12 2 0.203 5.974790 4.916383 7.033197 2.808372 9.141208 6.51 3 0.213 6.076441 5.071516 7.081366 2.927497 9.225386 6.29 4 0.213 6.076441 5.071516 7.081366 2.927497 9.225386 9.35 5 0.251 6.462716 5.639710 7.285722 3.367023 9.558409 5.43 6 0.251 6.462716 5.639710 7.285722 3.367023 9.558409 5.73 7 0.252 6.472881 5.654075 7.291687 3.378302 9.567460 3.33 8 0.252 6.472881 5.654075 7.291687 3.378302 9.567460 6.65 9 0.263 6.584697 5.809534 7.359861 3.501379 9.668016 4.62 10 0.304 7.001467 6.335857 7.667078 3.943852 10.059083 8.02 11 0.312 7.082788 6.426365 7.739212 4.027159 10.138417 7.83 12 0.332 7.286091 6.632504 7.939678 4.231070 10.341112 6.45 13 0.342 7.387742 6.724724 8.050760 4.330690 10.444794 6.75 14 0.370 7.672365 6.947803 8.396928 4.601377 10.743353 8.06 15 0.381 7.784182 7.023292 8.545071 4.704420 10.863943 8.18 16 0.411 8.089135 7.203329 8.974941 4.976158 11.202113 9.95 17 0.430 8.282273 7.302870 9.261675 5.141380 11.423165 9.08 18 0.430 8.282273 7.302870 9.261675 5.141380 11.423165 9.15 19 0.431 8.292438 7.307872 9.277003 5.149932 11.434944 7.56 20 0.471 8.699043 7.493636 9.904449 5.480506 11.917580 6.81 21 0.508 9.075152 7.648290 10.502014 5.767296 12.383008 9.42 resid 1 0.3485125 2 0.5352100 3 0.2135588 4 3.2735588 5 -1.0327158 6 -0.7327158 7 -3.1428810 8 0.1771190 9 -1.9646973 10 1.0185327 11 0.7472117 12 -0.8360907 13 -0.6377420 14 0.3876346 15 0.3958183 16 1.8608646 17 0.7977273 18 0.8677273 19 -0.7324378 20 -1.8890427 21 0.3448477 > matplot(tcell, cbind(mass,fitted.values), pch="*", type="pl") > matplot(tcell, cbind(residuals,0), pch="*", type="pl") > matplot(tcell, regr1_estimates[,2:6], lty = c(1,2,2,3,3), type = "l") > # estimate the parameters > # intercept > out_contrast1 <- data.frame(intercept=c(1), tcell=c(0)) > estimate_intercept <- predict(regression1, out_contrast1, se.fit=TRUE, interval="confidence") > estimate_intercept $fit fit lwr upr 1 3.91127 1.583651 6.23889 $se.fit [1] 1.112084 $df [1] 19 $residual.scale [1] 1.425826 > # mean > out_contrast2 <- data.frame(intercept=c(1), tcell=c(.32395)) > estimate_mean <- predict(regression1, out_contrast2, se.fit=TRUE, interval="confidence") > estimate_mean $fit fit lwr upr 1 7.204262 6.553036 7.855487 $se.fit [1] 0.3111408 $df [1] 19 $residual.scale [1] 1.425826 >