It is often helpful to start your R session by setting your working directory so you don't have to type the full path names to your data and other files
In [ ]:
# set the working directory
# setwd("~/Desktop/Rstatistics")
# setwd("C:/Users/dataclass/Desktop/Rstatistics")
You might also start by listing the files in your working directory
In [ ]:
getwd() # where am I?
list.files("dataSets") # files in the dataSets folder
The states.dta data comes from http://anawida.de/teach/SS14/anawida/4.linReg/data/states.dta.txt and appears to have originally appeared in Statistics with Stata by Lawrence C. Hamilton.
In [ ]:
# read the states data
states.data <- readRDS("dataSets/states.rds")
#get labels
states.info <- data.frame(attributes(states.data)[c("names", "var.labels")])
#look at last few labels
tail(states.info, 8)
In [ ]:
# summary of expense and csat columns, all rows
sts.ex.sat <- subset(states.data, select = c("expense", "csat"))
summary(sts.ex.sat)
# correlation between expense and csat
cor(sts.ex.sat)
In [ ]:
# scatter plot of expense vs csat
plot(sts.ex.sat)
In [ ]:
# Fit our regression model
sat.mod <- lm(csat ~ expense, # regression formula
data=states.data) # data set
# Summarize and print the results
summary(sat.mod) # show regression coefficients table
Many people find it surprising that the per-capita expenditure on students is negatively related to SAT scores. The beauty of multiple regression is that we can try to pull these apart. What would the association between expense and SAT scores be if there were no difference among the states in the percentage of students taking the SAT?
In [ ]:
summary(lm(csat ~ expense + percent, data = states.data))
In [ ]:
class(sat.mod)
names(sat.mod)
methods(class = class(sat.mod))[1:9]
In [ ]:
confint(sat.mod)
# hist(residuals(sat.mod))
In [ ]:
par(mar = c(4, 4, 2, 2), mfrow = c(1, 2)) #optional
plot(sat.mod, which = c(1, 2)) # "which" argument optional
In [ ]:
# fit another model, adding house and senate as predictors
sat.voting.mod <- lm(csat ~ expense + house + senate,
data = na.omit(states.data))
sat.mod <- update(sat.mod, data=na.omit(states.data))
# compare using the anova() function
anova(sat.mod, sat.voting.mod)
coef(summary(sat.voting.mod))
Use the states.rds data set. Fit a model predicting energy consumed per capita (energy) from the percentage of residents living in metropolitan areas (metro). Be sure to
summaryplot the model to look for deviations from modeling assumptionsSelect one or more additional predictors to add to your model and repeat steps 1-3. Is this model significantly better than the model with metro as the only predictor?
Interactions allow us assess the extent to which the association between one predictor and the outcome depends on a second predictor. For example: Does the association between expense and SAT scores depend on the median income in the state?
In [ ]:
#Add the interaction to the model
sat.expense.by.percent <- lm(csat ~ expense*income,
data=states.data)
#Show the results
coef(summary(sat.expense.by.percent)) # show regression coefficients table
In [ ]:
# make sure R knows region is categorical
str(states.data$region)
states.data$region <- factor(states.data$region)
#Add region to the model
sat.region <- lm(csat ~ region,
data=states.data)
#Show the results
coef(summary(sat.region)) # show regression coefficients table
anova(sat.region) # show ANOVA table
Again, make sure to tell R which variables are categorical by converting them to factors!
In the previous example we use the default contrasts for region. The default in R is treatment contrasts, with the first level as the reference. We can change the reference group or use another coding scheme using the C function.
In [ ]:
# print default contrasts
contrasts(states.data$region)
# change the reference group
coef(summary(lm(csat ~ C(region, base=4),
data=states.data)))
# change the coding scheme
coef(summary(lm(csat ~ C(region, contr.helmert),
data=states.data)))
See also ?contrasts, ?contr.treatment, and ?relevel.
Use the states data set.
Add on to the regression equation that you created in exercise 1 by generating an interaction term and testing the interaction.
Try adding region to the model. Are there significant differences across the four regions?
This far we have used the lm function to fit our regression models. lm is great, but limited--in particular it only fits models for continuous dependent variables. For categorical dependent variables we can use the glm() function.
For these models we will use a different dataset, drawn from the National Health Interview Survey. From the CDC website:
The National Health Interview Survey (NHIS) has monitored the health of the nation since 1957. NHIS data on a broad range of health topics are collected through personal household interviews. For over 50 years, the U.S. Census Bureau has been the data collection agent for the National Health Interview Survey. Survey results have been instrumental in providing data to track health status, health care access, and progress toward achieving national health objectives.
Load the National Health Interview Survey data:
In [ ]:
NH11 <- readRDS("dataSets/NatHealth2011.rds")
labs <- attributes(NH11)$labels
In [ ]:
str(NH11$hypev) # check stucture of hypev
levels(NH11$hypev) # check levels of hypev
# collapse all missing values to NA
NH11$hypev <- factor(NH11$hypev, levels=c("2 No", "1 Yes"))
# run our regression model
hyp.out <- glm(hypev~age_p+sex+sleep+bmi,
data=NH11, family="binomial")
coef(summary(hyp.out))
Generalized linear models use link functions, so raw coefficients are difficult to interpret. For example, the age coefficient of .06 in the previous model tells us that for every one unit increase in age, the log odds of hypertension diagnosis increases by 0.06. Since most of us are not used to thinking in log odds this is not too helpful!
One solution is to transform the coefficients to make them easier to interpret
In [ ]:
hyp.out.tab <- coef(summary(hyp.out))
hyp.out.tab[, "Estimate"] <- exp(coef(hyp.out))
hyp.out.tab
In addition to transforming the log-odds produced by glm to odds, we can use the predict() function to make direct statements about the predictors in our model. For example, we can ask "How much more likely is a 63 year old female to have hypertension compared to a 33 year old female?".
In [ ]:
# Create a dataset with predictors set at desired levels
predDat <- with(NH11,
expand.grid(age_p = c(33, 63),
sex = "2 Female",
bmi = mean(bmi, na.rm = TRUE),
sleep = mean(sleep, na.rm = TRUE)))
# predict hypertension at those levels
cbind(predDat, predict(hyp.out, type = "response",
se.fit = TRUE, interval="confidence",
newdata = predDat))
This tells us that a 33 year old female has a 13% probability of having been diagnosed with hypertension, while and 63 year old female has a 48% probability of having been diagnosed.
Instead of doing all this ourselves, we can use the effects package to compute quantities of interest for us (cf. the Zelig package).
In [ ]:
library(effects)
plot(allEffects(hyp.out))
Use the NH11 data set that we loaded earlier.
Note that the data is not perfectly clean and ready to be modeled. You will need to clean up at least some of the variables before fitting the model.
lmer function for liner mixed models, glmer for generalized mixed models
In [ ]:
library(lme4)
The Exam data set contans exam scores of 4,059 students from 65 schools in Inner London. The variable names are as follows:
| variable | Description |
|---|---|
| school | School ID - a factor. |
| normexam | Normalized exam score. |
| schgend | School gender - a factor. Levels are 'mixed', 'boys', and 'girls'. |
| schavg | School average of intake score. |
| vr | Student level Verbal Reasoning (VR) score band at intake - 'bottom 25%', 'mid 50%', and 'top 25%'. |
| intake | Band of student's intake score - a factor. Levels are 'bottom 25%', 'mid 50%' and 'top 25%'./ |
| standLRT | Standardised LR test score. |
| sex | Sex of the student - levels are 'F' and 'M'. |
| type | School type - levels are 'Mxd' and 'Sngl'. |
| student | Student id (within school) - a factor |
In [ ]:
Exam <- readRDS("dataSets/Exam.rds")
In [ ]:
# null model, grouping by school but not fixed effects.
Norm1 <-lmer(normexam ~ 1 + (1|school),
data=Exam, REML = FALSE)
summary(Norm1)
In [ ]:
Norm2 <-lmer(normexam~standLRT + (1|school),
data=Exam,
REML = FALSE)
summary(Norm2)
In [ ]:
anova(Norm1, Norm2)
In [ ]:
Norm3 <- lmer(normexam~standLRT + (standLRT|school), data=Exam,
REML = FALSE)
summary(Norm3)
In [ ]:
anova(Norm2, Norm3)
Use the dataset, bh1996: srcR[]{data(bh1996, package="multilevel")}
From the data documentation:
Variables are Cohesion (COHES), Leadership Climate (LEAD), Well-Being (WBEING) and Work Hours (HRS). Each of these variables has two variants - a group mean version that replicates each group mean for every individual, and a within-group version where the group mean is subtracted from each individual response. The group mean version is designated with a G. (e.g., G.HRS), and the within-group version is designated with a W. (e.g., W.HRS).
Use the states.rds data set.
In [ ]:
states <- readRDS("dataSets/states.rds")
Fit a model predicting energy consumed per capita (energy) from the percentage of residents living in metropolitan areas (metro). Be sure to
In [ ]:
states.en.met <- subset(states, select = c("metro", "energy"))
summary(states.en.met)
plot(states.en.met)
cor(states.en.met, use="pairwise")
summary
In [ ]:
mod.en.met <- lm(energy ~ metro, data = states)
summary(mod.en.met)
plot the model to look for deviations from modeling assumptions
In [ ]:
plot(mod.en.met)
Select one or more additional predictors to add to your model and repeat steps 1-3. Is this model significantly better than the model with metro as the only predictor?
In [ ]:
states.en.met.pop.wst <- subset(states, select = c("energy", "metro", "pop", "waste"))
summary(states.en.met.pop.wst)
plot(states.en.met.pop.wst)
cor(states.en.met.pop.wst, use = "pairwise")
mod.en.met.pop.waste <- lm(energy ~ metro + pop + waste, data = states)
summary(mod.en.met.pop.waste)
anova(mod.en.met, mod.en.met.pop.waste)
In [ ]:
mod.en.metro.by.waste <- lm(energy ~ metro * waste, data = states)
In [ ]:
mod.en.region <- lm(energy ~ metro * waste + region, data = states)
anova(mod.en.region)
Use the NH11 data set that we loaded earlier. Note that the data is not perfectly clean and ready to be modeled. You will need to clean up at least some of the variables before fitting the model.
In [ ]:
nh11.wrk.age.mar <- subset(NH11, select = c("everwrk", "age_p", "r_maritl"))
summary(nh11.wrk.age.mar)
NH11 <- transform(NH11,
everwrk = factor(everwrk,
levels = c("1 Yes", "2 No")),
r_maritl = droplevels(r_maritl))
mod.wk.age.mar <- glm(everwrk ~ age_p + r_maritl, data = NH11,
family = "binomial")
summary(mod.wk.age.mar)
In [ ]:
library(effects)
data.frame(Effect("r_maritl", mod.wk.age.mar))
In [ ]:
data(bh1996, package="multilevel")
From the data documentation:
Variables are Cohesion (COHES), Leadership Climate (LEAD), Well-Being (WBEING) and Work Hours (HRS). Each of these variables has two variants - a group mean version that replicates each group mean for every individual, and a within-group version where the group mean is subtracted from each individual response. The group mean version is designated with a G. (e.g., G.HRS), and the within-group version is designated with a W. (e.g., W.HRS).
Note that the group identifier is named "GRP".
In [ ]:
library(lme4)
mod.grp0 <- lmer(WBEING ~ 1 + (1 | GRP), data = bh1996)
summary(mod.grp0)
`> library(lme4)
mod.grp0 <- lmer(WBEING ~ 1 + (1 | GRP), data ` bh1996) > summary(mod.grp0) Linear mixed model fit by REML ['lmerMod'] Formula: WBEING ~ 1 + (1 | GRP) Data: bh1996
REML criterion at convergence: 19347
Scaled residuals: Min 1Q Median 3Q Max -3.322 -0.648 0.031 0.718 2.667
Random effects: Groups Name Variance Std.Dev. GRP (Intercept) 0.0358 0.189 Residual 0.7895 0.889 Number of obs: 7382, groups: GRP, 99
Fixed effects: Estimate Std. Error t value (Intercept) 2.7743 0.0222 125 > 2. [@2] Calculate the ICC for your null model
~ICC .0358/(.0358 + .7895) = .04~
In [ ]:
mod.grp1 <- lmer(WBEING ~ HRS + LEAD + (1 | GRP), data = bh1996)
summary(mod.grp1)
In [ ]:
mod.grp2 <- lmer(WBEING ~ HRS + LEAD + (1 + HRS | GRP), data = bh1996)
anova(mod.grp1, mod.grp2)
feedback form