সমাধান — অধ্যায় ৫.৬ · Mixed-Effects / Hierarchical Models¶

অধ্যায় part-5-modeling/05-06-mixed-effects-hierarchical-models.md-এর section ৭-এর সব অনুশীলনীর পূর্ণ সমাধান। চলমান dataset — exam scores, seed np.random.default_rng(20260619), $J=25$ স্কুল, unbalanced ($n_j$ ১০–৩০, মোট $N=533$); model $y_{ij}=50+2\cdot\text{hours}_{ij}+u_j+\varepsilon_{ij}$, $u_j\sim\mathcal N(0,6^2)$, $\varepsilon_{ij}\sim\mathcal N(0,8^2)$। প্রতিটি সংখ্যাগত ফল statsmodels (MixedLM, REML) দিয়ে যাচাই করা হয়েছে — §ঘ-এর runnable script এ-ই data-generation-ক্রমে নিচের canonical সংখ্যা হুবহু পুনরুৎপাদন করে। canonical সংখ্যা: pooled OLS $\hat\beta_1=1.800$ (SE $0.149$, intercept SE $0.867$); mixed $\hat\beta_0=50.82$ (SE $1.416$), $\hat\beta_1=1.878$ (SE $0.120$); $\hat\sigma_u=6.14$ ($\hat\sigma_u^2=37.64$), $\hat\sigma_\varepsilon=7.98$ ($\hat\sigma_\varepsilon^2=63.69$); ICC $\rho=0.371$; logLik $=-1894.5$; স্কুল $0$ ($n_0=20$) BLUP $\hat u_0=+6.05$।

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

৭.১ ★ — fixed effect বনাম random effect¶

(ক) hours-এর coefficient (সহগ) $\beta_1$ একটা fixed effect: পড়ার ঘণ্টা ও স্কোরের সম্পর্ক একটা population-জোড়া স্থির পরিমাণ যা সব স্কুলে অভিন্ন বলে ধরা হয় — আমরা এই একটিমাত্র সংখ্যা $\beta_1$-তেই আগ্রহী। school-এর প্রভাব random effect $u_j$: এই নির্দিষ্ট ২৫টা স্কুল নিজে আমাদের লক্ষ্য নয়, বরং এরা একটা বৃহত্তর জনগোষ্ঠী থেকে নমুনা — আমরা সেই population-এর বৈচিত্র্য ($\sigma_u^2$) সম্পর্কে অনুমান করতে চাই, প্রতিটা স্কুলের আলাদা মান নয়।

(খ) মাত্র $2$টা নির্দিষ্ট স্কুল এবং ঠিক ওই দুটোরই তুলনা লক্ষ্য হলে — fixed ধরাই যুক্তিসঙ্গত। random effect-এর পুরো ধারণা ("একটা জনগোষ্ঠী থেকে exchangeable নমুনা, যার variance estimate করি") মাত্র দুটো নির্দিষ্ট, পুনরায়-নমুনাযোগ্য-নয় level-এ অর্থহীন; দুটো level থেকে $\sigma_u^2$ নির্ভরযোগ্যভাবে estimate-ও করা যায় না।

(গ) parameter-সংখ্যা: - fixed ধরলে: প্রতিটা স্কুলের আলাদা intercept ⇒ $J=25$টা parameter (ছাড়াও $\beta_1$)। - random ধরলে: শুধু একটা variance parameter $\sigma_u^2$ — পৃথক $u_j$-গুলো estimate করা হয় না, বরং তাদের distribution ($\mathcal N(0,\sigma_u^2)$) estimate হয় (পরে BLUP হিসেবে $\hat u_j$ predict করা যায়, কিন্তু তারা স্বাধীন free parameter নয়)।

$25$ parameter → $1$ parameter: এই parsimony-ই random-effect ব্যবহারের প্রধান কারণ, বিশেষত $J$ বড় হলে।

৭.২ ★ — clustering-এ pooled OLS-এর SE কেন ভুল¶

