-
Notifications
You must be signed in to change notification settings - Fork 0
/
Lab2.R
455 lines (312 loc) · 14 KB
/
Lab2.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
####################################################
## ##
## 204B Lab 2 ##
## VIOLATION OF ASSUMPTIONS ##
## ##
## Kristine O'Laughlin ##
## Winter 2018 ##
## ##
####################################################
## REVIEW OF HOMEWORK 2
####################################################
## 1. Download the hw_w2_1.RDat file from canvas
# and load it into R
setwd('C:/Users/Kristine/Dropbox/Winter2018/PSC204B/Homeworks')
load(file = 'hw_w2_1.RDat')
head(hw_w2_1)
# Run a linear regression and check with there is
# suppression, and if so, what type?
# Start with a correlation matrix to determine
# relations among our set of variables:
cor(hw_w2_1) # So we see that there are positive
# relations between the predictors
# and outcome variables, while
# there is a negative association
# between the two predictors.
# Model 1:
summary(lm(y ~ x1, data = hw_w2_1)) #r^2_y,x1 = .1306
# Model 2:
summary(lm(y ~ x2, data = hw_w2_1)) #r^2_y,x2 = .1987
# Model 3:
summary(lm(y ~ x1 + x2, data = hw_w2_1)) # R^2_y,x1x2 = .4151
# Is the sum of the r^2_y,x1 and r^2_y,x2 < R^2_y,x1x2?
(.1306 + .1987) < .4151 # TRUE--taken as evidence of
# suppression
# Because correlation between x1 and x2 is < 0, and
# relation between predictors with outcome is > 0,
# this is an example of cooperative suppression
## 2. Download the hw_w2_2.RDat file from canvas and
# load it into R
load(file = 'hw_w2_2.RDat')
head(hw_w2_2)
# Fit a quadratic model and center at the beginning
# of the study, at the grand-mean and at the end.
# Centering:
# Beginning:
hw_w2_2$x_b = hw_w2_2$x - min(hw_w2_2$x)
# Grand-mean:
hw_w2_2$x_m = hw_w2_2$x - mean(hw_w2_2$x)
# End:
hw_w2_2$x_e = hw_w2_2$x - max(hw_w2_2$x)
# What is the effect of centering the data?
par(mfrow = c(3, 1))
with(hw_w2_2, plot(x_b, y, xlim = c(-60, 60),
main = 'Centering at Beginning'))
with(hw_w2_2, plot(x_m, y, xlim = c(-60, 60),
main = 'Centering at Grand-Mean'))
with(hw_w2_2, plot(x_e, y, xlim = c(-60, 60),
main = 'Centering at End'))
# Now fitting the quadratic model for each:
# Beginning:
fit1 = summary(lm(y ~ x_b + I(x_b^2), data = hw_w2_2))
# Grand-mean:
fit2 = summary(lm(y ~ x_m + I(x_m^2), data = hw_w2_2))
# End:
fit3 = summary(lm(y ~ x_e + I(x_e^2), data = hw_w2_2))
coef(fit1); coef(fit2); coef(fit3)
# Plot the figure with the fitted line and the simple
# slope or instantaneous rate of change at the three
# centering points.
# Start with bivariate scatterplot of raw data:
par(mfrow = c(1, 1))
with(hw_w2_2, plot(x, y))
# Creating squared variable (this will make it
# easier to plot the quadratic function later):
hw_w2_2$x2 = hw_w2_2$x^2
# Fit the quadratic model to the raw data:
fit = lm(y ~ x + x2, data = hw_w2_2)
# Create a sequence of x values at short intervals
# to get a smooth function:
times = seq(min(hw_w2_2$x), max(hw_w2_2$x), .01)
# Get predicted y-values for each of the values
# above
pred = predict(fit, list(x = times, x2 = times^2))
# Add the quadratic function:
lines(times, pred, col = 'red', lwd = 3)
# To get the instantaneous rate of change, we can
# pull this value from each of the fitted models
# we ran earlier:
# Beginning:
slope_b = coef(fit1)[2, 1]
# Grand-mean:
slope_m = coef(fit2)[2, 1]
# End:
slope_e = coef(fit3)[2, 1]
# Get predicted y-values using raw equation:
y_b = coef(fit1)[1]
y_m = coef(fit2)[1]
y_e = coef(fit3)[1]
# Get intercept for tangent lines:
b_b = y_b - slope_b*min(hw_w2_2$x)
b_m = y_m - slope_m*mean(hw_w2_2$x)
b_e = y_e - slope_e*max(hw_w2_2$x)
# We need to plot the instantaneous rate of change at each
# of the following points:
# Beginning:
points(min(hw_w2_2$x), y_b, pch = 8, cex = 1.8, col = 'blue',
lwd = 2)
# Grand-mean:
points(mean(hw_w2_2$x), y_m, pch = 8, cex = 1.8, col = 'green',
lwd = 2)
# End:
points(max(hw_w2_2$x), y_e, pch = 8, cex = 1.8, col = 'purple',
lwd = 2)
# Adding the tangent lines:
abline(b_b, slope_b, col = 'blue', lty = 2, lwd = 2)
abline(b_m, slope_m, col = 'green', lty = 2, lwd = 2)
abline(b_e, slope_e, col = 'purple', lty = 2, lwd = 2)
# Adding a legend:
legend('topright', c('Beginning', 'Mean', 'End'),
col = c('blue', 'green', 'purple'),
lty = c(2, 2, 2), lwd = c(2, 2, 2))
## 3. In this exercise you will simulate two variables that
# are statistically independent of each other to see what
# happens when we run a regression of one on the other.
## a. First generate 1000 data points from a normal
# distribution with mean 0 and standard deviation 1
# by typing var1 <- rnorm(1000,0,1) in R. Generate
# another variable in the same way (call it var2).
# Run a regression of one variable on the other. Is
# the slope coefficient statistically significant?
set.seed(012318)
var1 <- rnorm(1000,0,1)
var2 <- rnorm(1000,0,1)
summary(lm(var2 ~ var1)) # Slope is not statistically
# significant
## b. Now run a simulation repeating this process 100
# times. This can be done using a loop. From each
# simulation, save the z-score (the estimated
# coefficient of var1 divided by its standard error).
# If the absolute value of the z-score exceeds 2,
# the estimate is statistically significant. Here is
# code to perform the simulation:
# (If you do not have the following package installed,
# the simulation will not run).
# install.packages('arm')
library(arm)
z.scores <- rep (NA, 100)
set.seed(012618)
for (k in 1:100) {
var1 <- rnorm (1000,0,1)
var2 <- rnorm (1000,0,1)
fit <- lm (var2 ~ var1)
z.scores[k] <- coef(fit)[2]/se.coef(fit)[2]
}
# How many of these 100 z-scores are statistically
# significant?
# To answer this, we can look at how many of the results
# had an absolute value greater than 2:
sum(abs(z.scores) > 2)
# Our sum should be around 5...why?
## 4. Download the child.iq.dta dataset. You have
# access to children's test scores at age 3, mother's
# education, and the mother's age at the time she
# gave birth for a sample of 400 children. The data
# are a Stata file which you can read into R by saving
# in your working directory and then typing the
# following:
#install.packages('foreign')
library (foreign)
iq.data <- read.dta ("child.iq.dta")
# a. Fit a regression of child test scores on mother's
# age, display the data and fitted model, check
# assumptions, and interpret the slope coefficient.
# When do you recommend mothers should give birth?
# What are you assuming in making these
# recommendations?
head(iq.data)
iq.fit = lm(ppvt ~ momage, data = iq.data)
summary(iq.fit)
# We can get some useful diagnostic plots by using
# the plot() wrapper around the object from lm():
par(mfrow = c(2, 2)) # First rearranging the ploting window
plot(iq.fit)
par(mfrow = c(1, 1))
with(iq.data, plot(momage, ppvt, xlab = 'Mother\'s \n Age',
ylab = 'Child\'s Test Score'))
abline(lm(ppvt ~ momage, data = iq.data), col = 'red',
lwd = 3)
# So what is the interpretation of the slope here?
# When would you recommend that women have children?
# What assumptions do this recommendation depend on?
# b. Repeat this for a regression that further includes
# mother's education, interpreting both slope
# coefficients in this model. Have your conclusions
# about the timing of birth changed?
summary(lm(ppvt ~ momage + educ_cat, data = iq.data))
# How would you interpret the slope for mother's age?
# How would you interpret the slope for mother's education?
# Have your conclusions changed?
## 5. Download data from Hamermesh and Parker (2005) on
# student evaluations of instructors' beauty and
# teaching quality for several courses at the
# University of Texas. The teaching evaluations were
# conducted at the end of the semester, and the beauty
# judgments were made later, by six students who had
# not attended the classes and were not aware of the
# course evaluations.
load(file = './TeachingRatings.rda')
# a. Run a regression using beauty to predict course
# evaluations (eval), controlling for various other
# inputs. Display the fitted model graphically, and
# explaining the meaning of each of the coefficients.
# Plot the residuals versus fitted values.
head(TeachingRatings)
fit_eval = lm(eval ~ minority + age + gender + beauty,
data = TeachingRatings)
summary(fit_eval)
# We can plot the predicted regression line while controlling
# for covariates by holding values of the covariates constant
# Start by setting a range of values for beauty:
beauty_range = seq(min(TeachingRatings$beauty),
max(TeachingRatings$beauty), .01)
# Get values predicted for the fit_eval model, for white
# males, who are average age in the sample:
pred_eval_m = predict(fit_eval, list(minority = rep('no', length(beauty_range)),
age = rep(mean(TeachingRatings$age), length(beauty_range)),
gender = rep('male', length(beauty_range)),
beauty = beauty_range))
# Same as above, but now repeating for white females
# who are average age in the sample:
pred_eval_f = predict(fit_eval, list(minority = rep('no', length(beauty_range)),
age = rep(mean(TeachingRatings$age), length(beauty_range)),
gender = rep('female', length(beauty_range)),
beauty = beauty_range))
par(mfrow = c(1, 1))
with(TeachingRatings, plot(beauty, eval,
xlab = 'Beauty Rating',
ylab = 'Teacher Evaluation'))
lines(beauty_range, pred_eval_m, col = 'blue',
lwd = 3)
lines(beauty_range, pred_eval_f, col = 'red',
lwd = 3)
# Another approach--we can use visreg from the visreg()
# package to control for other variables in the plot:
# install.package('visreg')
library(visreg)
visreg(fit_eval)
# As we saw earlier, we can get plots for our fitted model,
# including a plot of residuals verus fitted values:
plot(fit_eval, which = 1) # We can use which = 1 to get
# only the first plot
# b. Fit some other models, including beauty and also
# other input variables. Consider at least one model
# with interactions. For each model, state what the
# predictors are, and what the inputs are (see Section
# 2.1), and explain the meaning of each of its
# coefficients.
summary(lm(eval ~ students + tenure + division + beauty,
data = TeachingRatings))
summary(lm(eval ~ gender + credits + tenure + beauty,
data = TeachingRatings))
summary(lm(eval ~ age + gender + age*gender + beauty,
data = TeachingRatings))
# Focusing on the final model above, how would we interpret
# each of the coefficients?
## 6. The data wfw90.dat is from the Work, Family, and Well-
# Being Survey (Ross, 1990). Pull out the data on earnings,
# sex, height, and weight.
dat = read.table('wfw90.dat', header = T)
head(dat)
# a. Fit a linear regression model predicting earnings
# from height. What transformation should you perform
# in order to interpret the intercept from this model
# as average earnings for people with average height?
summary(lm(earnings ~ height, data = dat))
# Need to mean center height to make the intercept
# interpretable as the average earnings for someone
# of average height.
dat$height_c = dat$height - mean(dat$height)
summary(lm(earnings ~ height_c, data = dat))
# Now, we interpret the intercept as the expected value
# for someone of average height. We see that this value
# is 2,400.76.
# b. Fit some regression models with the goal of
# predicting earnings from some combination of sex,
# height, and weight. Be sure to try various
# transformations and interactions that might make
# sense. Choose your preferred model and justify.
summary(lm(earnings ~ sex + height + sex*height, data = dat))
summary(lm(earnings ~ height + weight + height*weight, data = dat))
dat$bmi = (dat$weight*.45)/((dat$height*.025)^2)
summary(lm(earnings ~ sex + bmi + sex*bmi, data = dat))
# c. Interpret all model coefficients.
# How would we interpret the coefficients from the last
# model we fit above?
## 7. An economist runs a regression examining the
# relations between the average price of cigarettes,
# P, and the quantity purchased, Q, across a large
# sample of counties in the United States, assuming
# the following functional form, log Q = alpha+beta log P.
# Suppose the estimate for beta is 0.3. Interpret this
# coefficient.
# When the outcome is log-transformed, we know that we
# would interpret the slope as the percent change in the
# outcome for a one-unit change in x.
# When both the criterion and predictor variables are
# log-transformed (as is true in our example above), the
# relationship can be interpretted as the percent change
# in y for a one-percent increase in x.
# So for our current example, an increase in the price of
# cigarettes by 1% is associated with a 0.3% increase in the
# quantity of cigarettes purchased.