Statistics for Laboratory Scientists ( 140.615 )
Linear Regression
Example from class: David Sullivan’s heme data
h2o2 <- rep(c(0,10,25,50),each=3)
pf3d7 <- c(0.3399,0.3563,0.3538,
0.3168,0.3054,0.3174,
0.2460,0.2618,0.2848,
0.1535,0.1613,0.1525)Plot the data, and add the least squares fit line.
par(las=1)
plot(h2o2, pf3d7, xlab="H2O2 concentration", ylab="OD")
abline(lsfit(h2o2, pf3d7), col="red", lty=2)
The regression of optical density on hydrogen peroxide concentration, using the built-in ‘lm’ function
lm.out <- lm(pf3d7 ~ h2o2)
lm.out##
## Call:
## lm(formula = pf3d7 ~ h2o2)
##
## Coefficients:
## (Intercept) h2o2
## 0.353050 -0.003871lm.sum <- summary(lm.out)
lm.sum##
## Call:
## lm(formula = pf3d7 ~ h2o2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.013150 -0.007486 0.001275 0.003107 0.028525
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3530501 0.0049955 70.67 7.84e-15 ***
## h2o2 -0.0038710 0.0001759 -22.00 8.43e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01148 on 10 degrees of freedom
## Multiple R-squared: 0.9798, Adjusted R-squared: 0.9777
## F-statistic: 484.1 on 1 and 10 DF, p-value: 8.426e-10attributes(lm.sum)## $names
## [1] "call" "terms" "residuals" "coefficients" "aliased" "sigma"
## [7] "df" "r.squared" "adj.r.squared" "fstatistic" "cov.unscaled"
##
## $class
## [1] "summary.lm"The table of coefficients.
lm.sum$coef## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.353050073 0.0049955194 70.67335 7.844454e-15
## h2o2 -0.003870984 0.0001759324 -22.00268 8.426407e-10The parameter estimates.
lm.sum$coef[,1]## (Intercept) h2o2
## 0.353050073 -0.003870984The residual standard deviation.
lm.sum$sigma## [1] 0.01147782The coefficient of determination.
lm.sum$r.squared## [1] 0.9797619To see how to derive these statistics from scratch, please see the advanced code. It is simple arithmetic following the lecture notes and therefore not really advanced code, I just tucked it away there.
Also, to avoid data glitches, it is a good habit to make a data frame.
dat <- data.frame(conc=h2o2, od=pf3d7)
dat## conc od
## 1 0 0.3399
## 2 0 0.3563
## 3 0 0.3538
## 4 10 0.3168
## 5 10 0.3054
## 6 10 0.3174
## 7 25 0.2460
## 8 25 0.2618
## 9 25 0.2848
## 10 50 0.1535
## 11 50 0.1613
## 12 50 0.1525str(dat)## 'data.frame': 12 obs. of 2 variables:
## $ conc: num 0 0 0 10 10 10 25 25 25 50 ...
## $ od : num 0.34 0.356 0.354 0.317 0.305 ...The linear regression model.
lm.out <- lm(od ~ conc, data=dat)
summary(lm.out)##
## Call:
## lm(formula = od ~ conc, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.013150 -0.007486 0.001275 0.003107 0.028525
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3530501 0.0049955 70.67 7.84e-15 ***
## conc -0.0038710 0.0001759 -22.00 8.43e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01148 on 10 degrees of freedom
## Multiple R-squared: 0.9798, Adjusted R-squared: 0.9777
## F-statistic: 484.1 on 1 and 10 DF, p-value: 8.426e-10Plot of the data with least squares regression line..
par(las=1)
plot(dat$conc, dat$od, xlab="H2O2 concentration", ylab="OD")
abline(lm.out, col="red", lty=2)
Confidence intervals for the regression parameters
confint(lm.out)## 2.5 % 97.5 %
## (Intercept) 0.341919363 0.364180784
## conc -0.004262986 -0.003478982confint(lm.out,level=0.99)## 0.5 % 99.5 %
## (Intercept) 0.337217910 0.368882237
## conc -0.004428562 -0.003313406To see how to derive these from scratch, also see the advanced code.
Checking model assumptions
plot(lm.out,which=1)
plot(lm.out,which=2)