Differences

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

--- b:head_first_statistics:correlation_and_regression [2022/12/14 14:40] – [Correlation] hkimscil
+++ b:head_first_statistics:correlation_and_regression [2023/12/13 04:33] (current) – [exercise] hkimscil
@@ Line 33: / Line 33: @@
 {{:b:head_first_statistics:pasted:20221209-085250.png}}
-\begin{eqnarray*}
+\begin{align}
-b & = & \frac{\Sigma{(x-\overline(x))(y-\overline(y))}}{\Sigma{(x-\overline{x})^2}} \\
+b  = & \frac{\Sigma{(x-\overline{x})(y-\overline{y})}}{\Sigma{(x-\overline{x})^2}} \nonumber \\
-& = & \frac{\text{SP}}{\text{SS}_{x}} \\
+   = & \frac{SP}{SS_{x}} \\
-a & = & \overline{y} - b \overline{x}
+a  = & \overline{y} - b \; \overline{x} \;\;\; \because \; \overline{y} = a + b \; \overline{x}
-\end{eqnarray*}
+\end{align}
 <code>
+> # lm = linear modeling (lm function)
+> #    = regression
+> #    = correlation
 > mod <- lm(c~s, data=df)
 > summary(mod)
@@ Line 64: / Line 67: @@
 ===== Correlation =====
 Textbook says:
-$$ r = b * \frac {sd(x)}{sd(y)} $$
+\begin{align}
+r = & b * \frac {sd(x)}{sd(y)} \\
+= & b * \sqrt{\frac {SS_X} {SS_Y}} \;\;\;  \because \left(b = \frac {SP}{SS_X} \right) \nonumber \\
+= & \frac {SP}{SS_X} * \frac {\sqrt{SS_X}} {\sqrt{SS_Y}} \nonumber \\
+= & \frac {SP}{\sqrt{SS_X * SS_Y}} \nonumber \\
+= & \frac {COV}{SD(X) SD(Y)}
+\end{align}
 [[:Regression]] page says:
 \begin{eqnarray}
@@ Line 82: / Line 94: @@
 r & = & \frac {Cov(x,y)} {sd(x)*sd(y)} \\
 \end{eqnarray*}
+see [[:correlation#pearson_s_r]]
 아래는 위의 세 가지 방법으로 R에서 r 값을 구해본 것이다.
 <code>
 # check the coefficient from the lm analysis
-> summary(mod)$coefficients
+b <- summary(mod)$coefficients[2,1]
-             Estimate Std. Error  t value    Pr(>|t|)
+b
-(Intercept) 15.728319  4.4372270 3.544628 0.012150696
+df
-s            5.336411  0.9527289 5.601184 0.001379363
+sd.x <- sd(df$s)
+sd.y <- sd(df$c)
-> b <- 5.336411
+# 첫 번째
+r.value1 <- b * (sd.x/sd.y)
+r.value1
+# 두 번째
+# rsquared
+r.value2 <- sqrt(summary(mod)$r.squared)
+r.value2
+# or
+cov.sc <- cov(df$s, df$c)
+r.value3 <- cov.sc/(sd.x*sd.y)
+r.value3
+</code>
+아래는 위의 아웃풋
+<code>
+> # check the coefficient from the lm analysis
+> b <- summary(mod)$coefficients[2,1]
+> b
+[1] 5.336411
 > df
     s  c
@@ Line 102: / Line 138: @@
 5.9 42
 7.2 55
-> sd(df$s)
+> sd.x <- sd(df$s)
-[1] 1.813393
+> sd.y <- sd(df$c)
-> sd(df$c)
-[1] 10.56189
 >
-> r.value1 <- b * (sd(df$s)/sd(df$c))
+> # 첫 번째
+> r.value1 <- b * (sd.x/sd.y)
 > r.value1
-[1] 0.916219
+[1] 0.9162191
 >
+> # 두 번째
-# rsquared
+> # rsquared
 > r.value2 <- sqrt(summary(mod)$r.squared)
 > r.value2
 [1] 0.9162191
+>
-# or
+> # or
 > cov.sc <- cov(df$s, df$c)
-> sd.x <- sd(df$s)
-> sd.y <- sd(df$c)
 > r.value3 <- cov.sc/(sd.x*sd.y)
 >
 > r.value3
 [1] 0.9162191
+>
+</code>
+===== 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)
+mean.re
+mean.we
+# sum of square x (ssx) 와 ssy 구하기
+ss.re <- sum((re-mean.re)^2)
+ss.we <- sum((we-mean.we)^2)
+ss.re
+ss.we
+# df 값은
+df <- length(we)-1
+df
+# 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
+sd.re.1
+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
+cov.rewe.1
+cov.rewe
+# 교재에 기술된
+# 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
+# 위의 모든 정보가 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)
+> mean.re
+[1] 5.5
+> mean.we
+[1] 8.928571
+>
+> # sum of square x (ssx) 와 ssy 구하기
+> ss.re <- sum((re-mean.re)^2)
+> ss.we <- sum((we-mean.we)^2)
+> ss.re
+[1] 7
+> ss.we
+[1] 21.21429
+>
+> # 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
+>
+>
+> # 교재에 기술된
+> # 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
+[1] -1.428571
+> a <- mean.we - (b * mean.re)
+> a
+[1] 16.78571
+> # r 제곱값 (r^2)
+> r.sq <- r.val^2
+> r.sq
+[1] 0.6734007
+>
+> # 위의 모든 정보가 lm function의 결과에
+> m1 <- lm(we~re)
+> 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
+>
+> ## summary(m1)에서 se가 나타나는 부분 설명
+> # se of regression = sqrt of mse
+> # p = # of parameters (IVs)
+> #
+> p <- 1
+> n <- length(we)
+> n
+[1] 7
+> df.res <- n-p-1
+> # or
+> m1$df.residual
+[1] 5
+>
+> # 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)
+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
+>
+> 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>
+########################
+dat  <- read.csv("http://commres.net/wiki/_media/regression01-bankaccount.csv")
+# data
+# bankaccount = 통장갯수
+# income = 수입
+# famnum = 부양가족수
+# IV = 수입 = income
+# DV = 통장갯수 = bankaccount
+#
+# 컬럼 이름 바꾸기 (간단하게)
+colnames(dat) <- c("y", "x", "x2")
+dat
+attach(dat)
+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>