How to use this document. Run the code along with
me. Comments start with #. Each section mirrors the live
lab; code chunks show the command and its output together.
Today’s plan
- Clean environment + setup
- Load & inspect the
Prestige data
- Visualize distributions
- Linear regression:
prestige ~ education
- Confidence intervals for the coefficients
- Quick practice
- Wrap-up
Clean Environment +
Setup
rm(list = ls()) # clean environment
# Load the packages
library(ggplot2)
library(dplyr)
library(carData) # contains the Prestige dataset
Load & Inspect the
Prestige Data
data("Prestige")
# In RStudio you would run View(Prestige); here we preview the first rows.
knitr::kable(head(Prestige), caption = "First rows of the Prestige data")
First rows of the Prestige data
| gov.administrators |
13.11 |
12351 |
11.16 |
68.8 |
1113 |
prof |
| general.managers |
12.26 |
25879 |
4.02 |
69.1 |
1130 |
prof |
| accountants |
12.77 |
9271 |
15.70 |
63.4 |
1171 |
prof |
| purchasing.officers |
11.42 |
8865 |
9.11 |
56.8 |
1175 |
prof |
| chemists |
14.62 |
8403 |
11.68 |
73.5 |
2111 |
prof |
| physicists |
15.64 |
11030 |
5.13 |
77.6 |
2113 |
prof |
#> education income women prestige
#> Min. : 6.380 Min. : 611 Min. : 0.000 Min. :14.80
#> 1st Qu.: 8.445 1st Qu.: 4106 1st Qu.: 3.592 1st Qu.:35.23
#> Median :10.540 Median : 5930 Median :13.600 Median :43.60
#> Mean :10.738 Mean : 6798 Mean :28.979 Mean :46.83
#> 3rd Qu.:12.648 3rd Qu.: 8187 3rd Qu.:52.203 3rd Qu.:59.27
#> Max. :15.970 Max. :25879 Max. :97.510 Max. :87.20
#> census type
#> Min. :1113 bc :44
#> 1st Qu.:3120 prof:31
#> Median :5135 wc :23
#> Mean :5402 NA's: 4
#> 3rd Qu.:8312
#> Max. :9517
Variable notes
income = average income (in Canadian dollars)
education = average years of education
women = % of women in the occupation
prestige = social prestige score (our dependent
variable)
Visualize
Distributions
Histogram of prestige (default binning).
ggplot(Prestige, aes(x = prestige)) +
geom_histogram(color = "pink") +
labs(title = "Distribution of Occupational Prestige",
x = "Prestige", y = "Frequency")

Histogram of prestige (fixed bin width).
ggplot(Prestige, aes(x = prestige)) +
geom_histogram(binwidth = 5, fill = "skyblue", color = "black") +
labs(title = "Distribution of Occupational Prestige",
x = "Prestige", y = "Frequency")

Histogram of education.
ggplot(Prestige, aes(x = education)) +
geom_histogram(binwidth = 1, fill = "lightgreen", color = "black") +
labs(title = "Distribution of Education Years",
x = "Education (Years)", y = "Frequency")

