> library(openxlsx)
> data = read.xlsx("22_data.xlsx",sheet = 1)
> x = data[,-c(1,2)]
> x = cbind(rep(1,17),x)
> x_mat  = as.matrix(x)
> y =matrix(data[,2],ncol = 1)
> res = solve(t(x_mat)%*%x_mat)%*%t(x_mat)%*%y
> res[,1]
rep(1, 17) 19412.8597818
X1             0.2679605
X2            -0.2874013
X3          -297.3653736

> lm(y~x_mat)Call:
lm(formula = y ~ x_mat)Coefficients:(Intercept)  x_matrep(1, 17)          x_matX1
19412.859781545               NA      0.267960511  x_matX2          x_matX3  -0.287401287   -297.365373557

library(graphics)
pairs(data[,-1]pch = 21,bg = c('red','green3','blue'))
# pch参数是控制点的形状，bg是控制点的颜色

> summary(lm(y~x_mat+1))Call:
lm(formula = y ~ x_mat + 1)  #调用Residuals:      #残差统计量，残差第一四分位数(1Q)和第三分位数(3Q)有大约相同的幅度，意味着有较对称的钟形分布Min      1Q  Median      3Q     Max
-4397.9 -1102.4   153.8  1184.4  2934.6 Coefficients: (1 not defined because of singularities)        Estimate Std. Error t value Pr(>|t|)
(Intercept)      1.941e+04  3.524e+04   0.551    0.591
x_matrep(1, 17)         NA         NA      NA       NA
x_matX1          2.680e-01  4.466e-02   6.000 4.45e-05 ***
x_matX2         -2.874e-01  1.668e-01  -1.723    0.109
x_matX3         -2.974e+02  3.688e+02  -0.806    0.435
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#标记为Estimate的列包含由最小二乘法计算出来的估计回归系数。#标记为Std.Error的列是估计的回归系数的标准误差。#从理论上说，如果一个变量的系数是0，那么该变量将毫无贡献。然而，这里显示的系数只是估计，它们不会正好为0.#因此，我们不禁会问：从统计的角度而言，真正的系数为0的可能性有多大？这是t统计量和P值的目的，在汇总中被标记为t value和Pr(>|t|)
#P值估计系数不显著的可能性，有较大P值的变量是可以从模型中移除的候选变量

Residual standard error: 2013 on 13 degrees of freedom
Multiple R-squared:  0.9982,    Adjusted R-squared:  0.9977
F-statistic:  2348 on 3 and 13 DF,  p-value: < 2.2e-16

#Residual standard error 表示残差的标准差，F-statistic 表示F的统计量

x_mat = cbind(rep(1,17),data[,3])
y = data[,2]
res = lm(y~x_mat)
res
> resCall:
lm(formula = y ~ x_mat)Coefficients:
(Intercept)       x_mat1       x_mat2  -6213.0189           NA       0.1915

t = data[,1]
y = data[,2]
res = lm(y~t)
res> resCall:
lm(formula = y ~ t)Coefficients:
(Intercept)            t  -16428607         8213  

> newdata = data.frame(t = 2017)
> pre = predict(res,newdata,interval = "prediction",level = 0.95)
> prefit      lwr      upr
1 136337.8 116018.1 156657.4