Handout 1: Introduction and null model

• Working with R and R markdown
• The null model for iid normal data, \(y \sim N(\mu,\sigma^2)\).
• The null linear model gives a one-sample t-test of \(H_0: \mu=0\).
• Manually calculating the t-statistic
  • \(t-stat = (\hat{mean} - \mu_0)/se\)
• Manually calculating the 95% confidence interval for \(\mu\).
  • \(tmult = qt(1-.05/2, df=n-1)\)
  • \(se(y) = sd(y)/\sqrt{(n)}\)
  • \(CI = mean(y) + c(-1,+1) \times tmult \times se(y)\)
• The paired t-test is a one-sample test of within-pair differences
  • 配对t检验，其实就是将成对样本只差作为单样本的t检验。
• Our Executive Summary makes conclusions about the expected value, \(\mu\).
• expected mean，总体均值\(\mu\)
• 中央极限定理CLT：总体容量足够大的时候，我们能取到的所有样本的均值近似正态分布

• T型分布

  • \(\frac{\hat Y - \mu}{\frac{S}{\sqrt n}} \sim T_{n-1}\)
  • 自由度为n-1的t型分布 在概率论和统计学中，学生t-分布（Student’s t-distribution），可简称为t分布，用于根据小样本来估计呈正态分布且方差未知的总体的均值。如果总体方差已知（例如在样本数量足够多时），则应该用正态分布来估计总体均值。 它是对两个样本均值差异进行显著性测试的学生t检定的基础。学生t检定改进了Z检定（Z-test），因为Z检定以母体标准差已知为前提。虽然在样本数量大（超过30个）时，可以应用Z检定来求得近似值，但Z检定用在小样本会产生很大的误差，因此必须改用学生t检定以求准确。
• Null model summary 只能给出总体均值的期望值还有它为0的概率

• t.test 可以指定\(\mu\) 所以可以用来检测样本均值跟假定均值是否有差距。

Handout 2: Simple linear model

• A single numeric explanatory variable: \[ y_i = \beta_0 + \beta_1 x_i + \epsilon_i \]
• \(R^2\)
• Model checking:
  • Independence of observations
  • Constant scatter and no trend in residuals
  • Residuals approximately normal
  • Checking for (and possibly removing) any highly influential points
• Confidence and prediction intervals using the predict function
  • Confidence intervals are for estimating an expected value
  • Prediction intervals are for predicting a new \(y\) value.
• Writing the equation for the model
• Our Executive Summary makes conclusions about the expected change in \(y\) for a change in \(x\).
  • Typically not of interest to report about the intercept (\(\beta_0\)).

Handout 3: Fitting curves using linear models

• If the scatter plot of \(x\) and \(y\) suggests nonlinearity
  • the residual plot from a simple linear model will display stronger nonlinearity
  • 拟合新型模型后的残差图 do not appear to randomly scattered around 0.
  • 中间低两边高总之不是直线
• Add a quadratic term \[ y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \epsilon_i \]
• Determine whether the quadratic term is needed, by testing \(H_0: \beta_2=0\).
  • If not signficant, the simple linear model is preferred.
  • 如果它看起来像个曲线那就fit一个曲线，只要响应变量恒定的散布在曲线周围即可。

Handout 4: Categorical explanatory variable:

• A single categorical explanatory variable with two levels (i.e., two-sample t-test situation)
• Expressing the model using dummy variables (that takes the value 0 or 1): \[ y_i = \beta_0 + \beta_1 D_i + \epsilon_i \]
• Interpreting the coefficient \(\beta_1\) relative to the baseline \(\beta_0\).

Handout 5: Multiplicative models — working with logged response variable 乘法模型

• \(y=\alpha \cdot {(e^{\beta_1})}^x\)
• Symptoms that indicate we should workwith \(\log y\) include:
  • \(y\) values are right-skewed(右偏，又称为正偏) around trend line and variability increases with increasing \(y\).
  • Scatter plot of \(x\) vs \(y\) showing exponential increase or decrease
  • Our intuition telling us that effects will operate on the multiplicative scale
  • So a log-linear model was fitted
• Fit an appropriate linear model to \(\log y\) \[log(y_i)=\beta_0+\beta_1x_i+\epsilon_i\]
• Backtransform CIs using \(exp\).
• Interpretation is on the multiplicative scale — e.g., percentage change
• The percentage change is in the median of y.

