সমাধান — অধ্যায় ৫.২ · Regression Diagnostics, Inference & Model Selection¶

অধ্যায় ফাইল: part-5-modeling/05-02-regression-diagnostics-inference-selection.md (§৭ অনুশীলনী)-এর সব প্রশ্নের পূর্ণ সমাধান। সব সংখ্যাগত ফল NumPy/SciPy দিয়ে যাচাই ও statsmodels-এর সাথে মিলিয়ে দেখা হয়েছে; seeded কোড reproducible। চলমান উদাহরণ — FULL: $\text{price}\sim\text{area}+\text{age}+\text{rooms}$, true $\beta=[25,\,0.45,\,-0.8,\,6]$, seed default_rng(20260619), noise $\sigma=15$, $n=200$; predictor area$\sim U(50,200)$, age$\sim U(0,40)$, rooms$\in\{1,\dots,8\}$; অতিরিক্ত collinear predictor plot_size$\approx 1.15\cdot\text{area}$ ও row 0-এ একটি high-leverage point (area[0]=340, price[0]=250)। মূল সূত্র: hat matrix $H=X(X^\top X)^{-1}X^\top$, leverage $h_{ii}=H_{ii}$; $\operatorname{Var}(\hat\beta)=\sigma^2(X^\top X)^{-1}$, $t_j=\hat\beta_j/\widehat{\mathrm{se}}(\hat\beta_j)\sim t_{n-p}$; Cook's $D_i=\dfrac{r_i^2}{p}\cdot\dfrac{h_{ii}}{1-h_{ii}}$ (studentized $r_i$); $\text{VIF}_j=1/(1-R_j^2)$; $R^2_{\text{adj}}=1-\dfrac{\text{SSE}/(n-p)}{\text{SST}/(n-1)}$; $\text{AIC}=-2\ell+2K$, $\text{BIC}=-2\ell+(\ln n)K$।

স্মারক — fitted মান (seed 20260619, BASE মডেল price~area+age+rooms): $\hat\beta=[16.955,\,0.4956,\,-0.7203,\,6.498]$; $\widehat{\mathrm{se}}=[4.671,\,0.0258,\,0.1009,\,0.5032]$; $t=[3.63,\,19.21,\,-7.14,\,12.91]$; $R^2=0.7659$, $R^2_{\text{adj}}=0.7623$, $\hat\sigma=16.30$; $\text{SSE}=52092.4$, $\text{SST}=222476.5$, $df=196$; overall $F=213.7$ ($df=3,196$), $p\approx 1.6\times10^{-61}$।

ক · ধারণাগত (conceptual)¶

সমাধান ১ (★)¶

residual-vs-fitted plot পড়ার মূল নীতি: সুস্থ residual মানে fitted মান যাই হোক, residual-এর কেন্দ্র $0$ ও ছড়ানো একই (এলোমেলো অনুভূমিক ব্যান্ড)। যেকোনো গঠন কোনো-না-কোনো অনুমান-ভঙ্গের ইঙ্গিত।

কেস	যা দেখা যাচ্ছে	কোন অনুমান ভাঙছে	করণীয়
(ক)	এলোমেলো অনুভূমিক ব্যান্ড, প্যাটার্নহীন	কিছুই না — আদর্শ	মডেল গ্রহণযোগ্য; এগিয়ে inference
(খ)	funnel/মেগাফোন (ছড়ানো ডানে বাড়ে)	Equal variance (homoscedasticity)	response transform ($\log y$, $\sqrt y$), অথবা weighted least squares / robust (sandwich) SE
(গ)	স্পষ্ট U-আকৃতির বাঁক	Linearity	predictor-এ অরৈখিক term যোগ ($x^2$, $\log x$), অথবা মডেল-রূপ বদল
(ঘ)	বেশির ভাগ ব্যান্ডে, একটি বিন্দু বহু দূরে একা	অনুমান নয় — outlier / সম্ভাব্য influential point	Cook's $D$ ও leverage যাচাই; বিন্দুটি নিয়ে ও ছাড়া fit তুলনা, কারণ অনুসন্ধান (ডেটা-ভুল?)

মূল insight (অন্তর্দৃষ্টি): residual plot হলো অনুমানের "এক্স-রে"। funnel ⇒ E (variance), বাঁক ⇒ L (linearity); Normality আলাদা করে QQ-plot দিয়ে দেখা হয়, আর Independence সময়/ক্রমভিত্তিক residual plot বা Durbin–Watson দিয়ে। একটি plot-ই একসাথে কয়েকটি অনুমান যাচাই করে — তাই এটি diagnostics-এর প্রথম ও সবচেয়ে সস্তা পদক্ষেপ।

সমাধান ২ (★)¶

একই কথা নয়। তিনটি ধারণা স্পষ্টভাবে আলাদা:

leverage ($h_{ii}$): বিন্দুটি predictor-জায়গায় (অর্থাৎ $x$-অক্ষে) বাকিদের কেন্দ্র থেকে কত দূরে — কেবল $X$-এর ফাংশন, $y$-এর সাথে কোনো সম্পর্ক নেই। $h_{ii}=H_{ii}$, পরিসর $[0,1]$, গড় $p/n$।
outlier: বিন্দুটির residual বড় — অর্থাৎ মডেল-রেখা থেকে $y$-অক্ষে বহু দূরে।
influential: বিন্দুটি বাদ দিলে $\hat\beta$ (বা fitted মান) স্পষ্টভাবে বদলায় — Cook's $D$ এটাই পরিমাপ করে।

high-leverage কিন্তু small residual — হ্যাঁ, সম্ভব। একটি বিন্দু $x$-জায়গায় চরম প্রান্তিক (বড় $h_{ii}$) হতে পারে অথচ trend-রেখার ঠিক ওপরেই বসে — তখন residual ছোট, এমনকি বিন্দুটি fit-কে স্থিতিশীল করে। উল্টোভাবে $x$-কেন্দ্রের কাছের একটি বিন্দু বড় residual সত্ত্বেও fit-কে সামান্যই নাড়ায় (কম leverage)।

high-leverage হলেই fit নষ্ট করবে? না। leverage শুধু প্রভাবের সম্ভাবনা (potential) তৈরি করে। প্রকৃত ক্ষতি হয় তখনই যখন high leverage এবং বড় residual একসাথে থাকে — গাণিতিকভাবে Cook's $D_i=\dfrac{r_i^2}{p}\cdot\dfrac{h_{ii}}{1-h_{ii}}$-এ এই গুণফল দেখা যায়: একটি factor $r_i$ (residual), অন্যটি $\frac{h_{ii}}{1-h_{ii}}$ (leverage)। দুটোর একটিও শূন্যের কাছে হলে $D_i$ ছোট।

