title: "Analysis of Variance"
author: '01'
date: "2022-11-23"
output: html_document
editor_options:
  markdown:
    wrap: 72


library(dplyr)
library(ggplot2)
library(knitr)
library(palmerpenguins)
library(tidyr)
library(ISLR2)

Analysis of Variance 方差分析

Assumptions

1. 每个样本服从正态分布。

2. 每个样本方差\(\sigma^2\)相同。

3. 样本之间相互独立 levels：a factor has many levels.

   要做的事情：比较各总体（每个水平下的总体）的均值是否相等

举例：比较三种疗法的好坏 如果是两种，可以使用t test

这里，可以把三种疗法看作3个因子，（因子变量/分类变量），响应变量：治疗效果

Treatment is a between-groups factor with 3 levels.

This model is called one-way ANOVA, bacause there's only one

classification variable.

So this is a one-way between-groups ANOVA.

这里可以看作是独立的，对于其中两组可以做独立T检验。

对于某一种疗法，下面引入一个新的变量------时间。例如：在5周之后做测量，在6个月之后做测量。

Time is within-groups factor with 2 levels.

This model is called one-way within-groups ANOVA. （因为疗法都是一样的）

这里不是独立的，对于其中两组做T检验时，paired=TRUE

如果同时考虑疗法和时间： This model is a two-factor ANOVA.

对于两个因子，可以考虑它们时间是否有交互作用。

Eg. 5种药物

library(tidyverse)
data(cholesterol, package = "multcomp")
plotdata <- cholesterol %>%
  group_by(trt) %>%
  summarize(n = n(), #n()表示出现系数。
  mean = mean(response),
  sd = sd(response))
plotdata


ggplot(plotdata, aes(x = trt, y = mean, group = 1)) +
  geom_point(size = 3, color = "red") +
  geom_line(linetype = "dashed", color = "darkgrey") +
  geom_errorbar(aes(ymin = mean - qt(0.975, df = n - 1) * sd / sqrt(n),
ymax = mean + qt(0.975, df = n - 1) * sd / sqrt(n)), width = .1) +
 labs(x = "Treatment", y = "Response", title = "Mean plot with 95% confidence interval")

# geom_errorbar()是在算95%的置信区间

One-way ANOVA

fit <- aov(response ~ trt, data = cholesterol)
summary(fit)

原假设：所有方法都一样 备择假设：至少有一个方法不一样

所以此处可以拒绝原假设

\[Y_{ij} = \mu + \alpha_j + \epsilon+{ij}$$, j stands for class j, i
stands for sample i in class j.

$\mu$ = 总体的均值

$\alpha_j$ = $\mu_j - \mu$

$K$ = total number of class

$n_j$ = total number of samples in class j

$n$ = total number of samples

$\epsilon_{ij}$ 这里假设是正态的，均值为0，方差相等为$\sigma^2$

$$Y_{ij} = \mu + 0\alpha_1+ … + 1\alpha_j+…+0\alpha_K + \epsilon+{ij}\]

于是，我们对这样一个模型做linear regression：

\[Y = \mu + X_1\alpha_1+ … + X_j\alpha_j+…+X_K\alpha_K + \epsilon
\]

当Y属于第j类时，\(X_j\) = 1，否则\(X_j\) = 0

注意在这里\(\mu\)可以直接用总体来估计。 并且\(\sum{\alpha_{ij}} = 0\)

\(\hat b = ( X^T X) X^T Y\)

mean(cholesterol$response)
fit2 <- lm(response - mean(response) ~ 0 + trt, data = cholesterol) #这里的结果表示

fit2 <- lm(response  ~ 0 + trt, data = cholesterol) #这里的结果表示：trt1time的效果比平均效果差6.95

fit2 <- lm(response  ~ trt, data = cholesterol) #这里的结果表示：截距就是trt1time的疗效。trt1time的结果是5.78，trt2time的结果比trt2time的高3.443
fit2

要验证\(H_0: \alpha_1 = \alpha_2 = \alpha_3 = \alpha_4 = \alpha_5 = 0\)

\(H_1\) = 至少有一个\(\alpha_j\)不为0。

这里可以选择验证\(\sum {\alpha_i^2} = 0\)

\[\frac{1}{\sigma^2}\sum_{j = 1}^{K} n_j (\bar Y.,j - \bar Y ..)^2
\]

它服从卡方分布。 d.f. = K - 1

这里方差未知，需要估计$$\hat \sigma^2 = s^2 = \frac{1}{n-K} \sum \sum (Y_{ij} - \bar Y_{ij})^2 $$它服从卡方分布。

