Differences

This shows you the differences between two versions of the page.

--- b:head_first_statistics:correlation_and_regression [2023/12/12 23:39] – [exercise] hkimscil
+++ b:head_first_statistics:correlation_and_regression [2023/12/13 04:33] (current) – [exercise] hkimscil
@@ Line 162: / Line 162: @@
 ===== exercise =====
 <code>
+# correlation 에 대한 이해를 돕기 위한 연습
+# stat first 책에서 radiation exposure 양과
+# (독립변인) weight 간의 관계 (종속변인)
+#
+# 두 변인 (독립변인과 종속변인) 모두
+# 숫자로 측정된 변인
 re <- c(4, 4.5, 5, 5.5, 6, 6.5, 7)
 we <- c(12, 10, 8, 9.5, 8, 9, 6)
+# 손으로 variance, standard deviation 구하기
+# 우선 각각의 평균
 mean.re <- mean(re)
 mean.we <- mean(we)
-sd.re <- sd(re)
+mean.re
-sd.we <- sd(we)
+mean.we
-var.re <- var(re)
-var.we <- var(we)
+# sum of square x (ssx) 와 ssy 구하기
 ss.re <- sum((re-mean.re)^2)
 ss.we <- sum((we-mean.we)^2)
-sp <- sum((re-mean.re)*(we-mean.we))
+ss.re
-df <- 7-1
+ss.we
-cov.rewe <- sp/df
-cov.rewe2 <- cov(re,we)
-mean.re
+# df 값은
-mean.we
+df <- length(we)-1
-sd.re
+df
-sd.we
+# variance와 sd
+var.re.1 <- ss.re/df
+var.we.1 <- ss.we/df
+sd.re.1 <- sqrt(var.re.1)
+sd.we.1 <- sqrt(var.we.1)
+# R의 펑션을 이용해서 구하기
+var.re <- var(re)
+var.we <- var(we)
+sd.re <- sd(re)
+sd.we <- sd(we)
+var.re.1
+var.we.1
 var.re
 var.we
-ss.re
+sd.re.1
-ss.we
+sd.we.1
+sd.re
+sd.we
+# sum of product
+sp <- sum((re-mean.re)*(we-mean.we))
+cov.rewe.1 <- sp/df
+cov.rewe <- cov(re,we)
 sp
-df
+cov.rewe.1
 cov.rewe
-cov.rewe2
+# 교재에 기술된
 # r, b, r square
+# r = cov(x,y)/(sd(x)*sd(y))
+#   = sp / sqrt(ss(x)*ss(y))
 r.val <- sp / (sqrt(ss.re*ss.we))
+r.val2 <- cor(re,we)
 r.val
+r.val2
+# 기울기 b = sp / ss(x)
 b <- sp / ss.re
 b
 a <- mean.we - (b * mean.re)
 a
+# r 제곱값 (r^2)
 r.sq <- r.val^2
 r.sq
-cor(re,we)
+# 위의 모든 정보가 lm function의 결과에
 m1 <- lm(we~re)
 summary(m1)
+## summary(m1)에서 se가 나타나는 부분 설명
+# se of regression = sqrt of mse
+# p = # of parameters (IVs)
+#
+p <- 1
+n <- length(we)
+n
+df.res <- n-p-1
+# or
+m1$df.residual
+# residual 값들을 제곱해서 모두 더한 후
+# n으로 (df.residual값으로) 나누준 것
+# = 분산을 구하는 형식이다. 즉 아래는
+# regression line으로 예측하고 남은
+# 나머지 오차들에 (residuals) 대한 분산값
+# 이것을 mean square error라고 부른다
+# variance를 mean square라고 부르는 사람들이
+# 있다.
+# ss of square residuals = ssres
+ssres <- sum(m1$residuals^2)
+# ssres = sum of square error 값으로 부르기도
+# 한다.
+sse <- ssres
+mse <- sse/df.res
+mse
+# 그리고 그 값에 sqrt값을 씌워 sd값을
+# stanard deviation 구한다. 이것을
+# root mean square error라고 부른다.
+rmse <- sqrt(mse)
+rmse
+# summary(m1)에서 Residual standard error: 1.177
+# 에 해당하는 값
+summary(m1)
+# b coefficient 에 대한 standard error 값은
+# sqrt( (1/(n-2)) * (sse/ssx))
+se.b <- sqrt( (1/(n-2)) * (sse/ss.re) )
+se.b
+# 아래 아웃풋에서 Std. Error 참조
+summary(m1)
+c <- qt(.975, 5)
+se.b2 <- c*se.b
+b - se.b2
+b + se.b2
+confint(m1)
 </code>
 <code>
+> # correlation 에 대한 이해를 돕기 위한 연습
+> # stat first 책에서 radiation exposure 양과
+> # (독립변인) weight 간의 관계 (종속변인)
+> #
+> # 두 변인 (독립변인과 종속변인) 모두
+> # 숫자로 측정된 변인
 > re <- c(4, 4.5, 5, 5.5, 6, 6.5, 7)
 > we <- c(12, 10, 8, 9.5, 8, 9, 6)
 >