Handout 6: Categorical and numeric explanatory variables (a.k.a. ANCOVA model)

• One numeric explanatory variable and one 2-level categorical/factor explanatory variable (i.e., grouping variable).
• Fit interaction model
• Dummy variables for both intercept and slope: \[ y_i = \beta_0 + \beta_1 x_i + \beta_2 D_i + \beta_3 D_i x_i + \epsilon_i \]
• Interpreting the coefficients \(\beta_2\) and \(\beta_3\) relative to the coefficients of the baseline level (i.e, group 1) of the categorical variable:
  • \(\beta_2\) is change in intercept relative to baseline level
  • \(\beta_3\) is change in slope (i.e., change in effect of \(x\)) relative to baseline.
• Check p-value of \(\beta_3\) to determine if interaction is significant – remove if not.
• Model with interaction corresponds to different slopes in the two groups
  • Typically not of interest to report on intercept coefficients.

Handout 7: Categorical and numeric explanatory variables (a.k.a. ANCOVA model) 分类和数值变量

• Example where the categorical variable has 3 levels
• Use anova function to jointly test the statistical significance of the interaction terms.
  • 用anov 函数来一并测试全部交互项的统计学显著程度，是否有common slope。
  • 它可以同时检查两个或者更多个系数是否为零。
  • 如果交互项是显著的（i.e., 系数不为零）则即使该项的斜率为0，也不能去掉该项；
  • 与二次方模型相同，如果有二次方项，则不能去掉其1次方项，即使它的技术为0.
• Apply Occam’s razor (i.e., simplify the model by removing interaction term) if interaction not significant.
  • This gives the main-effects model (a.k.a. parallel lines).
    • 主效应模型，即不管属于哪个类别，随数字变量的变化效果都相同。
  • Further simplify model if common slope not significantly different from zero, or difference in intercepts not significant.
    • 如果这个共同的斜率common slope等于0，不管属于哪个分类，岁数字变量的变化都为0，也就是说分类变量已经没有作用了，见CS7.2 — 亦即，分类变量对响应变量没有影响，无须解释。
• For main-effects model:（有common slope \(\neq 0\)）
  • Report on the slope (which is common to all groups)
  • Report on the difference in intercepts (i.e, the vertical separation between the parallel lines)

Handout 8: Models with several explanatory variables 多重回归模型

• Use pairs20x plot to look at patterns and strength of correlations
  • 对于pairs20x
  • 对角线是概率密度分布图
  • 左下角是各变量的系数，表明他们的关系强度
  • 右上角是y相对x的散点图
• 本质上，解释变量的统计意义（P-Value）衡量的是：在响应变量的变化性被其他解释变量调整以后，该解释变量还能解释这个变异性的大小
• Using confounding (highly correlated) explanatory variables causes problems:
  • Loose statistical signficance even when explanatory variables have clear association with \(y\).
  • Can be hard to interpret the effects
• Avoid using highly correlated explanatory variables.
• Interpret the effect as conditional on all other explanatory terms in the model

Handout 9: Power law models 幂律模型

• Response increases/decreases to some power of the numeric explanatory variable: \(y=\alpha x^{\beta_1}\)
• 相应变量增加或减少到解释变量的n次幂
  • Typically, the variability increases with increasing \(y\).
• Take the log of both \(y\) and \(x\), and fit a linear model: \[ \log(y_i) = \beta_0 + \beta_1 \log(x_i) + \epsilon_i \]
• Interpret \(\beta_1\) as a percentage change in the median of \(y\) from a 1% increase in \(x\).
  • Or, could report percentage change due to doubling of \(x\), say.

Handout 10: 多级分类变量以及多级比较问题