এক বাক্যে আলাদা: leverage = $x$-জায়গায় কত প্রান্তিক (সম্ভাবনা); outlier = $y$-অক্ষে কত অমিল (residual); influential = বাদ দিলে answer কত বদলায় (leverage × residual = Cook's $D$)।

সমাধান ৩ (★★)¶

এই তিন লক্ষণের একত্র উপস্থিতি multicollinearity-র পাঠ্যবই-স্বাক্ষর — এখানে area ও plot_size প্রায় একই তথ্য বহন করে ($\text{corr}(\text{area},\text{plot\_size})\approx 0.932$, কারণ plot_size $\approx 1.15\cdot\text{area}$)।

কেন individual coefficient "অর্থহীন/অস্থির" হয়: regression coefficient $\hat\beta_j$-এর অর্থ "অন্য সব predictor স্থির রেখে predictor $j$-এর একক পরিবর্তনে $\hat y$-এর পরিবর্তন"। কিন্তু plot_size প্রায় হুবহু area-এর rescaled কপি — area স্থির রেখে plot_size-কে এককভাবে নাড়ানোর মতো প্রায় কোনো ডেটাই নেই। তাই মডেল area ও plot_size-এর মধ্যে total effect-কে কীভাবে ভাগ করবে তা প্রায় নির্ধারণ করতে পারে না: যৌথ অবদান মোটামুটি স্থির থাকে, কিন্তু প্রতিটি coefficient আলাদাভাবে বিশাল অনিশ্চয়তা নিয়ে দোলে। ফলে se ফুলে যায় (area: $0.0258\to 0.0697$, প্রায় ২.৭ গুণ), area-এর coefficient $0.496\to 0.710$-এ লাফ দেয়, আর plot_size-এর coefficient আসে ঋণাত্মক ($-0.207$) — অথচ DGP-তে বড় plot_size মানেই বড় area, তাই দাম বাড়ারই কথা। এই উল্টো sign পুরোপুরি spurious। প্রমাণ ১২-তে এর সঠিক রূপ: $\operatorname{Var}(\hat\beta_j)\propto 1/(1-R_j^2)=\text{VIF}_j$, আর এখানে $R_{\text{area}}^2\approx 0.870$ বলে variance প্রায় $7.7$ গুণ ফোলা।

একটি সূক্ষ্ম ফাঁদ — significance ≠ অর্থপূর্ণতা। এখানে plot_size-এর $t=-3.29$ ($p=0.0012$) এখনও "significant" দেখায়, কিন্তু তার মান ও sign-ই অর্থহীন — তাই এই significance বিভ্রান্তিকর। (collinearity-র মাত্রা আরও তীব্র হলে — যেমন $\text{corr}\to 1$ — se এত বড় হয়ে যায় যে $t$-ও ধসে $0$-এর কাছে নামে।) মূল বার্তা একই: collinear predictor-এর per-coefficient সংখ্যা বিশ্বাসযোগ্য নয়।

কেন $R^2$ পড়ে না: $R^2$ মাপে মডেল সম্মিলিতভাবে $y$-এর কতটা variation ব্যাখ্যা করে। area ও plot_size একই তথ্য বহন করে, তাই দুজনে মিলে আগের মতোই variation ব্যাখ্যা করে — predictive শক্তি প্রায় অটুট ($R^2:\ 0.766\to 0.778$)। collinearity স্বতন্ত্র coefficient-এর ব্যাখ্যা/নির্ভুলতা নষ্ট করে, সামগ্রিক fit/পূর্বাভাস নয়।

সারকথা: multicollinearity হলো "ভাগাভাগির সংকট" — মডেল মোট কতটা জানে ($R^2$) তা জানে, কিন্তু কোন predictor-কে কত credit দেবে তা জানে না (unstable, sign-flipped coefficients)। সংখ্যায় ধরার যন্ত্র: VIF (সমাধান ৮)।

সমাধান ৪ (★★)¶

(ক) কে কোন null পরীক্ষা করে:

overall $F$-test: $H_0:\beta_1=\beta_2=\cdots=\beta_{p-1}=0$ (intercept বাদে সব slope একসাথে শূন্য) বনাম $H_1:$ অন্তত একটি $\ne 0$। অর্থাৎ "মডেলটা কি আদৌ predictor-শূন্য (intercept-only) মডেলের চেয়ে ভালো?" পরিসংখ্যান $F=\dfrac{\text{SSR}/(p-1)}{\text{SSE}/(n-p)}\sim F_{p-1,\,n-p}$ under $H_0$। আমাদের FULL মডেলে $F=213.7$ ($df=3,196$), $p\approx 1.6\times10^{-61}$।
প্রতিটি $t$-test: $H_0:\beta_j=0$ বনাম $\beta_j\ne 0$, অন্য সব predictor মডেলে থাকা অবস্থায় (partial/marginal contribution)। পরিসংখ্যান $t_j=\hat\beta_j/\widehat{\mathrm{se}}(\hat\beta_j)\sim t_{n-p}$।

মূল পার্থক্য: $F$ একটি যৌথ (joint) প্রশ্ন; প্রতিটি $t$ একটি শর্তাধীন স্বতন্ত্র প্রশ্ন।

(খ) "F significant অথচ কোনো individual $t$ নয়" — কখন ও কেন বিরোধ নয়: ঘটে ঠিক যখন predictor-রা পরস্পর strongly collinear। তখন— - যৌথভাবে তারা $y$-এর বড় অংশ ব্যাখ্যা করে ⇒ SSR বড় ⇒ F বিশাল ও significant; - কিন্তু কোনো একটি predictor-কে আলাদা করে "একমাত্র দায়ী" বলা যায় না (অন্যেরা তার তথ্য নকল করে) ⇒ প্রতিটি se ফোলা ⇒ প্রতিটি $t$ ছোট, হয়তো কোনোটিই significant নয়।

এটি স্ববিরোধ নয়, কারণ দুই test আলাদা প্রশ্ন করে: $F$ বলে "এই predictor-গুচ্ছ একত্রে কাজে লাগছে — হ্যাঁ"; প্রতিটি $t$ বলে "এই নির্দিষ্ট predictor-টি, বাকিরা থাকতে, একক-অপরিহার্য কিনা — না"। দুটোই একসাথে সত্য হতে পারে। বাস্তবে এটিই multicollinearity-র একটি লাল-পতাকা: $F$ চমৎকার অথচ সব $t$ দুর্বল দেখলে VIF যাচাই করুন। (আমাদের FULL মডেলে তিন predictor-ই অবশ্য আলাদাভাবেও significant — কারণ base মডেলে কোনো collinearity নেই; "F significant অথচ t নয়" পরিস্থিতি দেখা দেয় তখনই যখন plot_size-এর মতো নকল predictor ঢোকানো হয় ও দুই collinear predictor একসাথে রাখা হয়।)

সমাধান ৫ (★★)¶

দুই criterion-ই গঠন $-2\ell+(\text{penalty})\cdot K$, যেখানে $\ell$ maximized log-likelihood, $K$ free parameter-সংখ্যা (Gaussian regression-এ $K=p$ — statsmodels-এর প্রচলিত গণনা; কিছু পাঠ্যবইয়ে $\sigma^2$ আলাদা গুনে $K=p+1$ ধরা হয়, যা AIC/BIC-এ একটা ধ্রুব $+2$/$+\ln n$ shift দেয় কিন্তু মডেল-তুলনার $\Delta$ অপরিবর্তিত রাখে)। পার্থক্য কেবল penalty গুণাঙ্কে।

(ক) $n=200$-এ কোন penalty কড়া? AIC-তে প্রতি parameter $2$; BIC-তে $\ln n=\ln 200\approx 5.298$। যেহেতু $5.30>2$, BIC-এর penalty প্রায় আড়াই গুণ কড়া ($5.30/2\approx 2.65$)। সাধারণভাবে যখনই $n>e^2\approx 7.39$, BIC-এর জরিমানা AIC-র চেয়ে বড়।

(খ) কে বেশি parsimonious: কড়া জরিমানা মানে প্রতিটি অতিরিক্ত parameter "কেনা" আরও দামি — তাই BIC ছোট (কম-predictor) মডেলের দিকে বেশি ঝোঁকে। AIC তুলনায় উদার, বড় মডেল মেনে নিতে রাজি। ($n$ যত বড়, ব্যবধান তত প্রকট, কারণ $\ln n$ বাড়তে থাকে অথচ AIC-র $2$ স্থির।)

(গ) দার্শনিক লক্ষ্য: - BIC ⟶ "সত্য মডেল আছে, সেটি খুঁজি"। BIC consistent: যদি সত্য (data-জনক) মডেল বিবেচিত সেটে থাকে, $n\to\infty$-এ BIC সম্ভাবনা-$1$-এ সেটিকেই বাছে। (Bayesian posterior-এর approximation থেকে উদ্ভূত।) - AIC ⟶ "ভবিষ্যৎ-পূর্বাভাসে সেরাটি খুঁজি"। AIC asymptotically efficient: প্রত্যাশিত out-of-sample prediction error (KL-divergence) minimize করার লক্ষ্য — সত্য মডেল সেটে না-ও থাকতে পারে।

ব্যবহারিক টীকা: prediction-ভিত্তিক কাজে (ML-ধাঁচে) AIC; "কোন variable সত্যিই কারণ/গুরুত্বপূর্ণ" ধরনের explanatory কাজে BIC প্রায়ই পছন্দ। তবে দুটোই in-sample $\ell$-এর ওপর penalty — cross-validation এর আধুনিক, অনুমান-মুক্ত বিকল্প। (এবং — সমাধান ৯/১৪-এ দেখব — AIC/BIC কখনোই একা শেষ কথা নয়: collinearity তারা ধরে না।)

খ · গণনামূলক (computational)¶

সমাধান ৬ (★)¶

দেওয়া: $\hat\beta_{\text{area}}=0.4956$, $\widehat{\mathrm{se}}=0.02580$, $df=n-p=200-4=196$, $t^\*=t_{0.975,196}=1.972$।

(ক) $t$-statistic: $$ t=\frac{\hat\beta_{\text{area}}}{\widehat{\mathrm{se}}}=\frac{0.4956}{0.02580}\approx \mathbf{19.2}. $$ এত বড় $t$ মানে coefficient তার নিজের অনিশ্চয়তার ১৯ গুণ — অত্যন্ত নির্ভরযোগ্যভাবে শূন্য থেকে আলাদা।

(খ) $95\%$ CI: $$ 0.4956\pm 1.972\times 0.02580 = 0.4956\pm 0.05088 = (\mathbf{0.4447,\ 0.5465}). $$

(গ) সিদ্ধান্ত (CI–test duality): CI $(0.4447,0.5465)$-তে $0$ নেই (ধারেকাছেও নয়)। CI–test duality অনুসারে "$95\%$ CI-তে $\beta_{\text{area}}=0$ না থাকা" $\Leftrightarrow$ "$\alpha=0.05$-এ $H_0:\beta_{\text{area}}=0$ বাতিল" — অর্থাৎ $H_0$ জোরালোভাবে বাতিল (সমতুল্যভাবে $\lvert t\rvert=19.2\gg t^\*=1.972$, p-value $\approx 6\times 10^{-47}$)। ব্যাখ্যা: অন্য predictor (age, rooms) স্থির রেখে area $1$ m² বাড়লে price গড়ে $\approx 0.50$ লাখ টাকা বাড়ে, এবং এই প্রভাব পরিসংখ্যানগতভাবে অত্যন্ত নিশ্চিত। (লক্ষণীয়, estimate $0.4956$ true মান $0.45$-এর খুব কাছে।)

সমাধান ৭ (★★)¶

দেওয়া: $\hat\beta_{\text{rooms}}=6.498$, $\widehat{\mathrm{se}}=0.5032$, $df=196$, $t^\*=1.972$।

(ক) $t$-stat: $t=6.498/0.5032\approx \mathbf{12.9}$ (vs $H_0:\beta=0$) — অত্যন্ত significant।

(খ) $95\%$ CI: $6.498\pm 1.972\times 0.5032 = 6.498\pm 0.9923=(\mathbf{5.506,\ 7.491})$।

(গ) ব্যাখ্যা: area ও age স্থির রেখে একটি অতিরিক্ত room থাকলে বাড়ির দাম গড়ে $\approx 6.5$ লাখ টাকা বেশি।

(ঘ) true-মান $6$ পরীক্ষা: এবার $H_0:\beta_{\text{rooms}}=6$ — $$ t=\frac{\hat\beta-6}{\widehat{\mathrm{se}}}=\frac{6.498-6}{0.5032}=\frac{0.498}{0.5032}\approx \mathbf{+0.99}. $$ $\lvert t\rvert=0.99<1.972=t^\*$, তাই $H_0:\beta=6$ বাতিল হয় না — data true-মান $6$-এর সাথে সম্পূর্ণ সঙ্গতিপূর্ণ। লক্ষ করুন এটি (খ)-এর CI $(5.506,7.491)$-এ $6$ থাকার ঠিক সমার্থক (আবার CI–test duality)। পাঠ: estimate ($6.498$) true ($6$) থেকে একটু সরে, কিন্তু সেই সরণ sampling-অনিশ্চয়তার মধ্যেই — সসীম $n$ ও noise-এর প্রত্যাশিত ফল, পক্ষপাত নয়।

সমাধান ৮ (★★)¶

দেওয়া: plot_size-যুক্ত মডেলে area-কে {age, rooms, plot_size} দিয়ে regress করে $R_{\text{area}}^2=0.870$।

(ক) VIF: $$ \text{VIF}{\text{area}}=\frac{1}{1-R. $$}}^2}=\frac{1}{1-0.870}=\frac{1}{0.130}\approx \mathbf{7.71

(খ) se কত গুণ ফোলা: $\operatorname{Var}(\hat\beta_j)\propto\text{VIF}_j$, তাই se $\propto\sqrt{\text{VIF}}=\sqrt{7.71}\approx \mathbf{2.78}$ গুণ — orthogonal হলে যা হতো তার প্রায় ২.৮ গুণ বড় standard error। (যাচাই: collinear মডেলে area-এর se $0.0258\to 0.0697$, প্রায় $2.70$ গুণ — base মডেলে area-এর VIF ঠিক $1$ নয় ($\approx 1.01$) বলে অনুপাত $\sqrt{7.71/1.01}\approx 2.76$-ই প্রত্যাশিত; মিলে যায়।)

(গ) থ্রেশহোল্ড ও দোষী: সাধারণ থাম্ব-রুল VIF $>5$ "মনোযোগ দাও", $>10$ "গুরুতর"। এখানে $7.71$ — $5$ ও $10$-এর মাঝে, অর্থাৎ "সতর্কতা" জোনে; গুরুতর নয়, কিন্তু উপেক্ষণীয়ও নয়। দোষী plot_size: এটি area-এর সাথে $\text{corr}\approx 0.93$, কার্যত rescaled নকল predictor। প্রতিকার সহজ — plot_size বা area-এর একটিকে বাদ দিন; তখন বাকি সবার VIF নেমে $\approx 1$ হয় (age, rooms-এর VIF এমনিতেই $\approx 1.04$ ও $1.00$, কারণ তারা স্বাধীন)।

সমাধান ৯ (★★)¶

মনে রাখুন (অধ্যায় ফাইলের প্রশ্নমতে): এই প্রশ্নের SSE-জোড়া নিছক দৃষ্টান্তমূলক, এই অধ্যায়ের dataset নয় — penalty-যুক্তিটা বিচ্ছিন্নভাবে দেখানোর hypothetical। চলমান উদাহরণের আসল AIC/BIC-চিত্র (যেখানে plot_size SSE যথেষ্ট কমায়) দেখুন সমাধান ১৪-এ।

দেওয়া (hypothetical): BASE ($p=4$, $\text{SSE}=26778.5$) বনাম $+z$ ($p=5$, $\text{SSE}=26776.9$, প্রায় অভিন্ন), $n=200$। এখানে $K=p$, তাই $K_{\text{BASE}}=4$, $K_{+}=5$, $\Delta K=1$।

Gaussian MLE log-likelihood $\ell=-\tfrac n2[\ln(2\pi)+\ln(\text{SSE}/n)+1]$। SSE প্রায় অভিন্ন বলে $\ell$ও প্রায় সমান, তাই $\Delta(-2\ell)\approx 0$।

(ক) $\Delta\text{AIC}$ ও $\Delta\text{BIC}$ ($+z$ − BASE), penalty-চালিত: $$ \Delta\text{AIC}=\underbrace{\Delta(-2\ell)}_{\approx 0}+2\,\Delta K\approx 0+2(1)=\mathbf{+2.0}, $$ $$ \Delta\text{BIC}=\Delta(-2\ell)+(\ln 200)\,\Delta K\approx 0+5.298(1)=\mathbf{+5.3}. $$

(খ) কে জেতে — ও সাধারণ নীতি: এই hypothetical-এ দুটোই ধনাত্মক — অর্থাৎ SSE না কমিয়ে $z$ যোগ করলে AIC ও BIC দুটোই বাড়ে (ছোটটাই ভালো, তাই বাড়া = খারাপ)। সুতরাং উভয় criterion-ই BASE মডেলকে বাছে, $z$-কে নাকচ করে; BIC বেশি জোরালোভাবে ($+5.3$ vs $+2.0$), কারণ তার per-parameter জরিমানা কড়া। সাধারণ নীতি: AIC/BIC একটি predictor কেবল তখনই রাখে যখন তার SSE-হ্রাস penalty-বৃদ্ধিকে পুষিয়ে দেয়।

সতর্কতা — চলমান উদাহরণে কিন্তু উল্টোটা। এই অধ্যায়ের আসল dataset-এ plot_size যোগ করলে SSE যথেষ্ট কমে ($52092\to 49347$, partial $F=10.85$, $p=0.0012$), তাই সেখানে AIC/BIC দুটোই নামে এবং mechanically plot_size রাখার পক্ষে যায় (সমাধান ১৪)। অর্থাৎ plot_size বাদ দেওয়ার যুক্তি AIC/BIC থেকে আসে না — আসে VIF + se-inflation + অস্থির/উল্টো-sign coefficient থেকে। এই প্রশ্নের hypothetical-টি কেবল "SSE প্রায় সমান হলে penalty-ই সিদ্ধান্ত নেয়" — এই বিচ্ছিন্ন নীতিটা বোঝায়।

সমাধান ১০ (★★)¶

দেওয়া: row 0-এর জন্য $h_{00}=0.1328$, $p=4$, $n=200$, এই বিন্দু বাদ দিলে intercept $16.95\to 22.0$, $D_0=0.529$ (পরের সর্বোচ্চ $D\approx 0.047$, threshold $4/n=0.02$)।

(ক) গড় leverage: $\bar h=p/n=4/200=\mathbf{0.02}$।

(খ) কত গুণ: $h_{00}/\bar h=0.1328/0.02\approx \mathbf{6.6}$ গুণ গড়ের ওপরে।

(গ) থ্রেশহোল্ড: $2p/n=2(0.02)=0.04$; $3p/n=0.06$। $h_{00}=0.1328$ দুটোকেই বহুগুণে ছাড়িয়েছে ($0.1328\gg 0.06$) — নিঃসন্দেহে high-leverage বিন্দু। (কারণ স্পষ্ট: area=340 বাকি area-পরিসর $[50,200]$-এর বহু বাইরে।)

(ঘ) influential কিনা — একটি সূক্ষ্ম পাঠ। Cook's $D_0=0.529$ কঠোর $D>1$ rule ছোঁয় না — তাই "নিঃসন্দেহে influential" তকমা $D>1$ মাপকাঠিতে দেওয়া যায় না। কিন্তু এটি আরও দুই দিক থেকে স্পষ্টভাবে প্রভাবশালী: (i) সচরাচর-ব্যবহৃত $4/n=0.02$ threshold বহুগুণে ছাড়ায় এবং $0.5$-ও পেরোয়; (ii) পরের সর্বোচ্চ $D$ (মাত্র $0.047$)-এর প্রায় ১১ গুণ — অর্থাৎ আপেক্ষিক মাত্রায় এটি অনন্য, dataset-এর সবচেয়ে প্রভাবশালী বিন্দু। প্রমাণ: একা এই বিন্দু intercept-কে $16.95$ থেকে $22.0$-এ টেনেছে (true মান $25$-এর কাছে চলে আসে!), area-slope-ও সামান্য বদলেছে ($0.4956\to 0.4600$), residual SE $16.30\to 15.76$-এ নেমেছে। সঠিক practice: বিন্দুটি নীরবে বাদ না দিয়ে আগে কারণ অনুসন্ধান (ডেটা-এন্ট্রি ভুল? সত্যিই ব্যতিক্রমী বাড়ি?); তারপর বিন্দুটি নিয়ে ও ছাড়া — দুটি fit রিপোর্ট করা, যাতে পাঠক প্রভাবের মাত্রা দেখতে পান। (পাঠ: কেবল একটি কঠোর threshold-এ আটকে না থেকে আপেক্ষিক মাত্রা দেখাই বেশি নির্ভরযোগ্য — $D>1$ পার না করেও একটি বিন্দু কার্যত influential হতে পারে।)

গ · প্রমাণভিত্তিক (proof-based)¶

সমাধান ১১ (★★)¶

$H=X(X^\top X)^{-1}X^\top$, $X\in\mathbb R^{n\times p}$ full-rank। সংক্ষেপে $M=(X^\top X)^{-1}$ (symmetric, কারণ $X^\top X$ symmetric ⇒ তার inverse-ও)।

(ক) Symmetric: $$ H^\top=\big(XMX^\top\big)^\top=(X^\top)^\top M^\top X^\top=X M X^\top=H, $$ যেহেতু $M^\top=M$। ∎

(খ) Idempotent: $$ H^2=XMX^\top\,XMX^\top=XM\,(X^\top X)\,MX^\top=XM\,(X^\top X)(X^\top X)^{-1}\,X^\top=XMX^\top=H, $$ কারণ $M(X^\top X)=(X^\top X)^{-1}(X^\top X)=I_p$। ∎ (জ্যামিতিক অর্থ: $H$ হলো column-space $\mathcal C(X)$-এর ওপর orthogonal projection — দুবার project করা একবারের সমান।)

(গ) $0\le h_{ii}\le 1$: $h_{ii}=e_i^\top H e_i$ ($e_i$ = $i$-তম standard basis vector)। symmetric+idempotent ⇒ $H=H^\top H$, তাই $$ h_{ii}=e_i^\top H^\top H e_i=\lVert He_i\rVert^2\ge 0. $$ উপরের সীমা: row $i$-তে $h_{ii}=\sum_{j}H_{ij}^2$ (idempotent+symmetric থেকে $H_{ii}=\sum_j H_{ij}H_{ji}=\sum_j H_{ij}^2$), যা $=h_{ii}^2+\sum_{j\ne i}H_{ij}^2$। অতএব $h_{ii}-h_{ii}^2=\sum_{j\ne i}H_{ij}^2\ge 0\Rightarrow h_{ii}(1-h_{ii})\ge 0\Rightarrow h_{ii}\le 1$। ∎

(ঘ) $\sum h_{ii}=\operatorname{tr}(H)=p$: trace-এর cyclic ধর্ম ($\operatorname{tr}(ABC)=\operatorname{tr}(CAB)$, ০.৫): $$ \operatorname{tr}(H)=\operatorname{tr}!\big(X(X^\top X)^{-1}X^\top\big)=\operatorname{tr}!\big((X^\top X)^{-1}X^\top X\big)=\operatorname{tr}(I_p)=p. $$ সুতরাং গড় leverage $\bar h=\frac1n\sum h_{ii}=p/n$ — এটাই "$2p/n$ থ্রেশহোল্ড" ও সমাধান ১০-এর ভিত্তি। ∎

সমাধান ১২ (★★★)¶

LINE: $y=X\beta+\varepsilon$, $\mathbb E[\varepsilon]=0$, $\operatorname{Var}(\varepsilon)=\sigma^2 I_n$, $X$ স্থির full-rank।

(ক) $\operatorname{Var}(\hat\beta)=\sigma^2(X^\top X)^{-1}$: $\hat\beta=(X^\top X)^{-1}X^\top y=Ay$ যেখানে $A=(X^\top X)^{-1}X^\top$ একটি স্থির matrix। স্থির $A$-র জন্য $\operatorname{Var}(Ay)=A\operatorname{Var}(y)A^\top$, আর $\operatorname{Var}(y)=\operatorname{Var}(\varepsilon)=\sigma^2 I$: $$ \operatorname{Var}(\hat\beta)=A(\sigma^2 I)A^\top=\sigma^2\,(X^\top X)^{-1}X^\top\,X(X^\top X)^{-1}=\sigma^2(X^\top X)^{-1}, $$ যেহেতু $(X^\top X)^{-1}X^\top X=I$ এবং $\big((X^\top X)^{-1}\big)^\top=(X^\top X)^{-1}$। ∎

(খ) $\operatorname{Var}(\hat\beta_j)\propto 1/(1-R_j^2)=\text{VIF}_j$: predictor-গুলো centered ধরি। $(X^\top X)^{-1}$-এর $j$-তম কর্ণ-উপাদান বের করতে block-inverse/partitioned-matrix সূত্র লাগে। predictor $j$-কে বাকি predictor-দের ওপর regress করলে যে residual sum of squares থাকে তা $\text{SS}_j(1-R_j^2)$, যেখানে $\text{SS}_j=\sum_i(x_{ij}-\bar x_j)^2$ এবং $R_j^2$ ওই auxiliary regression-এর $R^2$। partitioned-inverse সূত্র দেয় $$ (X^\top X)^{-1}{jj}=\frac{1}{\text{SS}_j\,(1-R_j^2)}. $$ অতএব $$ \operatorname{Var}(\hat\beta_j)=\sigma^2(X^\top X)^{-1}}=\frac{\sigma^2}{\text{SSj}\cdot\underbrace{\frac{1}{1-R_j^2}}. $$ প্রথম factor }_j$\sigma^2/\text{SS}_j$ হলো predictor $j$ যদি সম্পূর্ণ orthogonal হতো তখনকার variance; দ্বিতীয় factor $\text{VIF}_j$ হলো collinearity-জনিত গুণক বৃদ্ধি। $R_j^2\to 1$ (predictor $j$ অন্যদের দিয়ে প্রায় সম্পূর্ণ ব্যাখ্যাযোগ্য) ⇒ VIF $\to\infty$ ⇒ variance বিস্ফোরিত। ∎

