# Hinton's study time data again
studytimedata = data.frame(study.time=c(40,43,18,10,25,33,27,17,30,47),
exam.score=c(58,73,56,47,58,54,45,32,68,69))
# some guesses as to the regression line
plot(studytimedata[,1:2])
# we want a linear relation between study time and exam score.
# a line can be determined by its offset a (y-value at x=0)
# and its slope b
my.line = function(x, a, b) {
a + b*x
}
# here are some guesses
xvals = 0:50
lines(xvals, my.line(xvals, 25, 1.1), col="blue")
lines(xvals, my.line(xvals, 23, 1.2), col="yellowgreen")
lines(xvals, my.line(xvals, 24, 1.05), col="green")
lines(xvals, my.line(xvals, 27, 1.003), col="turquoise1")
# regression model Y = A + B*X
# finding the best line
lm(exam.score ~ study.time, data = studytimedata)
# closer inspection
lm.obj = lm(exam.score ~ study.time, data = studytimedata)
summary(lm.obj)
# we can plot the line using the function my.line that we defined earlier
lines(xvals, my.line(xvals, coef(lm.obj)[1], coef(lm.obj)[2]), col = "red")
# or, more simply, using abline(),
# which is able to extract the coefficients from the lm.obj object
abline(lm.obj, col = "tomato")
###################
# Reading the output of your lm model
summary(lm.obj)
# Call:
# lm(formula = exam.score ~ study.time, data = studytimedata)
# Residuals:
# Min 1Q Median 3Q Max
# -15.064 -5.888 2.288 6.218 11.255
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 34.4057 7.8925 4.359 0.00242 **
# study.time 0.7446 0.2532 2.941 0.01869 *
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# Residual standard error: 9.144 on 8 degrees of freedom
# Multiple R-squared: 0.5194, Adjusted R-squared: 0.4594
# F-statistic: 8.647 on 1 and 8 DF, p-value: 0.01869
# Accessing the information directly:
names(lm.obj)
# [1] "coefficients" "residuals" "effects" "rank" "fitted.values"
# [6] "assign" "qr" "df.residual" "xlevels" "call"
# [11] "terms" "model"
# coefficients: the beta values for intercept and slope.
# residuals: difference between observed and predicted values
# fitted.values: predicted values
###
# Coefficients: intercept and slope
coef(lm.obj)
# or equivalently
lm.obj$coefficients
# accessing the intercept:
lm.obj$coefficients[1]
# accessing the slope:
lm.obj$coefficients[2]
###
# Inspecting the predicted values:
# three ways to access them
lm.obj$fitted.values
fitted(lm.obj)
predict(lm.obj)
###
# Inspecting residuals
#
# Linear regression assumes that we have:
# observed values = fitted values + error
# where the error (the residuals) is normally distributed with mean of zero
# residuals = observed - fitted
# so the following two commands will yield the same results:
lm.obj$residuals
studytimedata$exam.score - lm.obj$fitted.values
# inspecting residuals: are they approximately normally distributed with a mean of zero?
quantile(lm.obj$residuals)
plot(density(lm.obj$residuals)
#####
# For each coefficient, a t-test is done
# Null hypothesis: The slope is really zero.
# If we can reject the null hypothesis, this means
# that the predictor in question is in fact correlated with the
# dependent variable.
#
# t value:
# t.value = (beta - nullhypothesis_value) / standarderror
# = (beta - 0) / standarderror
#
# In the table describing the coefficients,
# beta is the first column (Estimate), and the
# standard error is the second column (Std. Error)
lmsum.obj = summary(lm.obj)
# looking more into the summary object
names(lmsum.obj)
# [1] "call" "terms" "residuals" "coefficients" "aliased"
# [6] "sigma" "df" "r.squared" "adj.r.squared" "fstatistic"
# [11] "cov.unscaled"
#
# particularly useful: r.squared, adj.r.squared, sigma (the residual standard error)
# and fstatistic
beta = lmsum.obj$coefficients[2,1]
stderror = lmsum.obj$coefficients[2,2]
t.value = (beta - 0) / stderror
# degrees of freedom: number of data points minus 2
df = length(studytimedata$exam.score) - 2
# and do a two-tailed t-test, which here comes out as significant:
2 * pt(t.value, df = df, lower.tail = F)
# Below the coefficients, we have the residual standard error
# This is the variability in the dependent variable not
# handled through predictors. This is
residual.se = sqrt(sum(lm.obj$residuals**2) / df)
# Below that is R-squared: this is Pearson's r squared, indicating the fraction of
# variance in the data explained by the model.
# This is an indicator of the effect size for this linear regression model.
R.squared = cor(studytimedata$study.time, studytimedata$exam.score)**2
# Adjusted R-squared: adjusted for number of predictors.
# The more predictors we add, the better R-squared, but some of the
# improvement may be due to chance. Adjusted R-squared tries to
# correct for this.
###
# Below that, an F-test. The null hypothesis is that all coefficients (except the intercept)
# are equal to zero.
# F = systematic_variance / error_variance
# systematic_variance = model_sum_of_squares / model_df
# error_variance = error_sum_of_squares / error_df
# overall sum of squares: squared distances between observed values and mean observed value
SSQ.overall = sum((studytimedata$exam.score - mean(studytimedata$exam.score))^2)
# model sum of squares: squared distances between predicted (fitted) values and mean observed value
SSQ.model = sum((lm.obj$fitted.values - mean(studytimedata$exam.score))^2)
# error sum of squares: summed squared residuals
SSQ.error = sum(lm.obj$residuals^2)
# total df: number of datapoints minus one
# model df: num parameters (coefficient plus intercept) - 1
DFtotal = length(studytimedata$study.time) - 1
DFmodel = 1
DFerror = DFtotal - DFmodel
# variances: first model, then error
variance.model = SSQ.model / DFmodel
variance.error = SSQ.error / DFerror
# F statistic
F.value = variance.model / variance.error
# F test
1 - pf(F.value, DFmodel, DFerror, lower.tail = F)
##############
# Hinton, problems with regression: "smiling" example
hinton.smile = data.frame(smile.time = c(0.4, 0.8, 0.8, 1.2, 1.4, 1.8, 2.2, 2.6, 3.0),
sold = c(16, 12, 20, 16, 34, 30, 26, 22, 38))
plot(hinton.smile)
# a significant correlation when all items are taken into account
cor(hinton.smile$smile.time, hinton.smile$sold)
# no significant correlation when we leave out the last salesperson
cor(hinton.smile[-c(9),]$smile.time, hinton.smile[-c(9),]$sold)
# no correlation at all for the first 4 salespeople
cor(hinton.smile[1:4,]$smile.time, hinton.smile[1:4,]$sold)
# perfect negative correlation for participants 5-8
cor(hinton.smile[5:8,]$smile.time, hinton.smile[5:8,]$sold)