+> # 손으로 variance, standard deviation 구하기
+> # 우선 각각의 평균
 > mean.re <- mean(re)
 > mean.we <- mean(we)
-> sd.re <- sd(re)
-> sd.we <- sd(we)
-> var.re <- var(re)
-> var.we <- var(we)
-> ss.re <- sum((re-mean.re)^2)
-> ss.we <- sum((we-mean.we)^2)
-> sp <- sum((re-mean.re)*(we-mean.we))
-> df <- 7-1
-> cov.rewe <- sp/df
-> cov.rewe2 <- cov(re,we)
->
 > mean.re
 [1] 5.5
 > mean.we
 [1] 8.928571
-> sd.re
+>
-[1] 1.080123
+> # sum of square x (ssx) 와 ssy 구하기
-> sd.we
+> ss.re <- sum((re-mean.re)^2)
-[1] 1.88035
+> ss.we <- sum((we-mean.we)^2)
-> var.re
-[1] 1.166667
-> var.we
-[1] 3.535714
 > ss.re
 [1] 7
 > ss.we
 [1] 21.21429
-> sp
+>
-[1] -10
+> # df 값은
+> df <- length(we)-1
 > df
 [1] 6
+> # variance와 sd
+> var.re.1 <- ss.re/df
+> var.we.1 <- ss.we/df
+> sd.re.1 <- sqrt(var.re.1)
+> sd.we.1 <- sqrt(var.we.1)
+> # R의 펑션을 이용해서 구하기
+> var.re <- var(re)
+> var.we <- var(we)
+> sd.re <- sd(re)
+> sd.we <- sd(we)
+>
+> var.re.1
+[1] 1.166667
+> var.we.1
+[1] 3.535714
+> var.re
+[1] 1.166667
+> var.we
+[1] 3.535714
+> sd.re.1
+[1] 1.080123
+> sd.we.1
+[1] 1.88035
+> sd.re
+[1] 1.080123
+> sd.we
+[1] 1.88035
+>
+>
+> # sum of product
+> sp <- sum((re-mean.re)*(we-mean.we))
+> cov.rewe.1 <- sp/df
+> cov.rewe <- cov(re,we)
+> sp
+[1] -10
+> cov.rewe.1
+[1] -1.666667
 > cov.rewe
 [1] -1.666667
-> cov.rewe2
-[1] -1.666667
 >
+>
+> # 교재에 기술된
 > # r, b, r square
+> # r = cov(x,y)/(sd(x)*sd(y))
+> #   = sp / sqrt(ss(x)*ss(y))
 > r.val <- sp / (sqrt(ss.re*ss.we))
+> r.val2 <- cor(re,we)
 > r.val
 [1] -0.8206099
+> r.val2
+[1] -0.8206099
+>
+> # 기울기 b = sp / ss(x)
 > b <- sp / ss.re
 > b
@@ Line 259: / Line 380: @@
 > a
 [1] 16.78571
+> # r 제곱값 (r^2)
 > r.sq <- r.val^2
 > r.sq
 [1] 0.6734007
 >
-> cor(re,we)
+> # 위의 모든 정보가 lm function의 결과에
-[1] -0.8206099
 > m1 <- lm(we~re)
 > summary(m1)
@@ Line 280: / Line 401: @@
 re           -1.4286     0.4449  -3.211  0.02371 *
 ---
-Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+Signif. codes:
+‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 Residual standard error: 1.177 on 5 degrees of freedom
@@ Line 287: / Line 409: @@
 >
-> # regression formula
+> ## summary(m1)에서 se가 나타나는 부분 설명
-> # we hat = 16.7857 - 1.4286 re
+> # se of regression = sqrt of mse
+> # p = # of parameters (IVs)
 > #
-> m1$coefficients
+> p <- 1
-(Intercept)          re
+> n <- length(we)
-.785714   -1.428571
-> a <- m1$coefficients["(Intercept)"]
-> b <- m1$coefficients["re"]
-> predicted <- a + b*re
-> predicted
-[1] 11.071429 10.357143  9.642857  8.928571  8.214286  7.500000
-[7]  6.785714
-> m1$fitted.values
-         2         3         4         5         6
-.071429 10.357143  9.642857  8.928571  8.214286  7.500000
-.785714
-> t1 <- data.frame(predicted, m1$fitted.values)
-> t1
-  predicted m1.fitted.values
-11.071429        11.071429
-10.357143        10.357143
-  9.642857         9.642857
-  8.928571         8.928571
-  8.214286         8.214286
-  7.500000         7.500000
-  6.785714         6.785714
->
-> resid <- we-predicted
-> resid
-[1]  0.9285714 -0.3571429 -1.6428571  0.5714286 -0.2142857
-[6]  1.5000000 -0.7857143
-> m1$residuals
-          2          3          4          5
-.9285714 -0.3571429 -1.6428571  0.5714286 -0.2142857
-          7
-.5000000 -0.7857143
-> t2 <- data.frame(resid, m1$residuals)
-> t2
-       resid m1.residuals
-  0.9285714    0.9285714
--0.3571429   -0.3571429
--1.6428571   -1.6428571
-  0.5714286    0.5714286
--0.2142857   -0.2142857
-  1.5000000    1.5000000
--0.7857143   -0.7857143
->
-> sse <- sum(resid^2)
-> n <- length(re)
 > n
 [1] 7