(ক) positive within-cluster correlation মানে একই স্কুলের শিক্ষার্থীরা একটা সাধারণ $u_j$ ভাগ করে — তাই তাদের স্কোর পরস্পর-সম্পর্কিত, স্বাধীন নয়। একই স্কুলের দ্বিতীয় শিক্ষার্থী প্রথমজনের তুলনায় কম নতুন তথ্য যোগ করে (অনেকটা redundant)। ফলে $N=533$ পংক্তি থাকলেও কার্যকর (effective) স্বাধীন তথ্য-একক এর চেয়ে কম।

(খ) সবচেয়ে ক্ষতিগ্রস্ত between-cluster (স্কুল-স্তরের) পরিমাণ — যেমন intercept। তাদের জন্য কার্যকর নমুনা-আকার কাছাকাছি $J=25$ (স্বাধীন স্কুলের সংখ্যা), $N=533$ নয়; তাই OLS intercept SE মারাত্মক underestimated ($0.867$ vs সঠিক $1.416$)। hours একটা within-cluster predictor (প্রতিটা স্কুলের ভেতরেই বৈচিত্র্যময়), তাই অপেক্ষাকৃত কম প্রভাবিত — এমনকি mixed SE ($0.120$) OLS-এর ($0.149$)-এর চেয়ে একটু ছোটও হতে পারে (random-intercept noise সরিয়ে slope আরও সুনির্দিষ্ট হয়)।

(গ) underestimated SE ⇒ Wald $t$/$z=\hat\beta/\mathrm{SE}$ কৃত্রিমভাবে স্ফীত, $p$-value মিথ্যা-ছোট, confidence interval অতি-সংকীর্ণ ⇒ over-confident, anti-conservative inference (আসলে non-significant প্রভাবকে significant দেখানোর ঝুঁকি — বেশি false positive)।

৭.৩ ★ — complete / no / partial pooling¶

(ক) complete pooling। সব স্কুল মিলিয়ে একটাই regression (score ~ hours, school উপেক্ষা) — অর্থাৎ pooled OLS। হারায়: between-school বৈচিত্র্য পুরোপুরি ($u_j$ নেই), এবং clustering-জনিত SE-সংশোধন; ফলে §৭.২-এর over-confident inference।

(খ) no pooling। প্রতিটা স্কুলের সম্পূর্ণ আলাদা স্বাধীন regression ($25$টা পৃথক intercept, কোনো শেয়ার নেই)। সমস্যা: ছোট স্কুল ($n_j=10$)-কে শুধু তার নিজের $10$টা point দিয়ে fit করা হয় ⇒ estimate noisy, overfit; চরম (extreme) মান বেরিয়ে আসে যা শুধু নমুনা-শব্দ। তথ্য স্কুল-জুড়ে শেয়ার না-হওয়ায় ছোট গোষ্ঠীর দুর্বলতা সরাসরি estimate-এ পড়ে।

(গ) partial pooling। mixed model যা করে — প্রতিটা স্কুলের estimate-কে গ্র্যান্ড-গড়ের দিকে shrink করে, factor $\lambda_j=\frac{n_j\sigma_u^2}{n_j\sigma_u^2+\sigma_\varepsilon^2}$ অনুযায়ী। এটা দুই চরমের আপস: ছোট/noisy স্কুল বেশি টানা (no-pooling-এর overfit এড়ায়), বড়/তথ্যবহুল স্কুল কম টানা (complete-pooling-এর মতো সব এক করে ফেলে না)। ফলে no-pooling-এর চেয়ে কম variance, complete-pooling-এর চেয়ে কম bias — তথ্য "ধার" করে গোষ্ঠী-জুড়ে।

৭.৪ ★★ — ICC-এর ব্যাখ্যা¶

(ক) মোট variance $=\hat\sigma_u^2+\hat\sigma_\varepsilon^2=37.64+63.69=101.33$। স্কুল-ভেদ ভগ্নাংশ $$\rho=\frac{37.64}{101.33}=0.3714\approx0.371.$$ অর্থ: outcome-এর মোট পরিবর্তনশীলতার $\sim37\%$ আসে স্কুল-থেকে-স্কুল ভিন্নতা থেকে, বাকি $\sim63\%$ একই স্কুলের শিক্ষার্থী-থেকে-শিক্ষার্থী ভিন্নতা থেকে।

