首先先明確幾個概念(例子以SLR說明,MLR同理)

1. Y=β₀+β₁X₁+β₂X₂+...β_pX_p+e

此式代表著數(shù)據(jù)點的值,也就是觀測到的值, p代表predictor的個數(shù), X_p在矩陣中實際上為X_i,p, i代表第i個數(shù)據(jù)
這里β₀,β₁,e等都是未知的,只有觀測到的每個數(shù)據(jù)(Y_i,X₁,X₂,...X_p)是已知的

2. β₀,β₁與

β₀代表Y軸上的截距,可能并沒有含義(具體看定義)
β₁代表X₁每上升1,其它X固定,對應Y上升量為β₁
由于β₀,β₁是未知的,我們采用來模擬β₀,β₁的真實值獲得擬合直線.顯然當帶入1.的公式中,此時的Y也應該用來表示. 這里實際表示,也可寫成,表示X=x時的平均值
例如: 假設我們有10個病人,我們想知道吃完降壓藥后血壓下降和時間是否有線性關(guān)系,X₁可以取幾個值: 10,20,30,50min.... 通過公式得到的實際上代表10個人血壓在40min時的平均值.

3.

顯然僅能作為預測的樣本平均值.對于每一個個體如果我們選取X₀,對應Y用來表示.雖然我們不知道的具體值,但是我們能估計其分布:
(1)均值:均值為**,且為的無偏估計;
(2)95%CI: X=x₀時對應的95%CI 為
(3)方差:variance由兩部分組成: 個體的差異(如每個人血壓差異,同一個人每次量血壓也有差異) + 抽樣樣本對真實均值的偏離

具體怎么算可以自己推導,這里很容易理解:σ^2不會變化,但抽樣的樣本數(shù)量增加,樣本均值會更接近總體均值,會變小

使用隨機生成的5個數(shù)據(jù)作為示例

x1 <- rnorm(5)
x2 <- rnorm(5) #按標準正態(tài)分布生成5個數(shù)字
y <- c(1,2,3,4,5)
xy.lm <- lm(y~x1+x2) #選擇因變量和自變量

#直接進行分析
summary(xy.lm)
#或
anova(xy.lm)

結(jié)果解讀:

image.png

跳過summary先講Sum Sq和Mean Sq
SSEresidual(殘差)的平方和,即
SSR來自回歸模型(regression)
SSY/SST 來自Y的總體:
MSR=SSR/p
MSE=SSE/(n-1-p)
F=MSR/MSE 如果β₁=0,模型會是一條水平直線,導致MSR=0.也就是說越是β₁≠0,F會越大,就可以拒絕H₀,但不能說明擬合程度好
R²: R²=(SSY-SSE)/SSY=SSR/SSY 意思是X回歸的方差占Y的總方差的比例,代表在回歸直線上的集中程度 (簡單理解為誤差比例越小,擬合度越好)
Adjusted R²: R²_adj=1-(SSE/(n-p-1))/(SST/(n-1)) 可按預測自變量數(shù)目p進行調(diào)整. 當X很多時R²會過大導致無法代表擬合程度,但是會限制額外的自變量.
注意anova()使用的是Type I,和變量的順序有關(guān).更常使用anova(lm, type="III") 或anova(model1, model2) 研究每一個變量對模型的貢獻

分步計算過程

找找對應關(guān)系
想想怎么計算的均值,95%CI和標準差 (公式見開頭)

將XY轉(zhuǎn)化為矩陣,構(gòu)造design matrix

Y<-matrix(y)
X <- matrix(c(rep(1,5), x1, x2), nrow=5,byrow=F)
Xt <- t(X)

求H matrix;I,J輔助陣

矩陣相乘使用%*%
H <- X%*%solve(Xt%*%X)%*%Xt
J_n<-matrix(rep(1,5*5),nrow=5)
I_n <-diag(5)

Sum square

SStot <- t(Y)%*%(I_n - J_n/5)%*%Y #J/5 對應5個數(shù)據(jù)
SSreg <- t(Y)%*%(H - J_n/5)%*%Y
SSerr <- t(Y)%*%(I_n - H)%*%Y
print("SST    SSR    SSE");print(paste(SStot,SSreg,SSerr,sep="    "))

print("MSE");print(MSE<- as.numeric(SSerr/(5-3)))
varhate<- MSE*(I_n-H)

image.png

求, COV(β), SE of β

betahat <- solve((Xt%*%X))%*%Xt%*%Y 
print("betahat");print(betahat)
print("Cov of beta");print(cov_betas <- MSE*solve(t(X)%*%X))
print("se of beta");print(se_betas <- sqrt((diag(diag(cov_betas)))))

image.png

求

yhat<-X%*%betahat
    #or yhat<-H%*%Y
print("Yhat");print(yhat)

F statistic (總體)

print(MSreg <- SSreg/2) #2個X
Fstat <- MSreg/MSE
print(Fstat)
pvalF <- 1-pf(Fstat, 1, 5-3)
print(paste(Fstat,pvalF))

image.png