Linear Regression:
prestige ~ education
Predict the social prestige score from years of education.
model <- lm(prestige ~ education, data = Prestige)
summary(model)
#>
#> Call:
#> lm(formula = prestige ~ education, data = Prestige)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -26.0397 -6.5228 0.6611 6.7430 18.1636
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -10.732 3.677 -2.919 0.00434 **
#> education 5.361 0.332 16.148 < 2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 9.103 on 100 degrees of freedom
#> Multiple R-squared: 0.7228, Adjusted R-squared: 0.72
#> F-statistic: 260.8 on 1 and 100 DF, p-value: < 2.2e-16
- Intercept = predicted prestige when
education = 0.
- education = expected change in prestige for a
1-year increase in education.
Confidence Intervals
(95%)
confint(model, level = 0.95)
#> 2.5 % 97.5 %
#> (Intercept) -18.027220 -3.436744
#> education 4.702223 6.019533
# Lower and upper bounds for each coefficient
Quick Practice
P1) Add income as another predictor.
Include multiple predictors by connecting them with +.
Compare coefficients and p-values.
model2 <- lm(prestige ~ education + income, data = Prestige)
summary(model2)
#>
#> Call:
#> lm(formula = prestige ~ education + income, data = Prestige)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -19.4040 -5.3308 0.0154 4.9803 17.6889
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -6.8477787 3.2189771 -2.127 0.0359 *
#> education 4.1374444 0.3489120 11.858 < 2e-16 ***
#> income 0.0013612 0.0002242 6.071 2.36e-08 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 7.81 on 99 degrees of freedom
#> Multiple R-squared: 0.798, Adjusted R-squared: 0.7939
#> F-statistic: 195.6 on 2 and 99 DF, p-value: < 2.2e-16
P2) Using confint, which predictors have CIs
that include 0?
confint(model2, level = 0.95)
#> 2.5 % 97.5 %
#> (Intercept) -1.323493e+01 -0.460629799
#> education 3.445127e+00 4.829761535
#> income 9.162805e-04 0.001806051
P3) Try predicting prestige using the % of
women.
model_women <- lm(prestige ~ women, data = Prestige)
summary(model_women)
#>
#> Call:
#> lm(formula = prestige ~ women, data = Prestige)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -33.444 -12.391 -4.126 13.034 39.185
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 48.69300 2.30760 21.101 <2e-16 ***
#> women -0.06417 0.05385 -1.192 0.236
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 17.17 on 100 degrees of freedom
#> Multiple R-squared: 0.014, Adjusted R-squared: 0.004143
#> F-statistic: 1.42 on 1 and 100 DF, p-value: 0.2362
confint(model_women, level = 0.95)
#> 2.5 % 97.5 %
#> (Intercept) 44.114782 53.27121657
#> women -0.171008 0.04266233
Wrap-up
Key takeaways
confint() gives the CI for regression
coefficients.
- p-values tell us whether relationships are statistically
significant.
- Education and income often predict social prestige — a core social
science insight.
- A CI that includes 0 means the predictor is not statistically
distinguishable from “no effect.”
- Next week: one-sample t-test for group means.
LS0tCnRpdGxlOiAiTGFiIDYg4oCUIENvbmZpZGVuY2UgSW50ZXJ2YWxzICYgUmVncmVzc2lvbiBCYXNpY3MiCnN1YnRpdGxlOiAiUXVhbnRpdGF0aXZlIFJlYXNvbmluZyDCtyBMQUIgNDEzIgphdXRob3I6ICJJbnN0cnVjdG9yOiBTdWJpbiBOYSIKZGF0ZTogIk9jdG9iZXIgMTUsIDIwMjUiCm91dHB1dDoKICBodG1sX2RvY3VtZW50OgogICAgdGhlbWU6IGZsYXRseQogICAgaGlnaGxpZ2h0OiB0YW5nbwogICAgdG9jOiB0cnVlCiAgICB0b2NfZmxvYXQ6IHRydWUKICAgIHRvY19kZXB0aDogMgogICAgbnVtYmVyX3NlY3Rpb25zOiB0cnVlCiAgICBkZl9wcmludDogcGFnZWQKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKLS0tCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KAogIGVjaG8gPSBUUlVFLCBtZXNzYWdlID0gRkFMU0UsIHdhcm5pbmcgPSBGQUxTRSwKICBmaWcuYWxpZ24gPSAiY2VudGVyIiwgZmlnLndpZHRoID0gNywgZmlnLmhlaWdodCA9IDQuMiwKICBjb21tZW50ID0gIiM+IgopCmBgYAoKPiAqKkhvdyB0byB1c2UgdGhpcyBkb2N1bWVudC4qKiBSdW4gdGhlIGNvZGUgYWxvbmcgd2l0aCBtZS4gQ29tbWVudHMgc3RhcnQgd2l0aCBgI2AuCj4gRWFjaCBzZWN0aW9uIG1pcnJvcnMgdGhlIGxpdmUgbGFiOyBjb2RlIGNodW5rcyBzaG93IHRoZSBjb21tYW5kIGFuZCBpdHMgb3V0cHV0IHRvZ2V0aGVyLgoKKipUb2RheSdzIHBsYW4qKgoKMS4gQ2xlYW4gZW52aXJvbm1lbnQgKyBzZXR1cAoyLiBMb2FkICYgaW5zcGVjdCB0aGUgYFByZXN0aWdlYCBkYXRhCjMuIFZpc3VhbGl6ZSBkaXN0cmlidXRpb25zCjQuIExpbmVhciByZWdyZXNzaW9uOiBgcHJlc3RpZ2UgfiBlZHVjYXRpb25gCjUuIENvbmZpZGVuY2UgaW50ZXJ2YWxzIGZvciB0aGUgY29lZmZpY2llbnRzCjYuIFF1aWNrIHByYWN0aWNlCjcuIFdyYXAtdXAKCi0tLQoKIyBDbGVhbiBFbnZpcm9ubWVudCArIFNldHVwCgpgYGB7ciBjbGVhbi1zZXR1cH0Kcm0obGlzdCA9IGxzKCkpICAgIyBjbGVhbiBlbnZpcm9ubWVudAoKIyBMb2FkIHRoZSBwYWNrYWdlcwpsaWJyYXJ5KGdncGxvdDIpCmxpYnJhcnkoZHBseXIpCmxpYnJhcnkoY2FyRGF0YSkgICMgY29udGFpbnMgdGhlIFByZXN0aWdlIGRhdGFzZXQKYGBgCgojIExvYWQgJiBJbnNwZWN0IHRoZSBQcmVzdGlnZSBEYXRhCgpgYGB7ciBpbXBvcnR9CmRhdGEoIlByZXN0aWdlIikKCiMgSW4gUlN0dWRpbyB5b3Ugd291bGQgcnVuIFZpZXcoUHJlc3RpZ2UpOyBoZXJlIHdlIHByZXZpZXcgdGhlIGZpcnN0IHJvd3MuCmtuaXRyOjprYWJsZShoZWFkKFByZXN0aWdlKSwgY2FwdGlvbiA9ICJGaXJzdCByb3dzIG9mIHRoZSBQcmVzdGlnZSBkYXRhIikKCnN1bW1hcnkoUHJlc3RpZ2UpCmBgYAoKPiAqKlZhcmlhYmxlIG5vdGVzKioKPgo+IC0gYGluY29tZWAgPSBhdmVyYWdlIGluY29tZSAoaW4gQ2FuYWRpYW4gZG9sbGFycykKPiAtIGBlZHVjYXRpb25gID0gYXZlcmFnZSB5ZWFycyBvZiBlZHVjYXRpb24KPiAtIGB3b21lbmAgPSAlIG9mIHdvbWVuIGluIHRoZSBvY2N1cGF0aW9uCj4gLSBgcHJlc3RpZ2VgID0gc29jaWFsIHByZXN0aWdlIHNjb3JlIChvdXIgZGVwZW5kZW50IHZhcmlhYmxlKQoKIyBWaXN1YWxpemUgRGlzdHJpYnV0aW9ucwoKKipIaXN0b2dyYW0gb2YgcHJlc3RpZ2UgKGRlZmF1bHQgYmlubmluZykuKioKCmBgYHtyIGhpc3QtcHJlc3RpZ2UtZGVmYXVsdH0KZ2dwbG90KFByZXN0aWdlLCBhZXMoeCA9IHByZXN0aWdlKSkgKwogIGdlb21faGlzdG9ncmFtKGNvbG9yID0gInBpbmsiKSArCiAgbGFicyh0aXRsZSA9ICJEaXN0cmlidXRpb24gb2YgT2NjdXBhdGlvbmFsIFByZXN0aWdlIiwKICAgICAgIHggPSAiUHJlc3RpZ2UiLCB5ID0gIkZyZXF1ZW5jeSIpCmBgYAoKKipIaXN0b2dyYW0gb2YgcHJlc3RpZ2UgKGZpeGVkIGJpbiB3aWR0aCkuKioKCmBgYHtyIGhpc3QtcHJlc3RpZ2V9CmdncGxvdChQcmVzdGlnZSwgYWVzKHggPSBwcmVzdGlnZSkpICsKICBnZW9tX2hpc3RvZ3JhbShiaW53aWR0aCA9IDUsIGZpbGwgPSAic2t5Ymx1ZSIsIGNvbG9yID0gImJsYWNrIikgKwogIGxhYnModGl0bGUgPSAiRGlzdHJpYnV0aW9uIG9mIE9jY3VwYXRpb25hbCBQcmVzdGlnZSIsCiAgICAgICB4ID0gIlByZXN0aWdlIiwgeSA9ICJGcmVxdWVuY3kiKQpgYGAKCioqSGlzdG9ncmFtIG9mIGVkdWNhdGlvbi4qKgoKYGBge3IgaGlzdC1lZHVjYXRpb259CmdncGxvdChQcmVzdGlnZSwgYWVzKHggPSBlZHVjYXRpb24pKSArCiAgZ2VvbV9oaXN0b2dyYW0oYmlud2lkdGggPSAxLCBmaWxsID0gImxpZ2h0Z3JlZW4iLCBjb2xvciA9ICJibGFjayIpICsKICBsYWJzKHRpdGxlID0gIkRpc3RyaWJ1dGlvbiBvZiBFZHVjYXRpb24gWWVhcnMiLAogICAgICAgeCA9ICJFZHVjYXRpb24gKFllYXJzKSIsIHkgPSAiRnJlcXVlbmN5IikKYGBgCgojIExpbmVhciBSZWdyZXNzaW9uOiBgcHJlc3RpZ2UgfiBlZHVjYXRpb25gCgpQcmVkaWN0IHRoZSBzb2NpYWwgcHJlc3RpZ2Ugc2NvcmUgZnJvbSB5ZWFycyBvZiBlZHVjYXRpb24uCgpgYGB7ciBtb2RlbH0KbW9kZWwgPC0gbG0ocHJlc3RpZ2UgfiBlZHVjYXRpb24sIGRhdGEgPSBQcmVzdGlnZSkKc3VtbWFyeShtb2RlbCkKYGBgCgo+IC0gKipJbnRlcmNlcHQqKiA9IHByZWRpY3RlZCBwcmVzdGlnZSB3aGVuIGBlZHVjYXRpb25gID0gMC4KPiAtICoqZWR1Y2F0aW9uKiogPSBleHBlY3RlZCBjaGFuZ2UgaW4gcHJlc3RpZ2UgZm9yIGEgMS15ZWFyIGluY3JlYXNlIGluIGVkdWNhdGlvbi4KCiMgQ29uZmlkZW5jZSBJbnRlcnZhbHMgKDk1JSkKCmBgYHtyIGNvbmZpbnR9CmNvbmZpbnQobW9kZWwsIGxldmVsID0gMC45NSkKIyBMb3dlciBhbmQgdXBwZXIgYm91bmRzIGZvciBlYWNoIGNvZWZmaWNpZW50CmBgYAoKIyBRdWljayBQcmFjdGljZQoKKipQMSkgQWRkIGBpbmNvbWVgIGFzIGFub3RoZXIgcHJlZGljdG9yLioqCkluY2x1ZGUgbXVsdGlwbGUgcHJlZGljdG9ycyBieSBjb25uZWN0aW5nIHRoZW0gd2l0aCBgK2AuIENvbXBhcmUgY29lZmZpY2llbnRzIGFuZCAqcCotdmFsdWVzLgoKYGBge3IgcDF9Cm1vZGVsMiA8LSBsbShwcmVzdGlnZSB+IGVkdWNhdGlvbiArIGluY29tZSwgZGF0YSA9IFByZXN0aWdlKQpzdW1tYXJ5KG1vZGVsMikKYGBgCgoqKlAyKSBVc2luZyBgY29uZmludGAsIHdoaWNoIHByZWRpY3RvcnMgaGF2ZSBDSXMgdGhhdCBpbmNsdWRlIDA/KioKCmBgYHtyIHAyfQpjb25maW50KG1vZGVsMiwgbGV2ZWwgPSAwLjk1KQpgYGAKCioqUDMpIFRyeSBwcmVkaWN0aW5nIHByZXN0aWdlIHVzaW5nIHRoZSAlIG9mIHdvbWVuLioqCgpgYGB7ciBwM30KbW9kZWxfd29tZW4gPC0gbG0ocHJlc3RpZ2UgfiB3b21lbiwgZGF0YSA9IFByZXN0aWdlKQpzdW1tYXJ5KG1vZGVsX3dvbWVuKQpjb25maW50KG1vZGVsX3dvbWVuLCBsZXZlbCA9IDAuOTUpCmBgYAoKIyBXcmFwLXVwCgoqKktleSB0YWtlYXdheXMqKgoKLSBgY29uZmludCgpYCBnaXZlcyB0aGUgQ0kgZm9yIHJlZ3Jlc3Npb24gY29lZmZpY2llbnRzLgotICpwKi12YWx1ZXMgdGVsbCB1cyB3aGV0aGVyIHJlbGF0aW9uc2hpcHMgYXJlIHN0YXRpc3RpY2FsbHkgc2lnbmlmaWNhbnQuCi0gRWR1Y2F0aW9uIGFuZCBpbmNvbWUgb2Z0ZW4gcHJlZGljdCBzb2NpYWwgcHJlc3RpZ2Ug4oCUIGEgY29yZSBzb2NpYWwgc2NpZW5jZSBpbnNpZ2h0LgotIEEgQ0kgdGhhdCBpbmNsdWRlcyAwIG1lYW5zIHRoZSBwcmVkaWN0b3IgaXMgbm90IHN0YXRpc3RpY2FsbHkgZGlzdGluZ3Vpc2hhYmxlIGZyb20gIm5vIGVmZmVjdC4iCi0gTmV4dCB3ZWVrOiBvbmUtc2FtcGxlICp0Ki10ZXN0IGZvciBncm91cCBtZWFucy4K