(খ) একই স্কুলের দুই শিক্ষার্থী একই random intercept $u_j$ ভাগ করে। তাদের covariance $=\operatorname{Cov}(u_j+\varepsilon_{ij},\,u_j+\varepsilon_{i'j})=\sigma_u^2$ (শুধু ভাগ-করা অংশ), আর প্রতিটির variance $=\sigma_u^2+\sigma_\varepsilon^2$। তাই correlation $$\frac{\sigma_u^2}{\sigma_u^2+\sigma_\varepsilon^2}=\rho.$$ "variance-ভগ্নাংশ" আর "জোড়া-correlation" — দুটোই গাণিতিকভাবে একই অনুপাত, তাই একই সংখ্যা (পূর্ণ প্রমাণ §৭.১১)।

(গ) $\rho=0$ ⇒ $\sigma_u^2=0$, স্কুল-ভেদ নেই, observation কার্যত স্বাধীন ⇒ mixed model pooled OLS-এ ফিরে যায়। $\rho=1$ ⇒ $\sigma_\varepsilon^2=0$, within-school কোনো noise নেই — একই স্কুলের সবাই একই মান, স্কুলই সম্পূর্ণ নির্ধারক।

৭.৫ ★★ — random intercept বনাম random slope¶

(ক) random-intercept model-এ প্রতিটা স্কুলের regression-রেখা $y=(\beta_0+u_j)+\beta_1 x$ — অর্থাৎ সবার slope (ঢাল) অভিন্ন ($\beta_1$), শুধু উচ্চতা (intercept $\beta_0+u_j$) আলাদা। জ্যামিতিকভাবে: সব রেখা পরস্পর সমান্তরাল, কেবল উলম্বভাবে স্থানান্তরিত।

(খ) hours-এর প্রভাবও স্কুল-ভেদে বদলালে random slope যোগ করতে হবে: $$y_{ij}=\beta_0+\beta_1 x_{ij}+u_{0j}+u_{1j}x_{ij}+\varepsilon_{ij},$$ যেখানে স্কুল $j$-এর কার্যকর ঢাল $=\beta_1+u_{1j}$ এবং কার্যকর intercept $=\beta_0+u_{0j}$। এখন রেখাগুলো আর সমান্তরাল নয় — উচ্চতা ও ঢাল দুটোই স্কুল-ভেদে ভিন্ন।

(গ) নতুন parameter: intercept-random ($u_{0j}$) ও slope-random ($u_{1j}$)-এর মধ্যে covariance/correlation $\sigma_{01}$ (random-effects covariance matrix $\Sigma=\begin{pmatrix}\sigma_{u0}^2 & \sigma_{01}\\ \sigma_{01} & \sigma_{u1}^2\end{pmatrix}$-এর off-diagonal)। এটা গুরুত্বপূর্ণ কারণ এটা ধরে রাখে যেমন: উঁচু-intercept স্কুলে কি hours-এর রিটার্নও বেশি? — না ধরলে গোষ্ঠী-স্তরের dependence-কাঠামো ভুলভাবে মডেল হয়।

৭.৬ ★★ — REML বনাম ML¶

(ক) সাধারণ ML variance component estimate করার সময় কার্যত ধরে নেয় fixed effect $\hat\beta$ "জানা/নিখুঁত"। বাস্তবে $\beta$ estimate করতে কিছু degrees of freedom খরচ হয়; ML সেই খরচ হিসাবে না-নেওয়ায় residual-spread-কে কম দেখে — তাই $\hat\sigma^2$ biased-low (একটু ছোট)।

(খ) এটা ঠিক ৫.৩-এর degrees-of-freedom ধারণা: sample variance-এ $\frac{1}{n}\sum(x_i-\bar x)^2$ biased-low (কারণ $\bar x$ data থেকেই estimate), তাই $\frac{1}{n-1}$ দিয়ে ভাগ (Bessel সংশোধন)। REML কার্যত fixed-effect-এর df বাদ দিয়ে (residual-space-এ project করে) variance estimate করে — এটাই $n-1$-দিয়ে-ভাগের mixed-model অ্যানালগ; ML হলো $n$-দিয়ে-ভাগের অ্যানালগ।

(গ) ভিন্ন fixed-effect structure তুলনা (যেমন hours-সহ বনাম hours-ছাড়া) করতে ML লাগবে। কারণ REML-likelihood fixed-effect-কে integrate-out/project করে দেয় বলে তা fixed-design-নির্ভর — ভিন্ন fixed structure-এর দুটো REML-likelihood তুলনাযোগ্য নয়। বিপরীতে, fixed structure এক রেখে variance/random-structure তুলনা করতে REML ভালো (unbiased variance)। (চলমান fit REML — তাই canonical $\hat\sigma_u=6.14$, logLik $-1894.5$ সব REML-estimate।)

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

৭.৭ ★ — variance component থেকে ICC¶

(ক) $$\rho=\frac{\hat\sigma_u^2}{\hat\sigma_u^2+\hat\sigma_\varepsilon^2}=\frac{37.64}{37.64+63.69}=\frac{37.64}{101.33}=0.3714\approx0.371.\ \checkmark$$

(খ) $$\hat\sigma_u=\sqrt{37.64}=6.135\approx6.14,\qquad \hat\sigma_\varepsilon=\sqrt{63.69}=7.981\approx7.98.$$ true $\sigma_u=6$, $\sigma_\varepsilon=8$-এর কাছাকাছি — পার্থক্য শুধু সসীম-নমুনা (finite-sample) আন্দাজ-ভুল, যা $J=25$, $N=533$-এ স্বাভাবিক।

(গ) মোট outcome-variance $=37.64+63.69=101.33$; মোট SD $=\sqrt{101.33}\approx10.07$।

৭.৮ ★★ — shrinkage factor $\lambda_j$: ছোট বনাম বড় group¶

$\hat\sigma_u^2=37.64$, $\hat\sigma_\varepsilon^2=63.69$।

(ক) ছোট স্কুল $n_j=10$: $$\lambda_{10}=\frac{10\cdot37.64}{10\cdot37.64+63.69}=\frac{376.4}{376.4+63.69}=\frac{376.4}{440.09}=0.855.$$

(খ) বড় স্কুল $n_j=30$: $$\lambda_{30}=\frac{30\cdot37.64}{30\cdot37.64+63.69}=\frac{1129.2}{1129.2+63.69}=\frac{1129.2}{1192.89}=0.947.$$

(গ) ছোট স্কুলের $\lambda$ ছোট ($0.855<0.947$), তাই shrink বেশি ($1-\lambda=0.145$ বনাম মাত্র $0.053$)। কারণ কম data ($n_j=10$) মানে স্কুলের নিজস্ব গড় কম নির্ভরযোগ্য (বেশি sampling noise), তাই model সেই estimate-কে আরও বেশি population-গড়ের দিকে টানে — noisy চরম মানকে নিয়ন্ত্রণ করে। বড় স্কুলের গড় বেশি নির্ভরযোগ্য, তাই কম টানা হয়।

৭.৯ ★★ — BLUP ও shrink-করা স্কুল-গড়¶

(ক) স্কুল $0$, $n_0=20$: $$\lambda_0=\frac{20\cdot37.64}{20\cdot37.64+63.69}=\frac{752.8}{752.8+63.69}=\frac{752.8}{816.49}=0.922.$$

(খ) raw (no-pooling) intercept-বিচ্যুতি $\approx+6.56$ ধরে, partial-pooling BLUP $$\hat u_0\approx\lambda_0\times6.56=0.922\times6.56=6.05.\ \checkmark$$ এটাই shrinkage-এর সারমর্ম: BLUP হলো স্কুলের raw-বিচ্যুতিকে $\lambda_0$ দিয়ে গ্র্যান্ড-গড়ের দিকে টানা (এখানে মাত্র $\sim8\%$ টানা, কারণ $n_0=20$ বেশ তথ্যবহুল)। (নোট: §ঘ-script-এ unconditional raw group-mean deviation $\approx6.58$; hours-adjusted raw deviation এর কাছাকাছি, এবং $\lambda_0$ দিয়ে গুণে BLUP $6.05$ মেলে।)

(গ) স্কুল $0$-এর model-implied intercept (hours$=0$-এ প্রত্যাশিত গড় স্কোর): $$\hat\beta_0+\hat u_0=50.82+6.05=56.87.$$

৭.১০ ★★ — design effect¶

$\bar n=533/25=21.32$, $\rho=0.371$।

(ক) $$\text{Deff}=1+(\bar n-1)\rho=1+(21.32-1)\times0.371=1+20.32\times0.371=1+7.54=8.54.$$

(খ) between-cluster SE-গুণক: $$\sqrt{\text{Deff}}=\sqrt{8.54}\approx2.92.$$

(গ) intercept-এ প্রয়োগ: $0.867\times2.92\approx2.53$। এটা balanced-আনুমানিক একটা রক্ষণশীল অনুমান; data unbalanced ও mixed model partial-pooling করে বলে প্রকৃত mixed SE ($1.416$) এর চেয়ে ছোট, কিন্তু দিক একই ও দ্ব্যর্থহীন: naive OLS intercept SE ($0.867$) clustering উপেক্ষা করে মারাত্মক ছোট, আর সঠিক (mixed) SE অনেক বড়। এটাই §৭.২-এর গুণগত যুক্তির পরিমাণগত নিশ্চয়তা — design effect-ই বলে দেয় "কত গুণ ভুল"।

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

৭.১১ ★★ — compound symmetry: within-group covariance $=\sigma_u^2$¶

model: $y_{ij}=\underbrace{\beta_0+\beta_1 x_{ij}}_{\text{fixed (ধ্রুবক, দেওয়া } x)}+u_j+\varepsilon_{ij}$, যেখানে $u_j\sim\mathcal N(0,\sigma_u^2)$, $\varepsilon_{ij}\sim\mathcal N(0,\sigma_\varepsilon^2)$, এবং সব $u_j,\varepsilon_{ij}$ পরস্পর স্বাধীন।

(ক) fixed অংশ ধ্রুবক, তাই covariance-এ অবদান $0$ (ধ্রুবকের covariance শূন্য)। একই স্কুল $j$, দুই শিক্ষার্থী $i\ne i'$: $$ \operatorname{Cov}(y_{ij},y_{i'j})=\operatorname{Cov}(u_j+\varepsilon_{ij},\,u_j+\varepsilon_{i'j}). $$ bilinearity দিয়ে চারটা পদে বিস্তার (২.৬-এর সরাসরি প্রয়োগ): $$ =\underbrace{\operatorname{Cov}(u_j,u_j)}{=\sigma_u^2} +\underbrace{\operatorname{Cov}(u_j,\varepsilon}){=0} +\underbrace{\operatorname{Cov}(\varepsilon},u_j){=0} +\underbrace{\operatorname{Cov}(\varepsilon =\sigma_u^2. $$ শেষ তিনটা পদ শূন্য — কারণ },\varepsilon_{i'j})}_{=0\ (i\ne i')$u_j$ ও যেকোনো $\varepsilon$ স্বাধীন, এবং $i\ne i'$ হওয়ায় $\varepsilon_{ij},\varepsilon_{i'j}$ ভিন্ন ও স্বাধীন। তাই $\boxed{\operatorname{Cov}(y_{ij},y_{i'j})=\sigma_u^2}$ ∎

