Differences

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

--- b:head_first_statistics:correlation_and_regression [2023/12/13 02:30] – [Correlation] hkimscil
+++ b:head_first_statistics:correlation_and_regression [2023/12/13 12:51] – [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 67: / Line 70: @@
 \begin{align}
 r = & b * \frac {sd(x)}{sd(y)} \\
-= & b * \frac {SS_X} {SS_Y} \nonumber  \\
+= & b * \sqrt{\frac {SS_X} {SS_Y}} \;\;\;  \because \left(b = \frac {SP}{SS_X} \right) \nonumber \\
-\end{align}
+= & \frac {SP}{SS_X} * \frac {\sqrt{SS_X}} {\sqrt{SS_Y}} \nonumber \\
-\begin{align}
+= & \frac {SP}{\sqrt{SS_X * SS_Y}} \nonumber \\
-b = & b * \frac {sd(x)}{sd(y)} \\
+= & \frac {COV}{SD(X) SD(Y)}
-= & b * \frac {SS_X} {SS_Y} \nonumber  \\
+\end{align}
-\end{align}
@@ Line 158: / Line 160: @@
 >
 </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>
+<code>
+########################
+ss <- c(1.9, 2.5, 3.2, 3.8, 4.7, 5.5, 5.9, 7.2)
+at <- c(22, 33, 30, 42, 38, 49, 42, 55)
+mean.ss <- mean(ss)
+mean.at <- mean(at)
+ss.ss <- sum((ss-mean.sa)^2)
+ss.at <- sum((at-mean.at)^2)
+df <- 8-1
+var.ss <- ss.ss/df
+var.at <- ss.at/df
+sd.ss <- sqrt(var.ss)
+sd.at <- sqrt(var.at)
+sp.ssat <- sum((ss-mean.ss)*(at-mean.at))
+cov.ssat <- sp.ssat/df
+</code>