c:ms:2026:schedule:w10.lecture.note
Table of Contents
Correlation and Regression
- deriviation of a and b in a simple regression regression 에서 a 와 b 구하기
rs01
rm(list=ls())
rnorm2 <- function(n,mean,sd) {
mean + sd * scale(rnorm(n))
}
ss <- function(x) {
sum((x-mean(x))^2)
}
sp <- function(x, y) {
sum((x-mean(x))*(y-mean(y)))
}
set.seed(1)
n <- 40
x <- rnorm(n, 40, 4)
a <- 10
b <- 5
y <- a + rnorm(n, 1, 35) + b * x
dat <- data.frame(x, y)
# library(ggplot2)
# ggplot(datb, aes(x = xb, y = yb)) +
# geom_point() + # Adds scatter points
# geom_smooth(method = "lm", se = TRUE)
cor(dat$x, dat$y)
cor.test(dat$x, dat$y)
m.lma <- lm(y ~ x, data=dat)
summary(m.lma)
plot(dat)
abline(m.lma)
text(x=mean(x)+2*sd(x), y=mean(y)+2*sd(y), pos=1, col='black', labels = paste("r =", round(cor(x, y),5)))
abline(v=mean(x), col='blue')
abline(h=mean(y), col='red')
text(x=mean(x), y=mean(y)-sd(y), pos=4, col='blue', labels = paste("mean(x)\n=", round(mean(x),5)))
text(x=mean(x)-sd(x), y=mean(y)+sd(y), pos=4, col='red', labels = paste("mean(y)\n=", round(mean(y),5)))
cov(dat)==var(dat)
cov(dat$x, dat$y)
cor(dat$x, dat$y)
cor.test(dat$x, dat$y)
sp.xy <- sp(x, y)
sp.xy
df.xy <- length(x)-1
sp.xy/df.xy
cov(dat)
ss.x <- ss(x)
ss.y <- ss(y)
r.cal0 <- cov(x,y)/(sd(x)*sd(y))
r.cal <- sp(x,y)/sqrt(ss.x*ss.y)
r.cal
r.cal0
# commres.net
# what is b and a?
b.cal <- sp.xy / ss.x
a.cal <- mean(y) - b.cal * mean(x)
a.cal
b.cal
m.lma
summary(m.lma)
y.hat <- a.cal + b.cal * x
y
y.hat
res <- y - y.hat
ss.res <- sum(res^2)
ss.res
y.hat2 <- predict(m.lma, newdata=dat)
data.frame(y.hat, y.hat2)
y.m <- mean(y)
y.m
sum((y - y.m)^2)
ss(y)
res <- y - y.hat
res
m.lma$residuals
reg <- y.hat - y.m
ss.reg <- sum(reg^2)
ss.res <- sum(res^2)
ss.tot <- ss(y)
ss.tot
ss.res
ss.reg
ss.reg+ss.res
df.tot <- length(dat$y)-1
df.reg <- 2 - 1
df.res <- df.tot - df.reg
ms.tot <- ss.tot/df.tot
var(y)
ms.tot
ms.reg <- ss.reg / df.reg
ms.res <- ss.res / df.res
ms.reg
ms.res
f.cal <- ms.reg / ms.res
p.f <- pf(f.cal, df.reg, df.res, lower.tail = F)
cat(f.cal, df.reg, df.res, p.f)
r.sq <- ss.reg / ss.tot
cat("R-squared: ", r.sq)
cat("a: ", a.cal, "b: ", b.cal)
cat(min(res), max(res))
var.res <- ss.res / (n-2) # variance of residuals
sd.res <- sqrt(var.res) # standard deviation of residuals
sd.res # residual standard error
sigma(m.lma)
# residual standard error
summary(m.lma)
# 한편
# standard error for b
r.square <- summary(m.lma)$r.square
r.square
(r.cal)^2
ss.reg / ss.tot
1-r.square
ss.res/ss.tot
1-(ss.reg/ss.tot)
# sd(res)?
sd.res
var.res <- ss.res/(n-2)
var.res
sqrt(var.res)
ss.res # ss.res 값
ss.res/ss.tot # ss.res값을 ss.tot 비율로 본 것
sd.res <- ss.res/(n-2)
# 위의 ss.res 대신에
# 아래처럼 ss.res/ss.tot 을 쓴다면
# ss.res의 proportion값을 (전체 ss 값에 대한 비율값을)
# 사용한 것이 된다
sqrt((ss.res/ss.tot)/(n-2))
# 그리고 위는 아래와 같은 것
sqrt((1-r.square)/(n-2))
se.res <- sqrt((ss.res/ss.tot)/(n-2)) # se for residual
se.res
t.r <- r.cal/se.res
p.r <- pt(t.r, n-2, lower.tail = F)*2
cat(t.r, p.r)
cor.test(x,y)
# standard error for b1 (regression slope)
sd.res <- sqrt(ss.res/(n-2)) # residual sd (sd for regression)
sd.res
ss.x
se.b <- sd.res/sqrt(ss.x)
se.b
b.cal
t.b <- b.cal/se.b
p.b <- pt(t.b, n-2, lower.tail = F)*2
data.frame(b.cal, se.b, t.b, p.b)
summary(lm(y~x))
data.frame(ss.reg, ss.res, ss.tot, r.square, ss.res/ss.tot)
data.frame(f.cal, df.reg, df.res, p.f)
var.res <- ss.res/(n-2)
sd.res <- sqrt(var.res)
sd.res
ro01
> rm(list=ls())
>
> rnorm2 <- function(n,mean,sd) {
+ mean + sd * scale(rnorm(n))
+ }
> ss <- function(x) {
+ sum((x-mean(x))^2)
+ }
> sp <- function(x, y) {
+ sum((x-mean(x))*(y-mean(y)))
+ }
>
> set.seed(1)
> n <- 40
> x <- rnorm(n, 40, 4)
> a <- 10
> b <- 5
> y <- a + rnorm(n, 1, 35) + b * x
> dat <- data.frame(x, y)
>
> # library(ggplot2)
> # ggplot(datb, aes(x = xb, y = yb)) +
> # geom_point() + # Adds scatter points
> # geom_smooth(method = "lm", se = TRUE)
> cor(dat$x, dat$y)
[1] 0.6237667
> cor.test(dat$x, dat$y)
Pearson's product-moment correlation
data: dat$x and dat$y
t = 4.9195, df = 38, p-value = 1.707e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.3875647 0.7831105
sample estimates:
cor
0.6237667
> m.lma <- lm(y ~ x, data=dat)
> summary(m.lma)
Call:
lm(formula = y ~ x, data = dat)
Residuals:
Min 1Q Median 3Q Max
-65.304 -22.086 -2.542 14.342 72.912
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -69.434 58.453 -1.188 0.242
x 7.097 1.443 4.920 1.71e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 31.95 on 38 degrees of freedom
Multiple R-squared: 0.3891, Adjusted R-squared: 0.373
F-statistic: 24.2 on 1 and 38 DF, p-value: 1.707e-05
>
> plot(dat)
> abline(m.lma)
> text(x=mean(x)+2*sd(x), y=mean(y)+2*sd(y), pos=1, col='black', labels = paste("r =", round(cor(x, y),5)))
> abline(v=mean(x), col='blue')
> abline(h=mean(y), col='red')
> text(x=mean(x), y=mean(y)-sd(y), pos=4, col='blue', labels = paste("mean(x)\n=", round(mean(x),5)))
> text(x=mean(x)-sd(x), y=mean(y)+sd(y), pos=4, col='red', labels = paste("mean(y)\n=", round(mean(y),5)))
>
> cov(dat)==var(dat)
x y
x TRUE TRUE
y TRUE TRUE
> cov(dat$x, dat$y)
[1] 89.26937
> cor(dat$x, dat$y)
[1] 0.6237667
> cor.test(dat$x, dat$y)
Pearson's product-moment correlation
data: dat$x and dat$y
t = 4.9195, df = 38, p-value = 1.707e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.3875647 0.7831105
sample estimates:
cor
0.6237667
>
> sp.xy <- sp(x, y)
> sp.xy
[1] 3481.506
> df.xy <- length(x)-1
> sp.xy/df.xy
[1] 89.26937
> cov(dat)
x y
x 12.57885 89.26937
y 89.26937 1628.24418
> ss.x <- ss(x)
> ss.y <- ss(y)
>
> r.cal0 <- cov(x,y)/(sd(x)*sd(y))
> r.cal <- sp(x,y)/sqrt(ss.x*ss.y)
> r.cal
[1] 0.6237667
> r.cal0
[1] 0.6237667
>
> # commres.net
> # what is b and a?
> b.cal <- sp.xy / ss.x
> a.cal <- mean(y) - b.cal * mean(x)
> a.cal
[1] -69.43374
> b.cal
[1] 7.096781
> m.lma
Call:
lm(formula = y ~ x, data = dat)
Coefficients:
(Intercept) x
-69.434 7.097
> summary(m.lma)
Call:
lm(formula = y ~ x, data = dat)
Residuals:
Min 1Q Median 3Q Max
-65.304 -22.086 -2.542 14.342 72.912
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -69.434 58.453 -1.188 0.242
x 7.097 1.443 4.920 1.71e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 31.95 on 38 degrees of freedom
Multiple R-squared: 0.3891, Adjusted R-squared: 0.373
F-statistic: 24.2 on 1 and 38 DF, p-value: 1.707e-05
>
> y.hat <- a.cal + b.cal * x
> y
[1] 192.7126 205.8052 218.6811 262.3888 193.4837 169.8283 233.5089 252.6651 218.5835 235.7310 255.1693
[12] 197.3759 210.5144 127.1783 283.6544 279.4153 197.8234 193.3320 247.3646 218.1511 313.4362 225.2693
[23] 236.6322 172.1930 197.3820 216.4852 144.7105 232.8794 206.8009 295.4002 254.8164 184.0961 240.1289
[34] 177.2305 139.5816 212.9007 187.5990 209.8524 235.6025 205.6303
> y.hat
[1] 196.6543 219.6506 190.7164 259.7229 223.7913 191.1468 228.2742 235.3964 230.7823 205.7684 257.3526
[12] 225.5040 196.8023 151.5685 246.3711 213.1620 213.9779 241.2303 237.7496 231.2967 240.5246 236.6401
[23] 216.5542 157.9655 232.0326 212.8442 210.0149 172.6871 200.8642 226.3017 253.0065 211.5197 225.4424
[34] 212.9101 175.3467 202.6570 203.2447 212.7538 245.6641 236.1019
> res <- y - y.hat
> ss.res <- sum(res^2)
> ss.res
[1] 38794.04
>
> y.hat2 <- predict(m.lma, newdata=dat)
> data.frame(y.hat, y.hat2)
y.hat y.hat2
1 196.6543 196.6543
2 219.6506 219.6506
3 190.7164 190.7164
4 259.7229 259.7229
5 223.7913 223.7913
6 191.1468 191.1468
7 228.2742 228.2742
8 235.3964 235.3964
9 230.7823 230.7823
10 205.7684 205.7684
11 257.3526 257.3526
12 225.5040 225.5040
13 196.8023 196.8023
14 151.5685 151.5685
15 246.3711 246.3711
16 213.1620 213.1620
17 213.9779 213.9779
18 241.2303 241.2303
19 237.7496 237.7496
20 231.2967 231.2967
21 240.5246 240.5246
22 236.6401 236.6401
23 216.5542 216.5542
24 157.9655 157.9655
25 232.0326 232.0326
26 212.8442 212.8442
27 210.0149 210.0149
28 172.6871 172.6871
29 200.8642 200.8642
30 226.3017 226.3017
31 253.0065 253.0065
32 211.5197 211.5197
33 225.4424 225.4424
34 212.9101 212.9101
35 175.3467 175.3467
36 202.6570 202.6570
37 203.2447 203.2447
38 212.7538 212.7538
39 245.6641 245.6641
40 236.1019 236.1019
>
> y.m <- mean(y)
> y.m
[1] 217.0499
> sum((y - y.m)^2)
[1] 63501.52
> ss(y)
[1] 63501.52
>
> res <- y - y.hat
> res
[1] -3.941684 -13.845403 27.964735 2.665889 -30.307576 -21.318465 5.234736 17.268727 -12.198772
[10] 29.962596 -2.183295 -28.128092 13.712107 -24.390249 37.283390 66.253356 -16.154466 -47.898288
[19] 9.614998 -13.145540 72.911540 -11.370779 20.077988 14.227511 -34.650622 3.640985 -65.304380
[28] 60.192299 5.936666 69.098576 1.809915 -27.423536 14.686468 -35.679652 -35.765104 10.243725
[37] -15.645761 -2.901348 -10.061608 -30.471587
> m.lma$residuals
1 2 3 4 5 6 7 8 9
-3.941684 -13.845403 27.964735 2.665889 -30.307576 -21.318465 5.234736 17.268727 -12.198772
10 11 12 13 14 15 16 17 18
29.962596 -2.183295 -28.128092 13.712107 -24.390249 37.283390 66.253356 -16.154466 -47.898288
19 20 21 22 23 24 25 26 27
9.614998 -13.145540 72.911540 -11.370779 20.077988 14.227511 -34.650622 3.640985 -65.304380
28 29 30 31 32 33 34 35 36
60.192299 5.936666 69.098576 1.809915 -27.423536 14.686468 -35.679652 -35.765104 10.243725
37 38 39 40
-15.645761 -2.901348 -10.061608 -30.471587
>
> reg <- y.hat - y.m
> ss.reg <- sum(reg^2)
> ss.res <- sum(res^2)
> ss.tot <- ss(y)
> ss.tot
[1] 63501.52
> ss.res
[1] 38794.04
> ss.reg
[1] 24707.48
> ss.reg+ss.res
[1] 63501.52
>
> df.tot <- length(dat$y)-1
> df.reg <- 2 - 1
> df.res <- df.tot - df.reg
>
> ms.tot <- ss.tot/df.tot
> var(y)
[1] 1628.244
> ms.tot
[1] 1628.244
>
> ms.reg <- ss.reg / df.reg
> ms.res <- ss.res / df.res
> ms.reg
[1] 24707.48
> ms.res
[1] 1020.896
>
>
> f.cal <- ms.reg / ms.res
> p.f <- pf(f.cal, df.reg, df.res, lower.tail = F)
> cat(f.cal, df.reg, df.res, p.f)
24.20177 1 38 1.706524e-05>
> r.sq <- ss.reg / ss.tot
> cat("R-squared: ", r.sq)
R-squared: 0.3890849>
> cat("a: ", a.cal, "b: ", b.cal)
a: -69.43374 b: 7.096781> cat(min(res), max(res))
-65.30438 72.91154>
> var.res <- ss.res / (n-2) # variance of residuals
> sd.res <- sqrt(var.res) # standard deviation of residuals
> sd.res # residual standard error
[1] 31.95146
> sigma(m.lma)
[1] 31.95146
>
> # residual standard error
> summary(m.lma)
Call:
lm(formula = y ~ x, data = dat)
Residuals:
Min 1Q Median 3Q Max
-65.304 -22.086 -2.542 14.342 72.912
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -69.434 58.453 -1.188 0.242
x 7.097 1.443 4.920 1.71e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 31.95 on 38 degrees of freedom
Multiple R-squared: 0.3891, Adjusted R-squared: 0.373
F-statistic: 24.2 on 1 and 38 DF, p-value: 1.707e-05
>
> # 한편
> # standard error for b
> r.square <- summary(m.lma)$r.square
> r.square
[1] 0.3890849
> (r.cal)^2
[1] 0.3890849
> ss.reg / ss.tot
[1] 0.3890849
>
> 1-r.square
[1] 0.6109151
> ss.res/ss.tot
[1] 0.6109151
> 1-(ss.reg/ss.tot)
[1] 0.6109151
>
> # sd(res)?
> sd.res
[1] 31.95146
> var.res <- ss.res/(n-2)
> var.res
[1] 1020.896
> sqrt(var.res)
[1] 31.95146
> ss.res # ss.res 값
[1] 38794.04
> ss.res/ss.tot # ss.res값을 ss.tot 비율로 본 것
[1] 0.6109151
>
> sd.res <- ss.res/(n-2)
> # 위의 ss.res 대신에
> # 아래처럼 ss.res/ss.tot 을 쓴다면
> # ss.res의 proportion값을 (전체 ss 값에 대한 비율값을)
> # 사용한 것이 된다
> sqrt((ss.res/ss.tot)/(n-2))
[1] 0.126794
> # 그리고 위는 아래와 같은 것
> sqrt((1-r.square)/(n-2))
[1] 0.126794
>
> se.res <- sqrt((ss.res/ss.tot)/(n-2)) # se for residual
> se.res
[1] 0.126794
>
> t.r <- r.cal/se.res
> p.r <- pt(t.r, n-2, lower.tail = F)*2
> cat(t.r, p.r)
4.919529 1.706524e-05> cor.test(x,y)
Pearson's product-moment correlation
data: x and y
t = 4.9195, df = 38, p-value = 1.707e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.3875647 0.7831105
sample estimates:
cor
0.6237667
>
> # standard error for b1 (regression slope)
> sd.res <- sqrt(ss.res/(n-2)) # residual sd (sd for regression)
> sd.res
[1] 31.95146
> ss.x
[1] 490.5753
> se.b <- sd.res/sqrt(ss.x)
> se.b
[1] 1.442573
> b.cal
[1] 7.096781
> t.b <- b.cal/se.b
> p.b <- pt(t.b, n-2, lower.tail = F)*2
> data.frame(b.cal, se.b, t.b, p.b)
b.cal se.b t.b p.b
1 7.096781 1.442573 4.919529 1.706524e-05
>
> summary(lm(y~x))
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-65.304 -22.086 -2.542 14.342 72.912
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -69.434 58.453 -1.188 0.242
x 7.097 1.443 4.920 1.71e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 31.95 on 38 degrees of freedom
Multiple R-squared: 0.3891, Adjusted R-squared: 0.373
F-statistic: 24.2 on 1 and 38 DF, p-value: 1.707e-05
>
> data.frame(ss.reg, ss.res, ss.tot, r.square, ss.res/ss.tot)
ss.reg ss.res ss.tot r.square ss.res.ss.tot
1 24707.48 38794.04 63501.52 0.3890849 0.6109151
> data.frame(f.cal, df.reg, df.res, p.f)
f.cal df.reg df.res p.f
1 24.20177 1 38 1.706524e-05
> var.res <- ss.res/(n-2)
> sd.res <- sqrt(var.res)
> sd.res
[1] 31.95146
>
Gradient descent
- deriviation of a and b in a simple regression regression 에서 a 와 b 구하기
배운 것 확인
rs.02
############## # 배운 것 확인 ############## head(sam) tail(sam) mo.xy <- lm(y~x, data=sam) summary(mo.xy) plot(x,y) abline(mo.xy) abline(v=mean(x),col="red") abline(h=mean(y),col="blue") n <- length(y) n cov(x,y) sp(x,y)/(n-1) r.cal <- cor(x,y) r.cal cov(x,y)/(sd(x)*sd(y)) b.cal <- sp(x,y)/ss(x) a.cal <- mean(y)-b.cal*mean(x) cat(a.cal, b.cal) ss.tot <- ss(y) y.hat <- a.cal + b.cal*x mo.xy$fitted.values y.hat data.frame(mo.xy$fitted.values, y.hat) res <- y - y.hat reg <- y.hat - mean(y) res ss.res <- sum(res^2) ss.reg <- sum(reg^2) ss.tot ss.res ss.reg ss.res+ss.reg r.sq <- ss.reg/ss.tot r.sq 1-r.sq ss.res/ss.tot df.tot <- n - 1 df.reg <- 2 - 1 df.res <- df.tot - df.reg df.tot df.reg df.res ms.tot <- ss.tot / df.tot ms.tot var(y) ms.reg <- ss.reg / df.reg ms.res <- ss.res / df.res f.cal <- ms.reg / ms.res p.f <- pf(f.cal, df.reg, df.res, lower.tail = F) data.frame(f.cal, df.reg, df.res, p.f) summary(mo.xy) data.frame(head(res), head(summary(mo.xy)$residuals)) data.frame(min(res), max(res)) # residual의 분산값은 무엇인가? # residual의 집합인 res 변인의 ss(res) 값을 # ss.res 라고 하면 이를 n-2로 나눠준 값을 말한다. var.res <- ss.res / (n-2) # variance of residuals sd.res <- sqrt(var.res) # standard deviation of residuals # 위의 값을 standard deviation of residual이라고 부를 수 있다 data.frame(var.res, sd.res) sigma(mo.xy) summary(mo.xy)$sigma summary(mo.xy) r.cal r.sq r.cal^2 ss.reg / ss.tot 1-r.sq ss.res/ss.tot 1-(ss.reg/ss.tot) # 바로 위에서 한 residual 집합의 분산과 표준편차값 ss.res # ss.res 값 var.res <- ss.res/(n-2) var.res sqrt(var.res) sd.res ms.res <- var.res # 위의 var.res = ms.res 이라고도 부른다 # 그리고 이것은 ss.tot - ss.reg 방법 말고도 summary(mo.xy)$residuals res.1 <- summary(mo.xy)$residuals # 이것을 우리는 이미 위에서 # res <- y - y.hat 으로도 구했다. y.hat mo.xy$fitted.values data.frame(y.hat, mo.xy$fitted.values) sum((res.1-mean(res.1))^2) ss.res # 그런데 위의 식에서 ss.res 대신에 아래처럼 # ss.res 이 ss.total에서 차지하는 비율로 # 보면 ss.res/ss.tot 1-(ss.reg/ss.tot) 1-r.sq sd.res <- sqrt(ss.res/(n-2)) # 위의 ss.res 대신에 # 아래처럼 ss.res/ss.tot 을 쓴다면 # ss.res의 proportion값을 (전체 ss 값에 대한 비율값을) # 사용한 것이 된다. 즉, 일종의 표준화된 퍼센티지로 # 계산을 하는 것이 된다. sqrt((ss.res/ss.tot)/(n-2)) # 그리고 위는 아래와 같은 것 sqrt((1-r.sq)/(n-2)) # 이것을 se of residual 이라고 부른다 se.res <- sqrt((ss.res/ss.tot)/(n-2)) # se for residual se.res # 그리고 이 se.res은 r 값의 (correlation coefficient) # significance 정도를 가늠하는데 쓰인다. # 즉, r / se.res 이는 # t distribution을 따르기에 t.r <- r.cal/se.res p.r <- pt(t.r, n-2, lower.tail = F)*2 cat(r.cal, se.res, t.r, p.r) cor.test(x,y) # 그렇다면 b1의 significance한 정도는 어떻게 가늠하는가? # standard error for b1 (regression slope) # 위에서 우리는 ms.res <- ss.res/(n-2) ms.res sqrt(ms.res) sd.res <- sqrt(ms.res) # residual sd (sd for regression) sd.res ss.x <- ss(x) ss.x se.b <- sd.res/sqrt(ss.x) se.b sqrt(ms.res)/sqrt(ss.x) # 위의 값이 standard error of b 값이므로 .. .. # http://commres.net/regression#standard_error_of_b b.cal t.b <- b.cal/se.b p.b <- pt(t.b, n-2, lower.tail = F)*2 data.frame(b.cal, se.b, t.b, p.b) data.frame(r.cal, se.res, t.r, p.r ) # 왜 t.b와 t.r 값이 같은가? # summary(mo.xy) data.frame(ss.reg, ss.res, ss.tot, r.sq, ss.res/ss.tot) data.frame(f.cal, df.reg, df.res, p.f) anova(mo.xy) # remember ss.res? ms.res ms.reg
ro.02
> ##############
> # 배운 것 확인
> ##############
> head(sam)
x x2 y
1 100.7858 41.32553 3.192923
2 122.7128 51.70909 3.606455
3 115.7513 27.59380 3.299693
4 115.4523 36.08332 2.811801
5 113.1708 42.87865 3.247758
6 111.3423 54.87082 3.329180
> tail(sam)
x x2 y
95 110.09457 41.37388 3.252026
96 122.42522 35.62619 3.131744
97 119.04052 51.68809 3.281162
98 106.59454 40.61295 3.122031
99 97.07284 51.14235 2.578662
100 104.02543 41.19795 3.539551
> mo.xy <- lm(y~x, data=sam)
> summary(mo.xy)
Call:
lm(formula = y ~ x, data = sam)
Residuals:
Min 1Q Median 3Q Max
-0.92770 -0.20715 0.03368 0.20038 0.90512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.614428 0.344622 4.685 9.04e-06 ***
x 0.014476 0.003212 4.508 1.82e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3346 on 98 degrees of freedom
Multiple R-squared: 0.1717, Adjusted R-squared: 0.1633
F-statistic: 20.32 on 1 and 98 DF, p-value: 1.818e-05
> plot(x,y)
> abline(mo.xy)
> abline(v=mean(x),col="red")
> abline(h=mean(y),col="blue")
>
> n <- length(y)
> n
[1] 100
> cov(x,y)
[1] 1.587297
> sp(x,y)/(n-1)
[1] 1.587297
> r.cal <- cor(x,y)
> r.cal
[1] 0.414393
> cov(x,y)/(sd(x)*sd(y))
[1] 0.414393
>
> b.cal <- sp(x,y)/ss(x)
> a.cal <- mean(y)-b.cal*mean(x)
> cat(a.cal, b.cal)
1.614428 0.01447618>
> ss.tot <- ss(y)
> y.hat <- a.cal + b.cal*x
> mo.xy$fitted.values
. . . .
> y.hat
. . . .
> data.frame(mo.xy$fitted.values, y.hat)
mo.xy.fitted.values y.hat
1 3.073421 3.073421
2 3.390841 3.390841
3 3.290064 3.290064
4 3.285737 3.285737
5 3.252709 3.252709
6 3.226239 3.226239
> res <- y - y.hat
> reg <- y.hat - mean(y)
> res
. . . .
> ss.res <- sum(res^2)
> ss.reg <- sum(reg^2)
> ss.tot
[1] 13.24715
> ss.res
[1] 10.97233
> ss.reg
[1] 2.274822
> ss.res+ss.reg
[1] 13.24715
>
> r.sq <- ss.reg/ss.tot
> r.sq
[1] 0.1717215
> 1-r.sq
[1] 0.8282785
> ss.res/ss.tot
[1] 0.8282785
>
> df.tot <- n - 1
> df.reg <- 2 - 1
> df.res <- df.tot - df.reg
> df.tot
[1] 99
> df.reg
[1] 1
> df.res
[1] 98
>
> ms.tot <- ss.tot / df.tot
> ms.tot
[1] 0.1338096
> var(y)
[1] 0.1338096
> ms.reg <- ss.reg / df.reg
> ms.res <- ss.res / df.res
> f.cal <- ms.reg / ms.res
> p.f <- pf(f.cal, df.reg, df.res, lower.tail = F)
> data.frame(f.cal, df.reg, df.res, p.f)
f.cal df.reg df.res p.f
1 20.3177 1 98 1.817763e-05
>
> summary(mo.xy)
Call:
lm(formula = y ~ x, data = sam)
Residuals:
Min 1Q Median 3Q Max
-0.92770 -0.20715 0.03368 0.20038 0.90512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.614428 0.344622 4.685 9.04e-06 ***
x 0.014476 0.003212 4.508 1.82e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3346 on 98 degrees of freedom
Multiple R-squared: 0.1717, Adjusted R-squared: 0.1633
F-statistic: 20.32 on 1 and 98 DF, p-value: 1.818e-05
> data.frame(head(res), head(summary(mo.xy)$residuals))
head.res. head.summary.mo.xy..residuals.
1 0.119501750 0.119501750
2 0.215614401 0.215614401
3 0.009628661 0.009628661
4 -0.473936002 -0.473936002
5 -0.004950905 -0.004950905
6 0.102940633 0.102940633
> data.frame(min(res), max(res))
min.res. max.res.
1 -0.9276983 0.9051242
>
> # residual의 분산값은 무엇인가?
> # residual의 집합인 res 변인의 ss(res) 값을
> # ss.res 라고 하면 이를 n-2로 나눠준 값을 말한다.
> var.res <- ss.res / (n-2) # variance of residuals
> sd.res <- sqrt(var.res) # standard deviation of residuals
> # 위의 값을 standard deviation of residual이라고 부를 수 있다
> data.frame(var.res, sd.res)
var.res sd.res
1 0.1119626 0.3346081
> sigma(mo.xy)
[1] 0.3346081
> summary(mo.xy)$sigma
[1] 0.3346081
> summary(mo.xy)
Call:
lm(formula = y ~ x, data = sam)
Residuals:
Min 1Q Median 3Q Max
-0.92770 -0.20715 0.03368 0.20038 0.90512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.614428 0.344622 4.685 9.04e-06 ***
x 0.014476 0.003212 4.508 1.82e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3346 on 98 degrees of freedom
Multiple R-squared: 0.1717, Adjusted R-squared: 0.1633
F-statistic: 20.32 on 1 and 98 DF, p-value: 1.818e-05
>
> r.cal
[1] 0.414393
> r.sq
[1] 0.1717215
> r.cal^2
[1] 0.1717215
> ss.reg / ss.tot
[1] 0.1717215
>
> 1-r.sq
[1] 0.8282785
> ss.res/ss.tot
[1] 0.8282785
> 1-(ss.reg/ss.tot)
[1] 0.8282785
>
> # 바로 위에서 한 residual 집합의 분산과 표준편차값
> ss.res # ss.res 값
[1] 10.97233
> var.res <- ss.res/(n-2)
> var.res
[1] 0.1119626
> sqrt(var.res)
[1] 0.3346081
> sd.res
[1] 0.3346081
> ms.res <- var.res
>
> # 위의 var.res = ms.res 이라고도 부른다
> # 그리고 이것은 ss.tot - ss.reg 방법 말고도
> summary(mo.xy)$residuals
. . . . . .
> res.1 <- summary(mo.xy)$residuals
> # 이것을 우리는 이미 위에서
> # res <- y - y.hat 으로도 구했다.
> y.hat
. . . .
> mo.xy$fitted.values
. . . .
> data.frame(y.hat, mo.xy$fitted.values)
y.hat mo.xy.fitted.values
1 3.073421 3.073421
2 3.390841 3.390841
3 3.290064 3.290064
4 3.285737 3.285737
5 3.252709 3.252709
6 3.226239 3.226239
> sum((res.1-mean(res.1))^2)
[1] 10.97233
> ss.res
[1] 10.97233
>
>
> # 그런데 위의 식에서 ss.res 대신에 아래처럼
> # ss.res 이 ss.total에서 차지하는 비율로
> # 보면
> ss.res/ss.tot
[1] 0.8282785
> 1-(ss.reg/ss.tot)
[1] 0.8282785
> 1-r.sq
[1] 0.8282785
>
>
> sd.res <- sqrt(ss.res/(n-2))
> # 위의 ss.res 대신에
> # 아래처럼 ss.res/ss.tot 을 쓴다면
> # ss.res의 proportion값을 (전체 ss 값에 대한 비율값을)
> # 사용한 것이 된다. 즉, 일종의 표준화된 퍼센티지로
> # 계산을 하는 것이 된다.
> sqrt((ss.res/ss.tot)/(n-2))
[1] 0.09193379
> # 그리고 위는 아래와 같은 것
> sqrt((1-r.sq)/(n-2))
[1] 0.09193379
> # 이것을 se of residual 이라고 부른다
> se.res <- sqrt((ss.res/ss.tot)/(n-2)) # se for residual
> se.res
[1] 0.09193379
>
> # 그리고 이 se.res은 r 값의 (correlation coefficient)
> # significance 정도를 가늠하는데 쓰인다.
> # 즉, r / se.res 이는
> # t distribution을 따르기에
> t.r <- r.cal/se.res
> p.r <- pt(t.r, n-2, lower.tail = F)*2
> cat(r.cal, se.res, t.r, p.r)
0.414393 0.09193379 4.507515 1.817763e-05> cor.test(x,y)
Pearson's product-moment correlation
data: x and y
t = 4.5075, df = 98, p-value = 1.818e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.2372888 0.5648366
sample estimates:
cor
0.414393
>
> # 그렇다면 b1의 significance한 정도는 어떻게 가늠하는가?
> # standard error for b1 (regression slope)
>
> # 위에서 우리는 ms.res <- ss.res/(n-2)
> ms.res
[1] 0.1119626
> sqrt(ms.res)
[1] 0.3346081
> sd.res <- sqrt(ms.res) # residual sd (sd for regression)
> sd.res
[1] 0.3346081
>
> ss.x <- ss(x)
> ss.x
[1] 10855.25
> se.b <- sd.res/sqrt(ss.x)
> se.b
[1] 0.003211564
> sqrt(ms.res)/sqrt(ss.x)
[1] 0.003211564
> # 위의 값이 standard error of b 값이므로 .. ..
> # http://commres.net/regression#standard_error_of_b
> b.cal
[1] 0.01447618
> t.b <- b.cal/se.b
> p.b <- pt(t.b, n-2, lower.tail = F)*2
> data.frame(b.cal, se.b, t.b, p.b)
b.cal se.b t.b p.b
1 0.01447618 0.003211564 4.507515 1.817763e-05
> data.frame(r.cal, se.res, t.r, p.r )
r.cal se.res t.r p.r
1 0.414393 0.09193379 4.507515 1.817763e-05
> # 왜 t.b와 t.r 값이 같은가?
> #
> summary(mo.xy)
Call:
lm(formula = y ~ x, data = sam)
Residuals:
Min 1Q Median 3Q Max
-0.92770 -0.20715 0.03368 0.20038 0.90512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.614428 0.344622 4.685 9.04e-06 ***
x 0.014476 0.003212 4.508 1.82e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3346 on 98 degrees of freedom
Multiple R-squared: 0.1717, Adjusted R-squared: 0.1633
F-statistic: 20.32 on 1 and 98 DF, p-value: 1.818e-05
>
> data.frame(ss.reg, ss.res, ss.tot, r.sq, ss.res/ss.tot)
ss.reg ss.res ss.tot r.sq ss.res.ss.tot
1 2.274822 10.97233 13.24715 0.1717215 0.8282785
> data.frame(f.cal, df.reg, df.res, p.f)
f.cal df.reg df.res p.f
1 20.3177 1 98 1.817763e-05
> anova(mo.xy)
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x 1 2.2748 2.27482 20.318 1.818e-05 ***
Residuals 98 10.9723 0.11196
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> # remember ss.res?
> ms.res
[1] 0.1119626
> ms.reg
[1] 2.274822
>
multiple regression
rs_multiple_regression
############################
# multiple regression
############################
rm(list=ls())
# 1. mean of three variables
# iq, allowance, gpa
# x, x2, y
means <- c(108, 42, 3.2)
# allowance, iq, gpa
sds <- c(11, 8, 0.4) # standard deviations
# correlation matrix. note that correlation
# between v1 and v2 is near none (0.05)
corr_matrix <- matrix(c(1, 0.00, 0.3,
0.00, 1, 0.2,
0.3, 0.2, 1),
nrow = 3) # Target correlation matrix
# 2. covariance matrix
# diagonal matrix of SDs
sd_diag <- diag(sds)
sd_diag
# Calculate covariance matrix
cov_matrix <- sd_diag %*% corr_matrix %*% sd_diag
# 3. population data
# set.seed(32)
# set.seed(43211122222)
set.seed(1230)
n_sam <- 1000000
# n_sam <- 300
sam <- mvrnorm(n = n_sam, mu = means, Sigma = cov_matrix)
sam <- as.data.frame(sam)
colnames(sam) <- c("x", "x2", "y")
cor(sam)
y <- sam$y
x <- sam$x
x2 <- sam$x2
# y regression with x and x2
m.y.x1 <- lm(y~x, data=sam)
m.y.x2 <- lm(y~x2, data=sam)
# summary(m.y.x1)
# summary(m.y.x2)
explained.by.x1 <- summary(m.y.x1)$r.square
explained.by.x2 <- summary(m.y.x2)$r.square
round(explained.by.x1 + explained.by.x2,3)
rs.each.sum <- explained.by.x1 + explained.by.x2
m.y.x1x2 <- lm(y~x+x2, data=sam)
summary(m.y.x1x2)
rs.both <- summary(m.y.x1x2)$r.square
round(rs.both,3)
round(data.frame(rs.both, rs.each.sum, rs.both-rs.each.sum),3)
m.x1.x2 <- lm(x~x2, data=sam)
# R-squared = 0
# x2 = 0
# p-value = .7945
summary(m.x1.x2)
summary(m.x1.x2)$r.square
sqrt(summary(m.x1.x2)$r.square)
cor(x,x2)
cor.test(x,x2)
#################
# 그림으로 보기
#################
#########################
# multiple regression 2
#########################
rm(list=ls())
set.seed(1230)
# 1. mean of three variables
# iq, allowance, gpa
# x, x2, y
means <- c(108, 42, 3.2)
# allowance, iq, gpa
sds <- c(11, 8, 0.4) # standard deviations
# correlation matrix. note that correlation
# between v1 and v2 is near none (0.05)
corr_matrix <- matrix(c(1, 0.30, 0.4,
0.30, 1, 0.2,
0.4, 0.2, 1),
nrow = 3) # Target correlation matrix
# 2. covariance matrix
# diagonal matrix of SDs
sd_diag <- diag(sds)
sd_diag
# Calculate covariance matrix
cov_matrix <- sd_diag %*% corr_matrix %*% sd_diag
# 3. population data
# set.seed(32)
# set.seed(43211122222)
set.seed(1230)
n_sam <- 100
# n_sam <- 300
sam <- mvrnorm(n = n_sam, mu = means, Sigma = cov_matrix)
sam <- as.data.frame(sam)
colnames(sam) <- c("x", "x2", "y")
cor(sam)
y <- sam$y
x <- sam$x
x2 <- sam$x2
lm.y.x1 <- lm(y~x, data=sam)
lm.y.x2 <- lm(y~x2, data=sam)
summary(lm.y.x1)
summary(lm.y.x2)
bc <- summary(lm.y.x1)$r.square
cd <- summary(lm.y.x2)$r.square
bccd <- bc+cd
data.frame(bc, cd, bccd)
lm.y.x1x2 <- lm(y~x+x2, data=sam)
summary(lm.y.x1x2)
bcd <- summary(lm.y.x1x2)$r.square
# where does this diff come from?
b <- spcor.test(y,x,x2)$estimate^2
d <- spcor.test(y,x2,x)$estimate^2
b
d
bcd
bccd
c1 <- bccd-bcd
c2 <- bcd-(b+d)
b
c1
c2
d
b+c1+d
bcd
lm.x2.x1 <- lm(x2 ~ x, data=sam)
dg.res <- summary(lm.x2.x1)$residuals
lm.y.dg <- lm(y~dg.res, data=sam)
d1 <- summary(lm.y.dg)$r.square
d1
d
b
c1
d
b+c1+d
summary(lm(y~x+x2, data=sam))$r.square
ro_multiple_regression
> ############################
> # multiple regression
> ############################
> rm(list=ls())
> # 1. mean of three variables
> # iq, allowance, gpa
> # x, x2, y
> means <- c(108, 42, 3.2)
>
> # allowance, iq, gpa
> sds <- c(11, 8, 0.4) # standard deviations
> # correlation matrix. note that correlation
> # between v1 and v2 is near none (0.05)
> corr_matrix <- matrix(c(1, 0.00, 0.3,
+ 0.00, 1, 0.2,
+ 0.3, 0.2, 1),
+ nrow = 3) # Target correlation matrix
> # 2. covariance matrix
> # diagonal matrix of SDs
> sd_diag <- diag(sds)
> sd_diag
[,1] [,2] [,3]
[1,] 11 0 0.0
[2,] 0 8 0.0
[3,] 0 0 0.4
>
> # Calculate covariance matrix
> cov_matrix <- sd_diag %*% corr_matrix %*% sd_diag
>
> # 3. population data
> # set.seed(32)
> # set.seed(43211122222)
> set.seed(1230)
> n_sam <- 1000000
> # n_sam <- 300
> sam <- mvrnorm(n = n_sam, mu = means, Sigma = cov_matrix)
> sam <- as.data.frame(sam)
> colnames(sam) <- c("x", "x2", "y")
> cor(sam)
x x2 y
x 1.0000000000 -0.0002605277 0.3001789
x2 -0.0002605277 1.0000000000 0.1993134
y 0.3001788706 0.1993133956 1.0000000
>
> y <- sam$y
> x <- sam$x
> x2 <- sam$x2
>
> # y regression with x and x2
> m.y.x1 <- lm(y~x, data=sam)
> m.y.x2 <- lm(y~x2, data=sam)
> # summary(m.y.x1)
> # summary(m.y.x2)
> explained.by.x1 <- summary(m.y.x1)$r.square
> explained.by.x2 <- summary(m.y.x2)$r.square
> round(explained.by.x1 + explained.by.x2,3)
[1] 0.13
> rs.each.sum <- explained.by.x1 + explained.by.x2
>
> m.y.x1x2 <- lm(y~x+x2, data=sam)
> summary(m.y.x1x2)
Call:
lm(formula = y ~ x + x2, data = sam)
Residuals:
Min 1Q Median 3Q Max
-1.77065 -0.25170 0.00031 0.25145 1.83133
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.604e+00 4.168e-03 384.8 <2e-16 ***
x 1.091e-02 3.389e-05 321.9 <2e-16 ***
x2 9.956e-03 4.658e-05 213.8 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3732 on 999997 degrees of freedom
Multiple R-squared: 0.1299, Adjusted R-squared: 0.1299
F-statistic: 7.462e+04 on 2 and 999997 DF, p-value: < 2.2e-16
> rs.both <- summary(m.y.x1x2)$r.square
> round(rs.both,3)
[1] 0.13
>
> round(data.frame(rs.both, rs.each.sum, rs.both-rs.each.sum),3)
rs.both rs.each.sum rs.both...rs.each.sum
1 0.13 0.13 0
>
> m.x1.x2 <- lm(x~x2, data=sam)
> # R-squared = 0
> # x2 = 0
> # p-value = .7945
> summary(m.x1.x2)
Call:
lm(formula = x ~ x2, data = sam)
Residuals:
Min 1Q Median 3Q Max
-51.050 -7.432 -0.005 7.432 54.696
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.080e+02 5.877e-02 1838.108 <2e-16 ***
x2 -3.581e-04 1.374e-03 -0.261 0.794
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 11.01 on 999998 degrees of freedom
Multiple R-squared: 6.787e-08, Adjusted R-squared: -9.321e-07
F-statistic: 0.06787 on 1 and 999998 DF, p-value: 0.7945
> summary(m.x1.x2)$r.square
[1] 6.787467e-08
> sqrt(summary(m.x1.x2)$r.square)
[1] 0.0002605277
> cor(x,x2)
[1] -0.0002605277
> cor.test(x,x2)
Pearson's product-moment correlation
data: x and x2
t = -0.26053, df = 999998, p-value = 0.7945
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.002220491 0.001699438
sample estimates:
cor
-0.0002605277
> #################
> # 그림으로 보기
> #################
>
>
>
>
> #########################
> # multiple regression 2
> #########################
> rm(list=ls())
> set.seed(1230)
> # 1. mean of three variables
> # iq, allowance, gpa
> # x, x2, y
> means <- c(108, 42, 3.2)
>
> # allowance, iq, gpa
> sds <- c(11, 8, 0.4) # standard deviations
> # correlation matrix. note that correlation
> # between v1 and v2 is near none (0.05)
> corr_matrix <- matrix(c(1, 0.30, 0.4,
+ 0.30, 1, 0.2,
+ 0.4, 0.2, 1),
+ nrow = 3) # Target correlation matrix
> # 2. covariance matrix
> # diagonal matrix of SDs
> sd_diag <- diag(sds)
> sd_diag
[,1] [,2] [,3]
[1,] 11 0 0.0
[2,] 0 8 0.0
[3,] 0 0 0.4
>
> # Calculate covariance matrix
> cov_matrix <- sd_diag %*% corr_matrix %*% sd_diag
>
> # 3. population data
> # set.seed(32)
> # set.seed(43211122222)
> set.seed(1230)
> n_sam <- 100
> # n_sam <- 300
> sam <- mvrnorm(n = n_sam, mu = means, Sigma = cov_matrix)
> sam <- as.data.frame(sam)
> colnames(sam) <- c("x", "x2", "y")
> cor(sam)
x x2 y
x 1.0000000 0.4219157 0.3535740
x2 0.4219157 1.0000000 0.2700868
y 0.3535740 0.2700868 1.0000000
>
> y <- sam$y
> x <- sam$x
> x2 <- sam$x2
>
> lm.y.x1 <- lm(y~x, data=sam)
> lm.y.x2 <- lm(y~x2, data=sam)
> summary(lm.y.x1)
Call:
lm(formula = y ~ x, data = sam)
Residuals:
Min 1Q Median 3Q Max
-1.03671 -0.26187 0.03189 0.21676 0.71688
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.943438 0.347683 5.590 2.05e-07 ***
x 0.011705 0.003128 3.742 0.000308 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3753 on 98 degrees of freedom
Multiple R-squared: 0.125, Adjusted R-squared: 0.1161
F-statistic: 14 on 1 and 98 DF, p-value: 0.0003079
> summary(lm.y.x2)
Call:
lm(formula = y ~ x2, data = sam)
Residuals:
Min 1Q Median 3Q Max
-0.89560 -0.22360 -0.02491 0.25421 0.86636
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.695389 0.198771 13.560 < 2e-16 ***
x2 0.013381 0.004819 2.777 0.00658 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3863 on 98 degrees of freedom
Multiple R-squared: 0.07295, Adjusted R-squared: 0.06349
F-statistic: 7.711 on 1 and 98 DF, p-value: 0.006575
> bc <- summary(lm.y.x1)$r.square
> cd <- summary(lm.y.x2)$r.square
> bccd <- bc+cd
> data.frame(bc, cd, bccd)
bc cd bccd
1 0.1250146 0.0729469 0.1979615
>
> lm.y.x1x2 <- lm(y~x+x2, data=sam)
> summary(lm.y.x1x2)
Call:
lm(formula = y ~ x + x2, data = sam)
Residuals:
Min 1Q Median 3Q Max
-0.91766 -0.23439 0.01484 0.21300 0.72503
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.875581 0.349193 5.371 5.36e-07 ***
x 0.009650 0.003432 2.811 0.00597 **
x2 0.007287 0.005137 1.419 0.15921
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3734 on 97 degrees of freedom
Multiple R-squared: 0.1428, Adjusted R-squared: 0.1251
F-statistic: 8.08 on 2 and 97 DF, p-value: 0.0005682
> bcd <- summary(lm.y.x1x2)$r.square
> # where does this diff come from?
>
> b <- spcor.test(y,x,x2)$estimate^2
> d <- spcor.test(y,x2,x)$estimate^2
>
> b
[1] 0.06985246
> d
[1] 0.01778476
> bcd
[1] 0.1427994
> bccd
[1] 0.1979615
> c1 <- bccd-bcd
> c2 <- bcd-(b+d)
> b
[1] 0.06985246
> c1
[1] 0.05516214
> c2
[1] 0.05516214
> d
[1] 0.01778476
> b+c1+d
[1] 0.1427994
> bcd
[1] 0.1427994
>
> lm.x2.x1 <- lm(x2 ~ x, data=sam)
> dg.res <- summary(lm.x2.x1)$residuals
> lm.y.dg <- lm(y~dg.res, data=sam)
> d1 <- summary(lm.y.dg)$r.square
> d1
[1] 0.01778476
> d
[1] 0.01778476
>
> b
[1] 0.06985246
> c1
[1] 0.05516214
> d
[1] 0.01778476
> b+c1+d
[1] 0.1427994
> summary(lm(y~x+x2, data=sam))$r.square
[1] 0.1427994
>
c/ms/2026/schedule/w10.lecture.note.txt · Last modified: by hkimscil