Ridge regression ใน r (ทีละขั้นตอน)

การถดถอยแบบริดจ์ เป็นวิธีการที่เราสามารถใช้เพื่อให้พอดีกับแบบจำลองการถดถอยเมื่อมี พหุคอลลิเนียร์ ในข้อมูล

โดยสรุป การถดถอยกำลังสองน้อยที่สุดพยายามค้นหาการประมาณค่าสัมประสิทธิ์ที่ลดผลรวมที่เหลือของกำลังสอง (RSS):

RSS = Σ(ฉัน _ฉัน – ŷ _ฉัน )2

ทอง:

• Σ : สัญลักษณ์กรีกหมายถึง ผลรวม
• y _i : ค่าตอบสนองจริงสำหรับการสังเกต ^{ครั้งที่ 3}
• ŷ _i : ค่าตอบสนองที่คาดการณ์ไว้ตามแบบจำลองการถดถอยเชิงเส้นพหุคูณ

ในทางกลับกัน การถดถอยของสันจะพยายามลดสิ่งต่อไปนี้ให้เหลือน้อยที่สุด:

RSS + ΣΣβ _เจ ²

โดยที่ j ไปจาก 1 ถึง p ตัวแปรทำนายและ แล ≥ 0

เทอมที่สองในสมการนี้เรียกว่า การลงโทษการถอน ในการถดถอยสันเขา เราเลือกค่าสำหรับ γ ที่สร้างการทดสอบ MSE ต่ำที่สุดที่เป็นไปได้ (ค่าคลาดเคลื่อนกำลังสองเฉลี่ย)

บทช่วยสอนนี้ให้ตัวอย่างทีละขั้นตอนของวิธีดำเนินการ ridge regression ใน R

ขั้นตอนที่ 1: โหลดข้อมูล

สำหรับตัวอย่างนี้ เราจะใช้ชุดข้อมูลในตัวของ R ที่เรียกว่า mtcars เราจะใช้ hp เป็นตัวแปรตอบสนอง และตัวแปรต่อไปนี้เป็นตัวทำนาย:

• mpg
• น้ำหนัก
• อึ
• คิววินาที

เพื่อทำการถดถอยสัน เราจะใช้ฟังก์ชันจากแพ็คเกจ glmnet แพคเกจนี้ต้องการให้ ตัวแปรตอบสนอง เป็นเวกเตอร์และชุดของตัวแปรทำนายต้องเป็นของคลาส data.matrix

รหัสต่อไปนี้แสดงวิธีกำหนดข้อมูลของเรา:

 #define response variable
y <- mtcars$hp

#define matrix of predictor variables
x <- data.matrix(mtcars[, c('mpg', 'wt', 'drat', 'qsec')])

ขั้นตอนที่ 2: ติดตั้งโมเดลการถดถอยของสันเขา

ต่อไป เราจะใช้ฟังก์ชัน glmnet() เพื่อให้พอดีกับโมเดลการถดถอย Ridge และระบุ alpha=0

โปรดทราบว่าการตั้งค่าอัลฟ่าเท่ากับ 1 เทียบเท่ากับการใช้การถดถอยแบบ Lasso และการตั้งค่าอัลฟ่าเป็นค่าระหว่าง 0 ถึง 1 เทียบเท่ากับการใช้ตาข่ายแบบยืดหยุ่น

โปรดทราบว่าการถดถอยแบบสันจำเป็นต้องสร้างข้อมูลให้เป็นมาตรฐาน โดยตัวแปรทำนายแต่ละตัวมีค่าเฉลี่ยเป็น 0 และค่าเบี่ยงเบนมาตรฐานเป็น 1

โชคดีที่ glmnet() จะสร้างมาตรฐานนี้ให้กับคุณโดยอัตโนมัติ หากคุณได้กำหนดมาตรฐานให้กับตัวแปรแล้ว คุณสามารถระบุ standardize=False ได้

 library (glmnet)

#fit ridge regression model
model <- glmnet(x, y, alpha = 0 )

#view summary of model
summary(model)

          Length Class Mode   
a0 100 -none- numeric
beta 400 dgCMatrix S4     
df 100 -none- numeric
dim 2 -none- numeric
lambda 100 -none- numeric
dev.ratio 100 -none- numeric
nulldev 1 -none- numeric
npasses 1 -none- numeric
jerr 1 -none- numeric
offset 1 -none- logical
call 4 -none- call   
nobs 1 -none- numeric

ขั้นตอนที่ 3: เลือกค่าที่เหมาะสมที่สุดสำหรับ Lambda

ต่อไป เราจะระบุค่าแลมบ์ดาที่ทำให้เกิดการทดสอบค่าความคลาดเคลื่อนกำลังสองเฉลี่ย (MSE) ต่ำสุดโดยใช้ การตรวจสอบความถูกต้องแบบข้าม k-fold

โชคดีที่ glmnet มีฟังก์ชัน cv.glmnet() ที่ทำการตรวจสอบข้าม k-fold โดยอัตโนมัติโดยใช้ k = 10 ครั้ง

 #perform k-fold cross-validation to find optimal lambda value
cv_model <- cv. glmnet (x, y, alpha = 0 )

#find optimal lambda value that minimizes test MSE
best_lambda <- cv_model$ lambda . min
best_lambda

[1] 10.04567

#produce plot of test MSE by lambda value
plot(cv_model) 

ค่าแลมบ์ดาที่ย่อการทดสอบ MSE ให้เหลือน้อยที่สุดกลายเป็น 10.04567

ขั้นตอนที่ 4: วิเคราะห์แบบจำลองขั้นสุดท้าย

สุดท้ายนี้ เราสามารถวิเคราะห์แบบจำลองสุดท้ายที่สร้างจากค่าแลมบ์ดาที่เหมาะสมที่สุดได้

เราสามารถใช้รหัสต่อไปนี้เพื่อรับการประมาณค่าสัมประสิทธิ์สำหรับแบบจำลองนี้:

 #find coefficients of best model
best_model <- glmnet(x, y, alpha = 0 , lambda = best_lambda)
coef(best_model)

5 x 1 sparse Matrix of class "dgCMatrix"
                    s0
(Intercept) 475.242646
mpg -3.299732
wt 19.431238
drat -1.222429
qsec -17.949721

นอกจากนี้เรายังสามารถสร้าง Trace plot เพื่อให้เห็นภาพว่าการประมาณค่าสัมประสิทธิ์เปลี่ยนแปลงไปอย่างไรเนื่องจากการเพิ่มขึ้นของแลมบ์ดา:

 #produce Ridge trace plot
plot(model, xvar = " lambda ") 

สุดท้ายนี้ เราสามารถคำนวณ R-squared ของโมเดล จากข้อมูลการฝึกได้:

 #use fitted best model to make predictions
y_predicted <- predict (model, s = best_lambda, newx = x)

#find OHS and SSE
sst <- sum ((y - mean (y))^2)
sse <- sum ((y_predicted - y)^2)

#find R-Squared
rsq <- 1 - sse/sst
rsq

[1] 0.7999513

R กำลังสองกลายเป็น 0.7999513 นั่นคือแบบจำลองที่ดีที่สุดสามารถอธิบายความแปรผันของค่าตอบสนองของข้อมูลการฝึกอบรม ได้ 79.99%

คุณสามารถค้นหาโค้ด R แบบเต็มที่ใช้ในตัวอย่างนี้ ได้ ที่นี่