• \(lm(y ~ fact1)\)
• One-Way ANOVA 只有一个解释变量，不存在interaction
• 能够通过ANOVA 输出计算\(R^2\)
  • \(R^2 = 100 \times (1 - \frac{"Residuals \, Sum \, Sq"}{"Sum\,Sq列之和""})\)
• 多个levels要同时进行比较，否则:
  • We have 5% chance of erroneously making an untrue conclusion each time.
  • This is known as multiple comparison problem
  • There will be erroneous evidence of effect of …
  • So, we need to use founction multipleComp(data.fit)
• 步骤：
  • iid检查
  • anova(fit) 检查P-Value 得出各组数据是否存在差异，可计算\(R^2\)
  • summary1way(fit) 可以得出总体均值，Table of Effects：各组均值相对总体均值的偏差（reference是总体均值，不是G1）
  • multipleComp(fit) 得出所有比较结果， 提取其中P-value小于0.05 的项用来解释
  • 如果y值有取过对数，则记得做exp逆运算。
  • Executive Summary:

Handout 11: 两个分类变量

• ANOVA 的解释：ANOVA（方差分析）主要描述解释变量全部是分类的线性回归。
• 两个解释变量用 Two-Way ANOVA
• 步骤：
  • interactionPlots(y ~ x1 + x2, data = df)
    • 看不到平行线 parallel lines indicating interaction between two factors
  • fit with interaction lm(y~x1*x2, data = df)
  • iid check
  • anova(fit) 检查相关性是否确立
  • summary2way(fit, page="interaction") 如果有interaction
  • Executive Summary: 我们常常这样解释
    • We estimate that for someone who …,(因为互相影响，所以要先设定一个相同的前提下讨论另一个区别) those who with … can expect to …. more than those who ….
  • anova 还能确认其中一个factor 是否有影响（就是看这个factor里面不同组对应的响应变量是否有所不同。），如果P-Value小于0.05，则去掉一个解释变量，变成One-Way Anova，如果只有两个levels，又变成了two-sample t-test.

Handout 12: 分析计数值数据

dpois(12, 9.61)

## [1] 0.08685078

rpois(20, 10)

##  [1] 10 14 13  6 13  4 12  8  5  9 10  6 10  8 14  8  6  7 10 12

barplot(dpois(1:12, 3), names=1:12, las=1)

barplot(dpois(0:250,100), names=0:250, las=1)

AP.tb1 = matrix(c(17,83,27,19), nrow = 2, byrow = T, dimnames = list(c("attand", "not.attend"),c("fail","pass")))
# AP.tb1
# rowdistr(AP.tb1)
# show to show row proportions

# chisq.test(AP.tb1)
# chisq.test(AP.tb1, correct = F)

AP.chisq=chisq.test(AP.tb1, correct = F)
AP.chisq$expected

##                fail     pass
## attand     30.13699 69.86301
## not.attend 13.86301 32.13699

X = (17-30.14)^2/30.14 + (83-69.86)^2/69.86 + (27-13.86)^2/13.86 + (19-32.14)^2/32.14
X

## [1] 26.02961

chisq.test(AP.tb1, correct = F)

## 
##  Pearson's Chi-squared test
## 
## data:  AP.tb1
## X-squared = 26.016, df = 1, p-value = 3.386e-07

rowdistr(AP.tb1, comp = "between", plot = F)

## Row Proportions
##            fail pass Totals   n
## attand     0.17 0.83      1 100
## not.attend 0.59 0.41      1  46
## 
## 95% CIs for diffs between proportions with fac2 = fail 
## (rowname-colname) 
##        not.attend      
## attand (-0.577,-0.257)*
## 
## 95% CIs for diffs between proportions with fac2 = pass 
## (rowname-colname) 
##        not.attend    
## attand (0.257,0.577)*

rowdistr(AP.tb1, comp = "between", suppressText = T)

Handout 13： 比值比

AP.tb1

##            fail pass
## attand       17   83
## not.attend   27   19

1. 手动求CI for OR的公式 知道contingency table 可以交叉求出\(\hat{OR}\),
   注意排列顺序不同\(\hat{OR}\)可能不同，但也可能相同。
   
2. \(\hat{OR} = \frac{\hat{Odds_1}}{\hat{Odds_2}} = \frac{n_{11}n_{22}}{n_{12}n_{21}}\)
   • 这个\(\hat{OR}\) 算出来就可以知道同在列1组内不同行事件概率的比值了，而且可以推广到不同table排列方式的比值
   • \(se\left(\log(\hat{OR})\right) = \sqrt{\frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}}}\) 根据STATS703
     
   • 当样本比较大的时候\(\log(\hat{OR})\) 接近正态分布 is approximately normal,
   • so CI for \(\log(\hat{OR})\):
     • \(\Rightarrow \log(\hat{OR}) \pm 1.96 \times se\left(\log(\hat{OR})\right) = (a,b)\)
       
     • \(\Rightarrow\) 95% CI for \(\hat{OR}\):
       
     • \(\exp(a), \exp(b)\) 用以解释数据的行列比较
3. Test the association using the log odds: \(H_0\) means no row-column association \(\Rightarrow \hat{OR}=1 \Rightarrow log(\hat{OR})=0\) To test \(H_0\) we use Z statisic \(Z = \frac{\log(\hat{OR}) - 0}{\rm se(\log(\hat{OR}))}\) \(\Rightarrow P-value 进而判断 \hat{OR}=1 是否极端\) 其实这个P-Value就是测试\(H_0:log(\hat{OR}=0)\) 跟第15章的 summary(glm.fit) 出来的第四二个 \(P-value 测试 coefficient[4] \beta_3=0\) 是一样的，p值也一样 z-test 适用于大型数据，多余30个个体以上，如果知道标准差，那最好用z!

Handout 14： 费舍尔检验 FET

• 如果行列不关联，即选中观察值的完全是凭运气猜的。
• 能够选择到这样的表格的概率\(P_1\)+所有比这个表格更难猜到的表格概率\(P_2, P_3, \cdots\) 之和 要小于0.05 那么我们就说这个表格不是乱猜出来的。
• 这说明行列有关联，是因为这个关联性才得出这样的统计数字。

Handout 15： 推广线性统计模型到关联列表

• 多重线性模型提供了很多便利性：
  • simple linear regression
  • ANOVA two-way ANOVA summary2way()
  • ANCOVA
• 如何让多重线性模型在响应变量是count的时候仍然有效？
• log不是好办法
• 泊松分布
• 相关性：
  • 关联列表中没有相关性，第一行不同列的期望值差，与第二行内不同列期望值之差相同，就说他们有相关性，这是不对的。
  • 要看他们不同期望值之间的倍数关系而不是之差。
  • 这个倍数差别相同才能说他们有相关性。
  • 如果他们在倍数尺度上的差别相同
  • 这点关联列表跟two-way ANOVA 有所不同，关联列表里面的数是两分类变量组合出现的个数。要倍数差别相同才能说没有相关性，比如来听课的和不来听课的学生，开始通过的和不通过的人数都是8：2或者说是4倍，这样才能说没有相关性。但two-way ANOVA因为其响应变量不是出现个数而是领一个值，所以只要增加减少的幅度想等就可以说没有相关性了，所以看two-way ANOVA可以直接看斜率，而关联列表不行。
• 分析步骤：

1. fisher.test(tbl)$p.val \(H_0\) 是无相关 或者 `chisq.test(tbl)$expected 倍数差距相等 则无相关
2. 模型： \(n_{11} = n \times r_1 \times c_1\)
   \(log(n_{11}) = log\,n + log\, r_1 + log\, c_1 = \beta_0\)
   \(log(n_{12}) = \beta_0 + \beta_1 \Rightarrow \beta_1 = log c_2 - log c_1\)
   \(log(n_{21}) = \beta_0 + \beta_2 \Rightarrow \beta_2 = log r_2 - log r_1\)
   \(log(n_{22}) = \beta_0 + \beta_1 + \beta_2\) no interaction
   \(log(n_{22}) = \beta_0 + \beta_1 + \beta_2 + \beta_3\) interaction
   \(\beta_3 = log(\hat{OR})\)
   
3. fit = glm(y~fac1 * fac2, family = poisson, data = df)
   
4. summary(fit) residual deviance means residual sum-of-squares 4. summary 出来就coefficient \(\beta_0, \beta_1, \beta_2\) 将他们加起来再exp 就可以得到4个格子的期望值了。似乎没什么意思，估出来跟观测值相同，因为模型中我们加入了相关性参数。
   
5. 如果\(\beta_3 \neq 0\)
   • 其中\(\beta_3\) 的解释比较多 取exp 可以得到 \(\hat{OR}\)
   • 用起来解释不同组合的odds之比
6. 如果\(\beta_3 = 0\)
   • 没有interaction的情况下，我们可以说，独立的针对行内和列内进行比较
   • \(log(n_{12}) = \beta_0 + \beta_1 \Rightarrow \beta_1 = log c_2 - log c_1 \Rightarrow exp^{\beta_1} = \frac{c_2}{c_1}\)
     
   • \(log(n_{21}) = \beta_0 + \beta_2 \Rightarrow \beta_2 = log r_2 - log r_1 \Rightarrow exp^{\beta_2} = \frac{r_2}{r_1}\)
   • c2列出现的数量是c1列出现数量的多少倍
   • r2行出现的数量是r1行出现数量的多少倍

Handout 16： 推广线性统计模型到计数数据

• 这次数据类型不是列联表了，不过他的因变量仍然是一个count
• 数据类型，只不过可以说是数据个数少了吧。

• 如何将多重线性模型在响应变量是count的时候仍然有效？

• 之前已经学习了各种情况下拟合线性模型，现在y是一个count的时候只需要在R code里稍稍改一点点，重要的是如何解释输出。
• \(log(E[Y]) = linear model\)
• 重要的是我们不需要变换Y 因为Rcode 已经内建了变换
• 以前用lm 的倍数模型拟合count 数据的问题是，Y的值如果为0，如何log(0)呢，这样就会人为的造成outlier
• fit = glm(y~f2 * x2, family = poisson, data = df)
• 1-pchisq(deviance(gfit), df = df.residual(gfit))求P值用来假设泊松假设：方差等于均值
  • 如果模型正确，残差偏差应该是卡方分布的。
  • 如果p小于0.05
  • 说明我们的模型不适当not adequate
  • family = quasipoisson
• family = quasipoisson 说明数据的方差比他们满足泊松分布应有的方差要大。
  \(Var(Y) > E(Y) \quad or\quad Var(Y) = k\,E(Y)\)
  K > 1
  
• 其结果是我们估计效果的标准差变大了，但估计效果本身没有变
• z-value 的绝对值变小了，p值变大了。

• 用lm模型来拟合log(Y) 没有考虑到数据的方差随着均值的增大而增大。同时，被迫抛弃一些有效的数据点，因为他们的cook’s distance 太大。

Handout 17.5 逻辑回归

Handout 18： 当相应变量是百分比的时候

Handout 18.5： 垃圾邮件

• 要建模我们需要一些训练数据
• 逻辑回归可以称为监管学习方法
• fit = glm(spam~.) . 超级牛，涵盖所有解释变量除了spam本身。 sapply(email, function(x){if(is.factor(x)){table(email$spam,x)}}) 找出level 中没有case 或者control的category
• anova 去掉not significant variables
  • dropColumns 去掉不显著的columns
  • refit
  • 设定分类依据 classification rules
  • 其实就是设定判断垃圾邮件的threshold
  • 这个threshold affect the “true positive rate and false negative rate”
  • `predSpam = ifelse(fitted(spam.fit2) > 0.5, 1, 0)
  • 预测的正确率：
    • 通过一个混合矩阵：confusion matrix，就是一个two-way table，行是真实值，列是估值，对照一看就知道正确率了，这其实还是回到看关联列表的相关性和相关程度上来了。有意思
    • conf.matric = table(email$spam, predSpam)
    • 这个table的，对角线之和除以所有条码之和，就得到正确率了
    • pred.accuracy = sum(diag(conf.matrix))/sum(conf.matrix)
    • tpr = conf.matrix[2,2]/sum(conf.matrix[2,])
    • 这样得出的准确率是高，但是垃圾邮件的检测成功率却不高，这是因为模型把大部分精力用在了拟合费垃圾邮件上，必经费垃圾邮件数量庞大得多。
    • 给每个判断做一个权重，spam的权重为1，不是垃圾的权重为367/3554.
    • 重新refit
    • 这样准确率降低了一些，但垃圾邮件的检测成功率大幅提升了。
  • 交叉验证
  • 用同一组数据进行建模以及对模型的验证，多半会导致对建模的效果过分乐观。最好是有一组数据专门做建模，又叫训练模型。另一组数据用于测试，又叫测试模型。但通常我们没有机会或者这些信息。怎么办？
  • 交叉验证： 用一部分portion数据用来建模，一部分用来测试，这个比例proportion 有很多种：
    • 留一法，5倍法，10倍法，k倍法，
    • 10倍法 10-fold 就是90%用来训练模型，10%用来测试模型
    • 按这个比例随机抽取数据得到两个子集subset of the data，然后进行前面的计算。
    • 讲这个随机抽取和计算重复100遍，就可以得到这个模型的性能的均值，已经这个性能的变化范围（变异性）

Handout 时间序列

data("airpass.df")
passengers.ts = ts(airpass.df$passengers, start = 1949, frequency = 12)
plot(passengers.ts, main = "", ylab = "")

passengers.stl = stl(passengers.ts, s.window = "periodic")
plot(passengers.stl)
abline(v=1950:1960, lty=2)

plot(passengers.ts, main = "", ylab = "")

remainder = passengers.stl$time.series[,3]

data("mening.df")
mening.ts = ts(mening.df$mening, frequency=12, start=c(1990,1))
plot(mening.ts, main="")
abline(v=1991:2001, lty = 2)

# airpass.df = within(airpass.df, {month = factor(rep(month.abb, 12), levels = month.abb); time = 1:144})
# lines(passengers~month, data = airpass.df)

data(beer.df)
beer.ts = ts(beer.df, frequency = 12, start = c(1971,7))
plot(beer.ts, main = "US Beer Production, July 1970 to June 1978")
abline(v = 1971:1978, lty = 2)

passengers.ts = ts(airpass.df$passengers, start = 1949, frequency = 12)
plot(passengers.ts, main = "", ylab="")
abline(v=1950:1960, lty=2)

log.passengers.ts = log(passengers.ts)
plot(log.passengers.ts, main = "", ylab="")
abline(v=1950:1960, lty=2)

passengers.stl = stl(log.passengers.ts, s.window = "periodic")
plot(passengers.stl)
title("Decomposition of log(Passengers)")

# 代码段2 ----
passengers.ts = ts(log(airpass.df$passengers), start = 1949, frequency = 12)
plot(passengers.ts, main = "Data with Seasonal Component")

passengers.stl = stl(passengers.ts, s.window = "periodic")
sa.passengers = passengers.ts - passengers.stl$time.series[, "seasonal"]
plot(sa.passengers, main = "Data with seasonal component removed")

data("larain.df")
LArain.ts = ts(rev(larain.df$LA.Rain) * 25.4, frequency = 1, start = 1877)
plot(LArain.ts, xlab = "Season", ylab = "Total rainfall (mm)")

LArain.hw = HoltWinters(LArain.ts, alpha = 0.5, beta = FALSE, gamma = FALSE)
plot(LArain.hw, lwd = 2, col.predicted = "blue", main = "alpha = 0.5")

LArain.hw = HoltWinters(LArain.ts, alpha = 0.95, beta = FALSE, gamma = FALSE)
plot(LArain.hw, lwd = 2, col.predicted = "blue", main = "alpha = 0.95")

LArain.hw = HoltWinters(LArain.ts, alpha = 0.75, beta = FALSE, gamma = FALSE)
plot(LArain.hw, lwd = 2, col.predicted = "blue", main = "alpha = 0.75")

LArain.hw = HoltWinters(LArain.ts, alpha = 0.25, beta = FALSE, gamma = FALSE)
plot(LArain.hw, lwd = 2, col.predicted = "blue", main = "alpha = 0.25")

LArain.hw = HoltWinters(LArain.ts, alpha = 0.05, beta = FALSE, gamma = FALSE)
plot(LArain.hw, lwd = 2, col.predicted = "blue", main = "alpha = 0.05")

LArain.hw1 = HoltWinters(LArain.ts, beta = FALSE, gamma = FALSE)
LArain.pred = predict(LArain.hw1, n.ahead = 5)
plot(LArain.hw1, predicted.values = LArain.pred, lwd = 2, col.predicted = "blue")

LArain.pred = predict(LArain.hw1, n.ahead = 15, prediction.interval = TRUE)
plot(LArain.hw1, predicted.values = LArain.pred, lwd = 2, col.predicted = "blue")

hist(resid(LArain.hw1))

log.LArain.ts = log(LArain.ts)
LArain.hw2 = HoltWinters(log.LArain.ts, beta = FALSE, gamma = FALSE)
plot(LArain.hw2)

log.LArain.pred = predict(LArain.hw2, n.ahead = 15, prediction.interval = TRUE)
exp(log.LArain.pred)

## Time Series:
## Start = 1943 
## End = 1957 
## Frequency = 1 
##           fit      upr      lwr
## 1943 343.7816 929.1857 127.1928
## 1944 343.7816 931.5911 126.8644
## 1945 343.7816 933.9966 126.5377
## 1946 343.7816 936.4020 126.2126
## 1947 343.7816 938.8074 125.8893
## 1948 343.7816 941.2129 125.5675
## 1949 343.7816 943.6184 125.2474
## 1950 343.7816 946.0239 124.9289
## 1951 343.7816 948.4296 124.6121
## 1952 343.7816 950.8353 124.2968
## 1953 343.7816 953.2411 123.9831
## 1954 343.7816 955.6470 123.6709
## 1955 343.7816 958.0531 123.3604
## 1956 343.7816 960.4593 123.0513
## 1957 343.7816 962.8657 122.7438

log.pass.ts = ts(log(airpass.df$passengers), frequency = 12, start = 1949)
log.pass.hw = HoltWinters(log.pass.ts)
log.pass.pred = predict(log.pass.hw, n.ahead = 30)
plot(log.pass.hw, log.pass.pred, col.predicted = "blue", lwd = 2)

log.pass.pred = predict(log.pass.hw, n.ahead = 30, prediction.interval = TRUE)
plot(log.pass.hw, log.pass.pred, col.predicted = "blue", col.intervals = "red", lwd = 2)

plot(log.passengers.ts)

# 第二步：添加季节项 ----
airpass.df = within(airpass.df, { month = factor(rep(month.abb, 12), levels = month.abb); time = 1:144})
pass.fit = lm(log(airpass.df$passengers)~month, data = airpass.df)
plot(pass.fit$residuals)
# plot(pass.fit$residuals, type = "l")
lines(airpass.df$time, residuals(pass.fit)) # 用线将点连接起来
lines(lowess(airpass.df$time, residuals(pass.fit))) # 用离散点平滑后画出的曲线，用来了解趋势再好不过了。

# lowess returns a list containing components x and y which give the coordinates of the smooth. The smooth can be added to a plot of the original points with the function lines: see the examples. 在这里显然不用了，因为散点图已经可以看出明显的趋势了。

# 第三部：添加时间项 ----
pass.fit1 = lm( log(airpass.df$passengers)~time + month, data = airpass.df )
anova(pass.fit1)

## Analysis of Variance Table
## 
## Response: log(airpass.df$passengers)
##            Df  Sum Sq Mean Sq  F value    Pr(>F)    
## time        1 25.1233 25.1233 7143.582 < 2.2e-16 ***
## month      11  2.2843  0.2077   59.047 < 2.2e-16 ***
## Residuals 131  0.4607  0.0035                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot( residuals(pass.fit1)~fitted(pass.fit1))
lines(lowess(fitted(pass.fit1), residuals(pass.fit1)))

pass.fit2 = lm( log(airpass.df$passengers)~time + month + I(time^2), data = airpass.df )
anova(pass.fit2)

## Analysis of Variance Table
## 
## Response: log(airpass.df$passengers)
##            Df  Sum Sq Mean Sq   F value    Pr(>F)    
## time        1 25.1233 25.1233 10813.900 < 2.2e-16 ***
## month      11  2.2843  0.2077    89.386 < 2.2e-16 ***
## I(time^2)   1  0.1587  0.1587    68.307 1.409e-13 ***
## Residuals 130  0.3020  0.0023                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(residuals(pass.fit2)~fitted(pass.fit2))
lines(lowess(fitted(pass.fit2), residuals(pass.fit2)))

plot(residuals(pass.fit2), type = "l") #link the points lower case L
lines(lowess(fitted(pass.fit2), residuals(pass.fit2)))

acf(residuals(pass.fit2))

acf(residuals(pass.fit2))

pass.fit3 = lm( log(passengers)[-1] ~ time[-1] + month[-1] + I(time[-1]^2) + log(passengers)[-144], data = airpass.df)

anova(pass.fit3)

## Analysis of Variance Table
## 
## Response: log(passengers)[-1]
##                        Df  Sum Sq Mean Sq  F value    Pr(>F)    
## time[-1]                1 24.4515 24.4515 19260.14 < 2.2e-16 ***
## month[-1]              11  2.2733  0.2067   162.79 < 2.2e-16 ***
## I(time[-1]^2)           1  0.1617  0.1617   127.35 < 2.2e-16 ***
## log(passengers)[-144]   1  0.1362  0.1362   107.25 < 2.2e-16 ***
## Residuals             128  0.1625  0.0013                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

acf(residuals(pass.fit3))

总结：