d.f. = n - K

用这两个卡方分布去构造一个F分布的统计量。

说白了就是一开始就想知道$$\frac{1}{\sigma^2}\sum_{j = 1}^{K} n_j (\bar Y.,j - \bar Y ..)^2$$是不是接近0，所以去找了它的分布，然后发现它是一个卡方分布，然后就可以构造求p值

$$$$

我们发现在$$Y = \mu + \alpha_j + \epsilon_{ij}$$

中前两项是已知的，最后一项可以生成，所以我们可以生成Y，

\(sdf\)

Eg. two-way factorial ANOVA

data(ToothGrowth)
ToothGrowth$dose <- factor(ToothGrowth$dose) # what happens if we remove this line?
stats <- ToothGrowth %>%
  group_by(supp, dose) %>%
  summarize(n = n(),
  mean = mean(len),
  sd = sd(len),
  std.error = sd/sqrt(n))

stats


fit1 <- aov(len ~ supp * dose, data = ToothGrowth)
summary(fit1)
fit2 <- aov(len ~ supp + dose, data = ToothGrowth)
summary(fit2)

注意：把dose看作double和把dose看作一个因子，得到的是有差异的。

单因子H_0：所有的alpha都是0 多因子H_0：所有的alpha都是0，

\(y\) supp: \alpha_1 -> OJ, \alpha_2 -> VC dose: \beta_1 -> 0.5,

\beta_2 -> 1 \beta_3 -> 1

\(y_{ijk}\), i-th supp, j-th dose, k-th sample

\(y_{ijk} = \mu + \alpha_i + \beta_j + \epsilon_{ijk}\)

\(\mu + \alpha_1\) 表示OJ维生素的样本的均值

同样假设 \(\epsilon_{ijk} \sim N(0, \sigma^2)\)

\(n_{ij} = 10, i \in (0,2], j \in (0,3]\) the number of samples that take

i-th vitamin and j-th dose

\[\hat \mu = \frac{1}{n} \sum_{ijk} y_{ijk}
\]

\[\hat \alpha_i = \frac{1}{n_{11} + n_{12} + n_{13}} (\sum_{jk} y_{ijk}) - \hat \mu
\]

\(H_0: \alpha_1 = \alpha_2 = 0\) \(\alpha_1^2 + \alpha_2^2\) 小即可

即\(\sum_{i=1}^2 (y_i - \bar y)^2\) 就是上面表里的205，sum of square, d.f.

= 2-1 = 1 (\(\alpha_1 + \alpha_2 = 0\)) Residual sum of square:

\(\sum (y_{ijk} - \bar y_{ij.})^2}\)

total sum of square \(\sum_{ijk}(y_{ijk} - \bar y_{…})^2\) d.f. = 60 - 1

= 59，最后residual的那个d.f. = 59 - 1 - 2 = 56.

# validate 交互效应
#y ~ supp + dose + supp:dose
#也就是乘法

同样的剂量和不同的维生素之间是有差异的，

如果没有交互效应，两条线应该是平行的。

data(ToothGrowth)
lm(len ~ supp + dose, data = ToothGrowth)

ToothGrowth$dose = as.double(ToothGrowth$dose)
lm(len ~ supp + dose, data = ToothGrowth)

intercept表示当不用suppVC也不用dose就是9.272

suppVC表示，使用suppVC效果会差3.7

dose表示，固定其他的东西的情况下，多使用一个单位的剂量，效果会好9，764

这里因为factor类型的只有一个supp, 所有条件 if OJ y = 9.272 + 9.764_dose

if VC y = 9.272 - 3.700 + 9.764_dose

注意如果当dose是分类变量的时候会不一样。会出现dose1，dose2，相应的所有变量的解释会变得不一样。

同时，d.f.也会不一样。

summary： perdictor-response 1. quan-quan :regression 2. cate-quan :

anova 3. mix-quan: anova

when response is categorical.

\(E(Y|X) = \beta_0 + \beta_1 X\) \(Var(Y|X) = \sigma^2\)

P(Y = 1|X) = \(\beta_0 + \beta_1 X\) P(Y = 0|X)

用线性模型来拟合概率，但是线性模型求出的概率可能为负数。

这时，我们使用logit。 利用logistic function\(y = \frac{e^x}{1+e^x}\)

此时\(P(Y = 1|X) = \frac{e^{\beta_0+\beta_1X}}{1 + e^{\beta_0 + \beta_1 X}}\)