(গ) নাম-ব্যাখ্যা: $\widehat{\mathrm{se}}(\hat\beta_j)=\hat\sigma\sqrt{(X^\top X)^{-1}_{jj}}=\dfrac{\hat\sigma}{\sqrt{\text{SS}_j}}\sqrt{\text{VIF}_j}$। অর্থাৎ VIF ঠিক ততগুণ যত গুণে coefficient-এর variance বেড়েছে orthogonal-আদর্শের তুলনায় — তাই literal নাম "variance inflation factor", আর se বাড়ে $\sqrt{\text{VIF}_j}$ গুণ। (আমাদের উদাহরণে $\text{VIF}_{\text{area}}\approx 7.71$, তাই se $\sqrt{7.71}\approx 2.78$ গুণ ফোলা — সমাধান ৮।) ∎

সমাধান ১৩ (★★★)¶

(ক) $R^2$ অ-হ্রাসমান: ধরা যাক মডেল $\mathcal M_1$-এ predictor-সেট $S$, আর $\mathcal M_2$-এ $S\cup\{z\}$ (একটি predictor বেশি)। OLS করে $\text{SSE}=\min_\beta\lVert y-X\beta\rVert^2$, অর্থাৎ fitted space-এ minimum। $\mathcal M_2$-এর প্যারামিটার-জগৎ $\mathcal M_1$-কে ধারণ করে (নতুন coefficient $=0$ বসালেই $\mathcal M_1$ পুনরুৎপাদিত)। বৃহত্তর সেটে minimization-এর minimum কখনো বড় হতে পারে না, তাই $\text{SSE}_2\le \text{SSE}_1$। যেহেতু SST অপরিবর্তিত (একই $y$), $R^2=1-\text{SSE}/\text{SST}$ অ-হ্রাসমান — predictor যোগ করলে কখনো কমে না (সাধারণত একটু বাড়ে, এমনকি predictor-টি random noise হলেও)। ∎

