-
Notifications
You must be signed in to change notification settings - Fork 0
/
Sec 2.2: sex-height-weight
107 lines (92 loc) · 5.88 KB
/
Sec 2.2: sex-height-weight
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
# Dataset from a public site: 199 adults (111 females) with 3 variables: sex, height, and weight
# download data
dat <- read.table("https://www.mv.helsinki.fi/home/mjxpirin/HDS_course/material/Davis_height_weight.txt",
as.is = T, header=T)
dat$sex <- ifelse(dat$sex=='F', -0.5, 0.5) # dummy-coding sex
# 7 models in section 2.2
m1 <- lm(weight ~ sex, data = dat)
m2 <- lm(weight ~ height, data = dat)
m3 <- lm(weight ~ height + sex, data = dat)
m4 <- lm(weight ~ height * sex, data = dat)
m7 <- lm(height ~ weight * sex, data = dat)
m8 <- lm(height ~ sex, data = dat)
m9 <- lm(height ~ weight + sex, data = dat)
# m1 - sex difference in weight: point estimate & 95% interval
c(summary(m1)$coefficients['sex','Estimate'], summary(m1)$coefficients['sex','Estimate']-
qt(0.025, summary(m1)$df[2], lower.tail = F)*summary(m1)$coefficients['sex','Std. Error'],
summary(m1)$coefficients['sex','Estimate']+qt(0.025, summary(m1)$df[2], lower.tail = F)*
summary(m1)$coefficients['sex','Std. Error'])
#[1] 19.00584 16.35113 21.66054
# m3 - sex difference in weight: point estimate & 95% interval
c(summary(m3)$coefficients['sex','Estimate'], summary(m3)$coefficients['sex','Estimate']-
qt(0.025, summary(m3)$df[2], lower.tail = F)*summary(m3)$coefficients['sex','Std. Error'],
summary(m3)$coefficients['sex','Estimate']+qt(0.025, summary(m3)$df[2], lower.tail = F)*
summary(m3)$coefficients['sex','Std. Error'])
#[1] 8.216212 4.829524 11.602900
# m4 - sex difference in weight: point estimate & 95% interval
c(summary(m4)$coefficients['sex','Estimate'], summary(m4)$coefficients['sex','Estimate']-
qt(0.025, summary(m4)$df[2], lower.tail = F)*summary(m4)$coefficients['sex','Std. Error'],
summary(m4)$coefficients['sex','Estimate']+qt(0.025, summary(m4)$df[2], lower.tail = F)*
summary(m4)$coefficients['sex','Std. Error'])
#[1] -55.621635 -119.806446 8.563177
# m2 - height-weight relationship: point estimate & 95% interval
c(summary(m2)$coefficients['height', 'Estimate'], summary(m2)$coefficients['height', 'Estimate']-
qt(0.025, summary(m2)$df[2], lower.tail = F)*summary(m2)$coefficients['height','Std. Error'],
summary(m2)$coefficients['height', 'Estimate']+qt(0.025, summary(m2)$df[2], lower.tail = F)*
summary(m2)$coefficients['height','Std. Error'])
#[1] 1.149222 1.015734 1.282710
# m3 - height-weight relationship: point estimate & 95% interval
c(summary(m3)$coefficients['height', 'Estimate'], summary(m3)$coefficients['height', 'Estimate']-
qt(0.025, summary(m3)$df[2], lower.tail = F)*summary(m3)$coefficients['height','Std. Error'],
summary(m3)$coefficients['height', 'Estimate']+qt(0.025, summary(m3)$df[2], lower.tail = F)*
summary(m3)$coefficients['height','Std. Error'])
#[1] 0.8107219 0.6222912 0.9991527
# m4 - height-weight relationship: point estimate & 95% interval
c(summary(m4)$coefficients['height', 'Estimate'], summary(m4)$coefficients['height', 'Estimate']-
qt(0.025, summary(m4)$df[2], lower.tail = F)*summary(m4)$coefficients['height','Std. Error'],
summary(m4)$coefficients['height', 'Estimate']+qt(0.025, summary(m4)$df[2], lower.tail = F)*
summary(m4)$coefficients['height','Std. Error'])
#[1] 0.8092703 0.6221870 0.9963536
# m7 - sex difference in height: point estimate & 95% interval
c(summary(m7)$coefficients['sex','Estimate'], summary(m7)$coefficients['sex','Estimate']-
qt(0.025, summary(m7)$df[2], lower.tail = F)*summary(m7)$coefficients['sex','Std. Error'],
summary(m7)$coefficients['sex','Estimate']+qt(0.025, summary(m7)$df[2], lower.tail = F)*
summary(m7)$coefficients['sex','Std. Error'])
#[1] 15.248738 4.552696 25.944780
# m8 - sex difference in height: point estimate & 95% interval
c(summary(m8)$coefficients['sex','Estimate'], summary(m8)$coefficients['sex','Estimate']-
qt(0.025, summary(m8)$df[2], lower.tail = F)*summary(m8)$coefficients['sex','Std. Error'],
summary(m8)$coefficients['sex','Estimate']+qt(0.025, summary(m8)$df[2], lower.tail = F)*
summary(m8)$coefficients['sex','Std. Error'])
#[1] 13.30866 11.61144 15.00588
# m9 - sex difference in height: point estimate & 95% interval
c(summary(m9)$coefficients['sex','Estimate'], summary(m9)$coefficients['sex','Estimate']-
qt(0.025, summary(m9)$df[2], lower.tail = F)*summary(m9)$coefficients['sex','Std. Error'],
summary(m9)$coefficients['sex','Estimate']+qt(0.025, summary(m9)$df[2], lower.tail = F)*
summary(m9)$coefficients['sex','Std. Error'])
#[1] 7.010686 4.946643 9.074729
##### model in Section 4.3 - centered height
dd <- transform(dat, heightC=ave(height, sex,
FUN=function(x) scale(x, center = TRUE, scale = FALSE)))
# model with interaction
m5 <- lm(weight ~ heightC * sex, data = dd)
c(summary(m5)$coefficients['sex','Estimate'], summary(m5)$coefficients['sex','Estimate']-
qt(0.025, summary(m5)$df[2], lower.tail = F)*summary(m5)$coefficients['sex','Std. Error'],
summary(m5)$coefficients['sex','Estimate']+qt(0.025, summary(m5)$df[2], lower.tail = F)*
summary(m5)$coefficients['sex','Std. Error'])
#[1] 19.00584 16.74605 21.26563
# model without interaction
m6 <- lm(weight ~ heightC + sex, data = dd)
c(summary(m6)$coefficients['sex','Estimate'], summary(m6)$coefficients['sex','Estimate']-
qt(0.025, summary(m6)$df[2], lower.tail = F)*summary(m6)$coefficients['sex','Std. Error'],
summary(m6)$coefficients['sex','Estimate']+qt(0.025, summary(m6)$df[2], lower.tail = F)*
summary(m6)$coefficients['sex','Std. Error'])
#[1] 19.00584 16.72970 21.28197
##### estimate sex difference in weight with height of 190 cm
d190 <- dat; d190$height <- d190$height-190
m10 <- lm(weight ~ height * sex, data = d190)
c(summary(m10)$coefficients['sex','Estimate'], summary(m10)$coefficients['sex','Estimate']-
qt(0.025, summary(m10)$df[2], lower.tail = F)*summary(m10)$coefficients['sex','Std. Error'],
summary(m10)$coefficients['sex','Estimate']+qt(0.025, summary(m10)$df[2], lower.tail = F)*
summary(m10)$coefficients['sex','Std. Error'])
#[1] 15.182938 7.469974 22.895903