Eg. Default dataset

data(Default)
b <- ggplot(Default[1:1000,],aes(x=balance,y=income))
b + geom_point(aes(color=student))

Eg. IRIS dataset

b <- ggplot(iris,aes(x=Sepal.Length,y=Sepal.Width))
b + geom_point(aes(color=Species))

response variable Y = 1,2,3 它们之间实际上没有顺序关系，不能用线性回归，

Logistic regression

\[p(X)=\frac{e^{\beta_0+\beta_1X}}{1+e^{\beta_0+\beta_1X}}
\]

odds = $\frac{p(X)}{1 - p(X)} $

表示划分为一类的概率和划分为另一类的概率之比。

\[\frac{p(X)}{1 - p(X)} = e^{\beta_0+\beta_1X}
\]

X每增加一个单位，odds就增加\(e^{\beta_1}\)倍。

这里既然涉及概率，用Maximum Likelihood来求解。

联合概率是：\(P(Y_1 = y_1,… Y_n = y_n | X_1,…,X_n) = \Pi_{i=1}^n P(Y_i = y_i | X_i = x_i)\)

其中 if \(y_i = 0\)

\(P(Y_i = y_i | X_i = x_i) = P(Y_i = 0 | X_i = x_i) = 1 - P(Y_i = 1 | X_i = x_i) = 1 - \frac{e^{\beta_0+\beta_1 x_i}}{1+e^{\beta_0+\beta_1 x_i}} = \frac{1}{1+e^{\beta_0+\beta_1 x_i}}\)

if \(y_i = 1\)

\(P(Y_i = y_i | X_i = x_i) = P(Y_i = 0 | X_i = x_i) = 1 - P(Y_i = 1 | X_i = x_i) = 1 - \frac{e^{\beta_0+\beta_1 x_i}}{1+e^{\beta_0+\beta_1 x_i}} = \frac{1}{1+e^{\beta_0+\beta_1 x_i}}\)

\[\mathcal{l}(\beta_0,\beta_1)=\Pi_{i:y_i = 1}p(x_i)\Pi_{i \prime:y_{i \prime = 0}} (1 - p(x_{i \prime}))
\]

没有解析解，只能求数值解，

set.seed(1)
n <- nrow(Default)
relabel <- sample(n,n)
train.index <- relabel[1:9000]
test.index <- relabel[9001:10000]
train.data <- Default[train.index,]
glm.fit <- glm(student~income,data=Default,family="binomial")
glm.fit$coefficients

如果family = "gaussian"就是linear regression

Hypothesis tests

\(H_0: \beta_1 = 0\)

z-score = \(\frac{\hat \beta_1}{SE(\hat \beta_1)}\)

还是想知道\(\beta_1\)的分布，但是只能知道它的渐进分布为正态分布。

summary(glm.fit)



x.new=data.frame(c(10000,40000))
colnames(x.new) = 'income' #要把列名强调好
predict(glm.fit,newdata = x.new,type = "response")


test.data <- Default[test.index,]
glm.prob <- predict(glm.fit,test.data,type="response")
glm.pred <- rep("No",1000)
glm.pred[glm.prob > 0.5] <- "Yes"
table(glm.pred,test.data$student)

mean(glm.pred != test.data$student)

Multiple Logistic Regression

\[log ( \frac{p(X)}{1 - p(X)} ) = \beta_0 + \beta_1X_1+…+\beta_pX_p
\]

\(\beta_1\)表示固定其他的不变，在X_1增加一个单位时，左边会增大\(e^\{\beta_1}\)倍。

\[p(X) = \frac{e^{\beta_0 + \beta_1X_1+…+\beta_pX_p}}{1+e^{\beta_0 + \beta_1X_1+…+\beta_pX_p}}
\]

glm.fit2 <- glm(student~balance+income,data=train.data,family='binomial')
glm.fit2$coefficients

summary(glm.fit2)$coef


test.data <- Default[test.index,]
glm.prob2 <- predict(glm.fit2,test.data,type="response")
glm.pred2 <- rep("No",1000)
glm.pred2[glm.prob2 > 0.5] <- "Yes"
table(glm.pred2,test.data$student)

mean(glm.pred2 != test.data$student)

发现增加了一个predictor效果还变差了。需要找个一个评估模型的方法。例如：cross-validation

Multinomial Logistic Regression

library(nnet)
multinom.fit <- multinom(Species ~ ., data=iris,trace = FALSE)
summary(multinom.fit)$coefficients

table(iris$Species,predict(multinom.fit))