(খ) $R^2_{\text{adj}}$ কমতে পারে: লিখি $$ R^2_{\text{adj}}=1-\frac{\text{SSE}/(n-p)}{\text{SST}/(n-1)}=1-\frac{\hat\sigma^2}{\text{SST}/(n-1)},\qquad \hat\sigma^2=\frac{\text{SSE}}{n-p}. $$ denominator $\text{SST}/(n-1)$ predictor-সংখ্যায় স্থির। তাই $R^2_{\text{adj}}$ বাড়বে/কমবে শুধু $\hat\sigma^2=\text{SSE}/(n-p)$ কমলে/বাড়লে। predictor যোগ করলে দুটি বিপরীত প্রভাব: - numerator $\text{SSE}$ একটু কমে (≤, ক-অংশ থেকে); - denominator $n-p$ও কমে ($p$ ১ বাড়ে)।

যদি অতিরিক্ত predictor-টি প্রায় অর্থহীন হয়, SSE-র হ্রাস নগণ্য অথচ $n-p$-এর হ্রাস বাস্তব — তখন অনুপাত $\text{SSE}/(n-p)$ আসলে বাড়ে, তাই $R^2_{\text{adj}}$ কমে। অর্থাৎ adjusted $R^2$ কেবল তখনই বাড়ে যখন নতুন predictor "তার df-খরচ পুষিয়ে দেওয়ার মতো যথেষ্ট" SSE কমায় ($t$-statistic-এর শর্তে: মোটামুটি $\lvert t\rvert>1$ হলে $R^2_{\text{adj}}$ বাড়ে)। ∎