-> sse
+> df.res <- n-p-1
-[1] 6.928571
+> # or
-> ss.re
+> m1$df.residual
-[1] 7
+[1] 5
-> se.of.b1 <- sqrt( (1/(n-2))*(sse/ss.re) )
+>
-> se.of.b1
+> # residual 값들을 제곱해서 모두 더한 후
+> # n으로 (df.residual값으로) 나누준 것
+> # = 분산을 구하는 형식이다. 즉 아래는
+> # regression line으로 예측하고 남은
+> # 나머지 오차들에 (residuals) 대한 분산값
+> # 이것을 mean square error라고 부른다
+> # variance를 mean square라고 부르는 사람들이
+> # 있다.
+> # ss of square residuals = ssres
+> ssres <- sum(m1$residuals^2)
+> # ssres = sum of square error 값으로 부르기도
+> # 한다.
+> sse <- ssres
+> mse <- sse/df.res
+> mse
+[1] 1.385714
+>
+> # 그리고 그 값에 sqrt값을 씌워 sd값을
+> # stanard deviation 구한다. 이것을
+> # root mean square error라고 부른다.
+> rmse <- sqrt(mse)
+> rmse
+[1] 1.177164
+>
+> # summary(m1)에서 Residual standard error: 1.177
+> # 에 해당하는 값
+> summary(m1)
+Call:
+lm(formula = we ~ re)
+Residuals:
+       2       3       4       5       6       7
+.9286 -0.3571 -1.6429  0.5714 -0.2143  1.5000 -0.7857
+Coefficients:
+            Estimate Std. Error t value Pr(>|t|)
+(Intercept)  16.7857     2.4872   6.749  0.00108 **
+re           -1.4286     0.4449  -3.211  0.02371 *
+---
+Signif. codes:
+‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+Residual standard error: 1.177 on 5 degrees of freedom
+Multiple R-squared:  0.6734,	Adjusted R-squared:  0.6081
+F-statistic: 10.31 on 1 and 5 DF,  p-value: 0.02371
+>
+> # b coefficient 에 대한 standard error 값은
+> # sqrt( (1/(n-2)) * (sse/ssx))
+> se.b <- sqrt( (1/(n-2)) * (sse/ss.re) )
+> se.b
 [1] 0.444926
+> # 아래 아웃풋에서 Std. Error 참조
 > summary(m1)
@@ Line 360: / Line 490: @@
 re           -1.4286     0.4449  -3.211  0.02371 *
 ---
-Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+Signif. codes:
+‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 Residual standard error: 1.177 on 5 degrees of freedom
@@ Line 366: / Line 497: @@
 F-statistic: 10.31 on 1 and 5 DF,  p-value: 0.02371
+>
+> c <- qt(.975, 5)
+> se.b2 <- c*se.b
+> b - se.b2
+[1] -2.57229
+> b + se.b2
+[1] -0.2848526
+> confint(m1)
+.5 %     97.5 %
+(Intercept) 10.39213 23.1792968
+re          -2.57229 -0.2848526
 >
 </code>
+===== exercise 2: 직접해보기 =====
+아래에서처럼 데이터를 R에서 다운로드 받아서 정리한 후에 위에서처럼 sp, ss.x, ss.y, b, a, r, r squared, 등등을 구해 보시오.
 <code>
 ########################
-ss <- c(1.9, 2.5, 3.2, 3.8, 4.7, 5.5, 5.9, 7.2)
+dat  <- read.csv("http://commres.net/wiki/_media/regression01-bankaccount.csv")
-at <- c(22, 33, 30, 42, 38, 49, 42, 55)
+# data
-mean.ss <- mean(ss)
+# bankaccount = 통장갯수
-mean.at <- mean(at)
+# income = 수입
-ss.ss <- sum((ss-mean.sa)^2)
+# famnum = 부양가족수
-ss.at <- sum((at-mean.at)^2)
+# IV = 수입 = income
-df <- 8-1
+# DV = 통장갯수 = bankaccount
-var.ss <- ss.ss/df
+#
-var.at <- ss.at/df
+# 컬럼 이름 바꾸기 (간단하게)
-sd.ss <- sqrt(var.ss)
+colnames(dat) <- c("y", "x", "x2")
-sd.at <- sqrt(var.at)
+dat
-sp.ssat <- sum((ss-mean.ss)*(at-mean.at))
+attach(dat)
-cov.ssat <- sp.ssat/df
+dat
+# ss.y
+# ss.x
+# df.y
+# df.x
+# sp.xy
+# sd.x
+# sd.y
+# b coefficient
+# b
+# a intercept
+# a
+# pred <- a + b*x
+# resid <- y - pred
+# ss.pred
+# ss.resid
+# ss.y
+# r
+# r.sq
 </code>