(খ) variance: $$ \operatorname{Var}(y_{ij})=\operatorname{Var}(u_j+\varepsilon_{ij})=\operatorname{Var}(u_j)+\operatorname{Var}(\varepsilon_{ij})=\sigma_u^2+\sigma_\varepsilon^2 $$ (cross-term $2\operatorname{Cov}(u_j,\varepsilon_{ij})=0$ স্বাধীনতায়)। তাই within-group correlation: $$ \operatorname{Corr}(y_{ij},y_{i'j})=\frac{\operatorname{Cov}(y_{ij},y_{i'j})}{\sqrt{\operatorname{Var}(y_{ij})\operatorname{Var}(y_{i'j})}}=\frac{\sigma_u^2}{\sigma_u^2+\sigma_\varepsilon^2}=\rho. $$ এটাই ICC-র দ্বিতীয় ব্যাখ্যা — পরিমাণগতভাবে variance-ভগ্নাংশের অভিন্ন।

(গ) ভিন্ন স্কুল $j\ne j'$: $u_j,u_{j'}$ স্বাধীন, এবং কোনো $\varepsilon$ শেয়ার নেই, তাই সব covariance-পদ শূন্য ⇒ $\operatorname{Cov}(y_{ij},y_{i'j'})=0$। ফলে গোষ্ঠী $j$-এর covariance matrix (আকার $n_j\times n_j$): $$ \Sigma_j=\sigma_\varepsilon^2 I_{n_j}+\sigma_u^2\,\mathbf 1\mathbf 1^\top =\begin{pmatrix}\sigma_u^2+\sigma_\varepsilon^2 & \sigma_u^2 & \cdots & \sigma_u^2\ \sigma_u^2 & \sigma_u^2+\sigma_\varepsilon^2 & \cdots & \sigma_u^2\ \vdots & & \ddots & \vdots\ \sigma_u^2 & \sigma_u^2 & \cdots & \sigma_u^2+\sigma_\varepsilon^2\end{pmatrix}. $$ এই গঠনকে compound symmetry বলে: প্রতিটা diagonal সমান ($\sigma_u^2+\sigma_\varepsilon^2$), প্রতিটা off-diagonal সমান ($\sigma_u^2$) — অর্থাৎ গোষ্ঠীর সব জোড়ার covariance অভিন্ন, যেকোনো অদলবদলে (exchangeable) matrix একই থাকে।

৭.১২ ★★★ — shrinkage factor $\lambda_j$-এর উৎপত্তি¶

সরল model $y_{ij}=\mu+u_j+\varepsilon_{ij}$ (predictor নেই), স্কুল-গড় $\bar y_j=\frac1{n_j}\sum_i y_{ij}$, BLUP $\hat u_j=\lambda_j(\bar y_j-\mu)$।

(ক) গোষ্ঠী-গড় বিস্তার: $$ \bar y_j=\mu+u_j+\bar\varepsilon_j,\qquad \bar\varepsilon_j=\frac1{n_j}\sum_{i=1}^{n_j}\varepsilon_{ij}. $$ $\varepsilon$-গুলো i.i.d. $\mathcal N(0,\sigma_\varepsilon^2)$, তাই $$ \operatorname{Var}(\bar\varepsilon_j)=\frac{\sigma_\varepsilon^2}{n_j}. $$ সুতরাং $\bar y_j-\mu$-কে দেখা যায় signal $u_j$ (prior variance $\sigma_u^2$) আর noise $\bar\varepsilon_j$ (variance $\sigma_\varepsilon^2/n_j$)-এর যোগফল হিসেবে।

(খ) Normal–Normal conjugacy। prior $u_j\sim\mathcal N(0,\sigma_u^2)$ এবং likelihood $\bar y_j-\mu\mid u_j\sim\mathcal N(u_j,\ \sigma_\varepsilon^2/n_j)$। দুই Normal-এর posterior mean হলো inverse-variance-weighted: $$ \mathbb E[u_j\mid \bar y_j]=\frac{1/(\sigma_\varepsilon^2/n_j)}{1/\sigma_u^2+1/(\sigma_\varepsilon^2/n_j)}(\bar y_j-\mu) =\frac{\sigma_u^2}{\sigma_u^2+\sigma_\varepsilon^2/n_j}(\bar y_j-\mu). $$ লব ও হরকে $n_j$ দিয়ে গুণ: $$ =\frac{n_j\sigma_u^2}{n_j\sigma_u^2+\sigma_\varepsilon^2}(\bar y_j-\mu)=\lambda_j(\bar y_j-\mu),\qquad \boxed{\lambda_j=\frac{n_j\sigma_u^2}{n_j\sigma_u^2+\sigma_\varepsilon^2}}\ \ \blacksquare $$ স্বজ্ঞা: কম-noise উৎস (বেশি $n_j$ ⇒ ছোট $\sigma_\varepsilon^2/n_j$) বেশি weight পায়, তাই estimate স্কুলের নিজস্ব গড়ের দিকে বেশি ঝোঁকে; বেশি-noise হলে prior (গ্র্যান্ড-গড় $0$)-এর দিকে বেশি টানে।

(গ) সীমা: - $n_j\to\infty\Rightarrow\sigma_\varepsilon^2/n_j\to0\Rightarrow\lambda_j\to1$ — no shrink: বড় স্কুলের নিজস্ব গড়কেই বিশ্বাস (no-pooling-এর কাছাকাছি)। - $n_j\to0$ (বা $\sigma_u^2\to0$) $\Rightarrow\lambda_j\to0$ — পূর্ণ shrink: তথ্য নেই, BLUP $\to0$, অর্থাৎ গ্র্যান্ড-গড়ই estimate (complete-pooling)।

এটাই §৭.৮-এর ছোট-বনাম-বড় আচরণের গাণিতিক ভিত্তি: $\lambda_j$ monotone (একঘেয়ে)-ভাবে $n_j$-তে বাড়ে, তাই ছোট স্কুল কম $\lambda$ ⇒ বেশি টানা। partial pooling এভাবে স্বয়ংক্রিয়ভাবে প্রতিটা গোষ্ঠীর data-পরিমাণ অনুযায়ী বিশ্বাস বণ্টন করে — no-pooling ও complete-pooling-এর মধ্যে data-চালিত আপস।

ঘ · কোডিং (Python)¶

৭.১৩ ★★★ — পূর্ণ mixed-model pipeline¶

নিচের script এ-ই data-generation-ক্রমে (প্রথমে সব $n_j$, তারপর সব $u_j$, তারপর প্রতিটা গোষ্ঠীর hours→eps) উপরের canonical সংখ্যাগুলো হুবহু পুনরুৎপাদন করে। (random-call-এর ক্রম বদলালে seed-আউটপুট বদলায়, তাই এই ক্রম রাখা জরুরি।)

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

# ---- 1. seeded data generation (canonical order) ----
rng = np.random.default_rng(20260619)
J = 25
n_j = rng.integers(10, 31, size=J)         # group sizes 10..30  -> N = 533
u   = rng.normal(0, 6, size=J)             # school random intercepts, sigma_u = 6
rows = []
for j in range(J):
    hours = rng.uniform(0, 10, size=n_j[j])
    eps   = rng.normal(0, 8, size=n_j[j])  # residual, sigma_eps = 8
    score = 50 + 2*hours + u[j] + eps      # true: b0=50, b1=2
    for k in range(n_j[j]):
        rows.append({"school": j, "hours": hours[k], "score": score[k]})
df = pd.DataFrame(rows)
print(f"J={J}  N={len(df)}  group sizes {n_j.min()}-{n_j.max()}  n_0={n_j[0]}")
# -> J=25  N=533  group sizes 10-30  n_0=20

# ---- 2. pooled OLS (complete pooling) — SEs are WRONG under clustering ----
ols = smf.ols("score ~ hours", data=df).fit()
print(f"OLS : b1={ols.params['hours']:.3f} (SE {ols.bse['hours']:.3f}), "
      f"intercept SE={ols.bse['Intercept']:.3f}")
# -> OLS : b1=1.800 (SE 0.149), intercept SE=0.867

# ---- 3. mixed model (random intercept; partial pooling), REML ----
mix = smf.mixedlm("score ~ hours", data=df, groups=df["school"]).fit(reml=True)
b0, b1 = mix.params["Intercept"], mix.params["hours"]
se0, se1 = mix.bse["Intercept"], mix.bse["hours"]
print(f"MIX : b0={b0:.2f} (SE {se0:.3f}), b1={b1:.3f} (SE {se1:.3f})")
# -> MIX : b0=50.82 (SE 1.416), b1=1.878 (SE 0.120)
print(f"intercept SE: OLS {ols.bse['Intercept']:.3f}  vs  mixed {se0:.3f}  "
      f"(clustering inflates it)")
# -> intercept SE: OLS 0.867 vs mixed 1.416

# ---- 4. variance components ----
var_u = float(mix.cov_re.iloc[0, 0])   # between-school variance  (sigma_u^2)
var_e = float(mix.scale)               # residual variance        (sigma_eps^2)
print(f"var_u={var_u:.2f} (sigma_u={var_u**0.5:.2f}), "
      f"var_e={var_e:.2f} (sigma_eps={var_e**0.5:.2f})")
# -> var_u=37.64 (sigma_u=6.14), var_e=63.69 (sigma_eps=7.98)

# ---- 5. ICC ----
icc = var_u / (var_u + var_e)
print(f"ICC = {icc:.3f}")          # -> ICC = 0.371
print(f"logLik = {mix.llf:.1f}")   # -> logLik = -1894.5

# ---- 6. BLUPs (random effects) ----
re = mix.random_effects                       # dict: school -> Series(['Group'])
blup0 = float(re[0].iloc[0])
print(f"BLUP school 0 (n_0=20) = {blup0:+.2f}")   # -> +6.05

# shrinkage check for school 0
lam0 = n_j[0]*var_u / (n_j[0]*var_u + var_e)
grand = df["score"].mean()
raw_dev0 = df.loc[df.school == 0, "score"].mean() - grand   # ~ +6.58 (uncond.)
print(f"lambda_0={lam0:.3f}, raw dev~{raw_dev0:+.2f}, "
      f"lambda_0*raw~{lam0*raw_dev0:+.2f}  (vs BLUP {blup0:+.2f})")
# -> lambda_0=0.922 ; BLUP = shrink of raw group deviation toward grand mean

# ---- 7. shrinkage is monotone in n_j: small groups pulled MORE ----
for nj in (10, 20, 30):
    lam = nj*var_u/(nj*var_u+var_e)
    print(f"  n_j={nj:2d} -> lambda={lam:.3f}  (shrink {1-lam:.3f})")
# -> n=10 lambda=0.855 ; n=20 lambda=0.922 ; n=30 lambda=0.947

প্রত্যাশিত আউটপুট (canonical):

পরিমাণ	মান
pooled OLS slope $\hat\beta_1$ (SE)	$1.800$ ($0.149$)
OLS intercept SE	$0.867$
mixed $\hat\beta_0$ (SE)	$50.82$ ($1.416$)
mixed $\hat\beta_1$ (SE)	$1.878$ ($0.120$)
$\hat\sigma_u^2$ ($\hat\sigma_u$)	$37.64$ ($6.14$)
$\hat\sigma_\varepsilon^2$ ($\hat\sigma_\varepsilon$)	$63.69$ ($7.98$)
ICC $\rho$	$0.371$
logLik (REML)	$-1894.5$
BLUP স্কুল $0$ ($n_0=20$)	$+6.05$
$\lambda_0$	$0.922$

মূল পর্যবেক্ষণ। 1. SE-সংশোধন: intercept SE OLS $0.867\to$ mixed $1.416$ — clustering ধরায় প্রায় $1.6$ গুণ বড় (§৭.১০-এর design effect $\sqrt{8.5}\approx2.9$-এর সাথে দিক-সঙ্গতিপূর্ণ)। OLS clustering উপেক্ষা করে over-confident। 2. shrinkage: BLUP-গুলো প্রতিটা স্কুলের raw-deviation থেকে গ্র্যান্ড-গড়ের দিকে $\lambda_j$-গুণ টানা; (group_mean − grand_mean) বনাম BLUP scatter করলে সব point origin-মুখী একটা $<1$-ঢালের রেখায় বসে — ছোট স্কুল বেশি টানা। 3. variance components: REML-estimate ($37.64$, $63.69$) true ($36$, $64$)-এর কাছাকাছি; ICC $0.371$ মানে মোট পরিবর্তনশীলতার $\sim37\%$ স্কুল-ভেদ থেকে।

(ঐচ্ছিক যাচাই: mix.summary() দিয়ে পূর্ণ table; mix.cov_re, mix.scale, mix.llf সরাসরি variance/logLik দেয়। random-slope সম্প্রসারণ: `smf.mixedlm("score ~ hours", df, groups=df["sc

পরিমাণ	মান
pooled OLS slope \(\hat\beta_1\) (SE)	\(1.800\) (\(0.149\))
OLS intercept SE	\(0.867\)
mixed \(\hat\beta_0\) (SE)	\(50.82\) (\(1.416\))
mixed \(\hat\beta_1\) (SE)	\(1.878\) (\(0.120\))
\(\hat\sigma_u^2\) (\(\hat\sigma_u\))	\(37.64\) (\(6.14\))
\(\hat\sigma_\varepsilon^2\) (\(\hat\sigma_\varepsilon\))	\(63.69\) (\(7.98\))
ICC \(\rho\)	\(0.371\)
logLik (REML)	\(-1894.5\)
BLUP স্কুল \(0\) (\(n_0=20\))	\(+6.05\)
\(\lambda_0\)	\(0.922\)