(গ) দৃষ্টান্তে যাচাই (illustrative, এই অধ্যায়ের dataset নয়): ধরা যাক $\text{SST}=222476$, $n=200$, BASE-এ ($p=4$) $\text{SSE}=52092$, আর একটি অর্থহীন predictor যোগে ($p=5$) $\text{SSE}=52060$ (হ্রাস মাত্র $\approx 32$, df-খরচের তুলনায় নগণ্য)। তাহলে $$ \hat\sigma^2:\ \frac{52092}{196}=265.78\ \longrightarrow\ \frac{52060}{195}=266.97\ (\textbf{বেড়েছে}), $$ ফলে $R^2_{\text{adj}}=1-\hat\sigma^2/(\text{SST}/199)$ কমে ($\approx 0.7623\to 0.7613$) — ঠিক যেমন (খ) বলে। নগণ্য SSE-হ্রাস df-হ্রাসকে পোষাতে পারল না।

চলমান উদাহরণে সতর্কতা: এই অধ্যায়ের আসল dataset-এ plot_size কিন্তু SSE যথেষ্ট কমায় ($52092\to 49347$), তাই সেখানে $\hat\sigma^2$ $265.78\to 253.06$-এ নামে এবং $R^2_{\text{adj}}$ সামান্য বাড়ে ($0.7623\to 0.7736$)। অর্থাৎ adjusted $R^2$ একা plot_size-কে "অর্থহীন" বলে না — উল্টো সামান্য সমর্থন করে। plot_size-এর সমস্যা (collinearity) ধরা পড়ে VIF + se-inflation + sign-flip-এ, adjusted $R^2$ বা AIC/BIC-তে নয় (সমাধান ১৪, §৫.E4)। এটাই দেখায় কেন একটিমাত্র মাপকাঠির ওপর নির্ভর করা বিপজ্জনক। ∎

ঘ · কোডিং (Python)¶

সমাধান ১৪ (★★★)¶

নিচের একটি self-contained script সব অংশ (১–৪) করে; from-scratch ফল statsmodels-এর সাথে মিলিয়ে যাচাই করা হয়েছে।

import numpy as np
from scipy import stats

rng = np.random.default_rng(20260619)
n = 200
area      = rng.uniform(50, 200, n)
age       = rng.uniform(0, 40, n)
rooms     = rng.integers(1, 9, n)              # 1..8 inclusive
price     = (0.45*area - 0.8*age + 6*rooms + 25 + rng.normal(0, 15, n))
plot_size = 1.15*area + rng.normal(0, 10, n)   # collinear with area (corr ~ 0.93)

# high-leverage / influential point at ROW 0  (injected into the live arrays)
area[0], price[0] = 340.0, 250.0

# ---------- from-scratch OLS + diagnostics ----------
def fit_diag(X, y):
    n, p = X.shape
    XtX_inv = np.linalg.inv(X.T @ X)
    beta = XtX_inv @ X.T @ y                    # = np.linalg.solve(X.T@X, X.T@y)
    yhat = X @ beta
    resid = y - yhat
    SSE = resid @ resid
    SST = np.sum((y - y.mean())**2)
    R2  = 1 - SSE/SST
    R2adj = 1 - (SSE/(n-p)) / (SST/(n-1))
    s2  = SSE/(n-p)                             # unbiased sigma^2
    se  = np.sqrt(np.diag(s2 * XtX_inv))        # se(beta_j) = sigma_hat * sqrt((X'X)^-1_jj)
    t   = beta / se
    H   = X @ XtX_inv @ X.T
    h   = np.diag(H)                            # leverages
    stud = resid / np.sqrt(s2 * (1 - h))        # (internally) studentized residuals
    cook = (stud**2 / p) * (h / (1 - h))        # Cook's distance
    return dict(beta=beta, se=se, t=t, R2=R2, R2adj=R2adj, s2=s2,
                SSE=SSE, SST=SST, h=h, cook=cook, resid=resid, yhat=yhat, n=n, p=p)

# ===== 1. BASE model: price ~ area + age + rooms  (leverage point is IN the data) =====
Xb = np.column_stack([np.ones(n), area, age, rooms])
b = fit_diag(Xb, price)
labels = ["intercept", "area", "age", "rooms"]
df = b['n'] - b['p']                            # 196
tcrit = stats.t.ppf(0.975, df)                  # ~1.972
print("=== BASE: coefficient inference ===")
for j, lab in enumerate(labels):
    lo = b['beta'][j] - tcrit*b['se'][j]
    hi = b['beta'][j] + tcrit*b['se'][j]
    print(f"{lab:9s} beta={b['beta'][j]:8.4f}  se={b['se'][j]:7.4f}  "
          f"t={b['t'][j]:8.3f}  95%CI=({lo:8.4f},{hi:8.4f})")
print(f"R2={b['R2']:.4f}  R2adj={b['R2adj']:.4f}  sigma_hat={np.sqrt(b['s2']):.3f}")
# overall F-test
SSR, p_pred = b['SST'] - b['SSE'], b['p'] - 1
F = (SSR/p_pred) / (b['SSE']/df)
print(f"overall F = {F:.2f}  (df={p_pred},{df})  p = {stats.f.sf(F, p_pred, df):.2e}")

# ===== 2. Diagnostics: leverage + Cook's D (row 0 is the injected point) =====
i_h = np.argmax(b['h']); i_c = np.argmax(b['cook'])
print("\n=== diagnostics (leverage point at row 0) ===")
print(f"max leverage  : idx {i_h}, h={b['h'][i_h]:.4f}  "
      f"(mean h = p/n = {b['p']/n:.4f}, 2p/n = {2*b['p']/n:.4f})")
print(f"max Cook's D   : idx {i_c}, D={b['cook'][i_c]:.4f}  "
      f"(D>1? {'yes' if b['cook'][i_c] > 1 else 'NO'}; 4/n={4/n:.3f}; "
      f"next largest = {np.sort(b['cook'])[-2]:.4f})")
# refit WITHOUT row 0
keep = np.arange(n) != 0
b0 = fit_diag(Xb[keep], price[keep])
print("with row 0   : beta =", np.round(b['beta'], 4),
      f" R2={b['R2']:.4f}  sigma={np.sqrt(b['s2']):.3f}")
print("without row 0: beta =", np.round(b0['beta'], 4),
      f" R2={b0['R2']:.4f}  sigma={np.sqrt(b0['s2']):.3f}")

# ===== 3. Multicollinearity: VIF in +plot_size model =====
preds = {"area": area, "age": age, "rooms": rooms.astype(float), "plot_size": plot_size}
P = np.column_stack(list(preds.values()))
print("\n=== VIF (+plot_size model) ===")
for j, name in enumerate(preds):
    others = np.delete(P, j, axis=1)
    A = np.column_stack([np.ones(n), others])
    coef, *_ = np.linalg.lstsq(A, P[:, j], rcond=None)
    R2j = 1 - np.sum((P[:, j] - A@coef)**2) / np.sum((P[:, j] - P[:, j].mean())**2)
    vif = 1/(1-R2j)
    zone = "<-- >10 severe" if vif > 10 else ("<-- 5-10 warning" if vif > 5 else "")
    print(f"  {name:10s} VIF = {vif:7.2f}   (R2_j = {R2j:.4f})  {zone}")
print(f"  corr(area, plot_size) = {np.corrcoef(area, plot_size)[0,1]:.4f}")

# ===== 4. Model selection: BASE vs +plot_size =====
Xc = np.column_stack([np.ones(n), area, age, rooms, plot_size])
c = fit_diag(Xc, price)
def aic_bic(SSE, n, p):
    # statsmodels OLS convention: K = p (intercept + slopes). ll = Gaussian MLE
    # log-likelihood. (Counting sigma^2 too -> K=p+1 shifts both by a constant
    # +2 / +ln n that cancels in the DELTA; statsmodels uses K=p, matched here.)
    ll = -0.5*n*(np.log(2*np.pi) + np.log(SSE/n) + 1)
    K = p
    return -2*ll + 2*K, -2*ll + np.log(n)*K
aic_b, bic_b = aic_bic(b['SSE'], n, b['p'])
aic_c, bic_c = aic_bic(c['SSE'], n, c['p'])
print("\n=== selection: BASE vs +plot_size ===")
print(f"BASE     : R2adj={b['R2adj']:.4f}  SSE={b['SSE']:.1f}  AIC={aic_b:8.2f}  BIC={bic_b:8.2f}")
print(f"+plot_sz : R2adj={c['R2adj']:.4f}  SSE={c['SSE']:.1f}  AIC={aic_c:8.2f}  BIC={bic_c:8.2f}")
print(f"dAIC={aic_c-aic_b:+.2f}  dBIC={bic_c-bic_b:+.2f}  "
      f"-> AIC/BIC LOWER with plot_size (mechanically keep it)")
# why we still drop plot_size:
print(f"plot_size coef = {c['beta'][4]:+.4f}  se = {c['se'][4]:.4f}  "
      f"t = {c['t'][4]:.2f}  (sign FLIPPED vs true +1.15 -> spurious)")
print(f"se(area): base {b['se'][1]:.4f} -> +plot {c['se'][1]:.4f}  "
      f"({c['se'][1]/b['se'][1]:.2f}x inflated)")
# partial F for plot_size
Fp = ((b['SSE'] - c['SSE'])/1) / (c['SSE']/(n - c['p']))
print(f"partial F (plot_size) = {Fp:.3f}  p = {stats.f.sf(Fp, 1, n-c['p']):.4f}")

# ---------- cross-check with statsmodels (authoritative AIC/BIC) ----------
try:
    import statsmodels.api as sm
    m  = sm.OLS(price, Xb).fit()
    mc = sm.OLS(price, Xc).fit()
    assert np.allclose(m.params, b['beta']) and np.allclose(m.bse, b['se'])
    infl = m.get_influence()
    assert np.allclose(infl.hat_matrix_diag, b['h'])
    assert np.allclose(infl.cooks_distance[0], b['cook'])
    print("\n[statsmodels cross-check passed: params, se, leverage, Cook's D all match]")
    print(f"  sm BASE     AIC={m.aic:.2f}  BIC={m.bic:.2f}")
    print(f"  sm +plot_sz AIC={mc.aic:.2f}  BIC={mc.bic:.2f}  "
          f"(dAIC={mc.aic-m.aic:+.2f}, dBIC={mc.bic-m.bic:+.2f})")
except ImportError:
    print("\n[statsmodels not installed -- skipped cross-check]")

প্রত্যাশিত আউটপুট (মূল সংখ্যা):

=== BASE: coefficient inference ===
intercept beta= 16.9547  se= 4.6712  t=   3.630  95%CI=(  7.7424, 26.1670)
area      beta=  0.4956  se= 0.0258  t=  19.208  95%CI=(  0.4447,  0.5465)
age       beta= -0.7203  se= 0.1009  t=  -7.138  95%CI=( -0.9193, -0.5213)
rooms     beta=  6.4982  se= 0.5032  t=  12.915  95%CI=(  5.5059,  7.4905)
R2=0.7659  R2adj=0.7623  sigma_hat=16.303
overall F = 213.69  (df=3,196)  p = 1.59e-61

=== diagnostics (leverage point at row 0) ===
max leverage  : idx 0, h=0.1328  (mean h = p/n = 0.0200, 2p/n = 0.0400)
max Cook's D   : idx 0, D=0.5289  (D>1? NO; 4/n=0.020; next largest = 0.0473)
with row 0   : beta = [16.9547  0.4956 -0.7203  6.4982]  R2=0.7659  sigma=16.303
without row 0: beta = [21.9973  0.46   -0.7668  6.4966]  R2=0.7561  sigma=15.758

=== VIF (+plot_size model) ===
  area       VIF =    7.71   (R2_j = 0.8703)  <-- 5-10 warning
  age        VIF =    1.04   (R2_j = 0.0371)
  rooms      VIF =    1.00   (R2_j = 0.0043)
  plot_size  VIF =    7.80   (R2_j = 0.8717)  <-- 5-10 warning
  corr(area, plot_size) = 0.9315

=== selection: BASE vs +plot_size ===
BASE     : R2adj=0.7623  SSE=52092.4  AIC= 1688.07  BIC= 1701.26
+plot_sz : R2adj=0.7736  SSE=49347.4  AIC= 1679.24  BIC= 1695.73
dAIC=-8.83  dBIC=-5.53  -> AIC/BIC LOWER with plot_size (mechanically keep it)
plot_size coef = -0.2065  se = 0.0627  t = -3.29  (sign FLIPPED vs true +1.15 -> spurious)
se(area): base 0.0258 -> +plot 0.0697  (2.70x inflated)
partial F (plot_size) = 10.847  p = 0.0012

[statsmodels cross-check passed: params, se, leverage, Cook's D all match]
  sm BASE     AIC=1688.07  BIC=1701.26
  sm +plot_sz AIC=1679.24  BIC=1695.73  (dAIC=-8.83, dBIC=-5.53)

AIC/BIC convention নোট। ওপরের aic_bic statsmodels-এর মতো $K=p$ (intercept + slopes) ধরে, তাই from-scratch ও statsmodels মান হুবহু মেলে (BASE $1688.07/1701.26$, +plot_size $1679.24/1695.73$ — cross-check assert-এ নিশ্চিত)। কিছু পাঠ্যবই $\sigma^2$-কে আলাদা parameter গুনে $K=p+1$ ধরে; তাতে প্রতিটি AIC একটা ধ্রুব $+2$ (BIC $+\ln n$) shift পায়, কিন্তু সেটা সব মডেলে এক বলে মডেল-তুলনার $\Delta\text{AIC}/\Delta\text{BIC}$ — এবং তাই চূড়ান্ত সিদ্ধান্ত — অপরিবর্তিত থাকে।

ব্যাখ্যা ও পাঠ (অংশ-ভিত্তিক):

Inference: তিন predictor-ই অত্যন্ত significant ($\lvert t\rvert>7$, CI-তে $0$ নেই)। overall $F=213.7$ ($p\approx1.6\times10^{-61}$) — মডেল intercept-only-র চেয়ে বিপুলভাবে ভালো। $\hat\beta=[16.95,0.4956,-0.7203,6.498]$ true $[25,0.45,-0.8,6]$-এর কাছাকাছি (intercept সবচেয়ে কম-নিশ্চিত — তার se বড়, $4.67$ — এবং row 0 leverage point তাকে নিচে টেনে রেখেছে)।
Diagnostics: row 0-ই একসাথে সর্বোচ্চ leverage ($h=0.1328$, গড়ের ৬.৬ গুণ, $2p/n$-র বহু ওপরে) ও সর্বোচ্চ Cook's $D$ ($0.529$) — অর্থাৎ ইচ্ছাকৃত high-leverage point নিখুঁতভাবে ধরা পড়ল। লক্ষণীয় — $D_0=0.529$ কঠোর $D>1$ rule পার করে না, কিন্তু $4/n=0.02$ threshold বহুগুণে ছাড়ায় এবং পরের সর্বোচ্চ $D$ ($0.047$)-এর ~১১ গুণ; তাই আপেক্ষিক মাত্রায় এটি স্পষ্টভাবে সবচেয়ে প্রভাবশালী। এই বিন্দু বাদ দিলে intercept $16.95\to 22.0$ (true $25$-এর কাছে), area-slope $0.4956\to 0.4600$, $\hat\sigma$ $16.30\to 15.76$ — একটিমাত্র বিন্দু গোটা fit-কে DGP থেকে সরিয়ে রেখেছিল।
Multicollinearity: area ও plot_size-এর VIF $\approx 7.7$ — "$5$–$10$ সতর্কতা" জোনে, "$>10$ গুরুতর" নয়; age, rooms-এর VIF $\approx 1$ (স্বাধীন বলে)। corr$(\text{area},\text{plot\_size})\approx 0.93$। সংকেত স্পষ্ট — plot_size area-এর সাথে redundant।
Selection — এবং একটি গুরুত্বপূর্ণ মোড়। plot_size যোগে SSE যথেষ্ট কমে ($52092\to 49347$; partial $F=10.85$, $p=0.0012$), তাই $R^2_{\text{adj}}$ সামান্য বাড়ে ($0.7623\to 0.7736$) এবং AIC ($-8.83$), BIC ($-5.53$) দুটোই নামে — অর্থাৎ AIC/BIC/adjusted-$R^2$ যান্ত্রিকভাবে plot_size রাখার পক্ষে! তবু আমরা plot_size বাদ দিই, কারণ (i) area↔plot_size-এর VIF $\approx 7.7$ (collinearity), (ii) area-এর se $0.0258\to 0.0697$-এ (২.৭ গুণ) ফুলে যায়, এবং (iii) plot_size-এর coefficient আসে ঋণাত্মক ($-0.207$) — অথচ DGP-তে তার প্রভাব ধনাত্মক ($+1.15$ scale), তাই sign পুরোপুরি spurious। মূল শিক্ষা: AIC/BIC কম মানেই multicollinearity নেই, এমন নয় — এরা কেবল fit-vs-complexity দেখে, coefficient-এর ব্যাখ্যাযোগ্যতা বা stability নয়। ব্যাখ্যাযোগ্য, স্থিতিশীল মডেল চাইলে M3 (BASE)-ই কাম্য; model selection কখনো একটিমাত্র সংখ্যার দাস হওয়া উচিত নয় — VIF, se-stability, sign ও domain-জ্ঞান একসাথে বিচার করতে হয়। আর row 0 leverage point দেখায়, diagnostics উপেক্ষা করলে একটিমাত্র বিন্দুই উপসংহার বদলে দিতে পারে।

প্লট কোড (ঐচ্ছিক, §২-এর দুই চিত্র):

import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
ax[0].scatter(b['yhat'], b['resid'], s=14)
ax[0].axhline(0, color='k', lw=1)
ax[0].set_xlabel("fitted values"); ax[0].set_ylabel("residuals")
ax[0].set_title("Residual vs Fitted (row 0 stands apart)")
ax[1].stem(np.arange(n), b['cook'])
ax[1].axhline(4/n, color='r', ls='--', label="D = 4/n threshold")
ax[1].set_xlabel("observation index"); ax[1].set_ylabel("Cook's distance")
ax[1].set_title("Cook's distance (row 0 dominates)"); ax[1].legend()
plt.tight_layout(); plt.show()

মূল take-away: একটি সম্পূর্ণ regression-রিপোর্ট কখনো শুধু $\hat\beta$ ও $R^2$ নয় — তার সাথে চাই (i) inference (se, $t$, CI, $F$: প্রতিটি দাবি কতটা নিশ্চিত), (ii) diagnostics (residual plot, leverage, Cook's $D$: অনুমান ও influential point), এবং (iii) selection (VIF, $R^2_{\text{adj}}$, AIC/BIC: কোন predictor রাখব — কিন্তু কোনো একটিকে অন্ধভাবে নয়)। এই pipeline-ই ৫.১-এর "মডেল বানানো"-কে বিশ্বাসযোগ্য বিজ্ঞানে পরিণত ক