\chapter{تحلیل بیزی تقریبی داده‌های بقا}
\section{مقدمه}
در این فصل ابتدا توضیحاتی در خصوص رهیافت بیزی تقریبی معرفی شده توسط رو و مارتینو
\cite{RUE}
ارائه می‌گردد و مدل کاکس فضایی با تابع خطر پایه نمایی به‌صورت یک مدل گاوسی پنهان بیان می‌شود و سپس با استفاده از رهیافت بیزی تقریبی به تحلیل بیزی این مدل پرداخته می‌شود.
\section{مدل‌های گاوسی پنهان}
در مدل‌های جمعی ساختاری فرض می‌کنیم متغیر پاسخ $ y_{i} $ متعلق به خانواده نمایی با میانگین $ \mu_{i} $ باشد که به پیش‌گوی جمعی $ \eta_{i} $ از طریق تابع پیوند $ g(\cdot) $ مرتبط است. که این ارتباط به‌صورت $ g(\mu_{i})=\eta_{i} $ قابل تعریف است.
پیش‌گوی جمعی یعنی، $ \eta_{i}  $ اثرات متغیر‌های تبیینی مختلف را با فرم
\begin{eqnarray}
\eta_{i}=\alpha+\sum_{j=1}^{n_{f}}f^{(j)}(u_{ji})+\sum_{k=1}^{n_{\beta}}\beta_{k}z_{ki}+\epsilon_{i}
\end{eqnarray}
بیان می‌کند. که مجموعه $\lbrace f^{(j)}(\cdot)\rbrace $ توابعی نامعلوم از متغیرهای تبیینی $ \bm u_{i}=(u_{1i},\cdots,u_{n_{f_{i}}}) $، می‌باشد که برای بیان ارتباط خطی میان متغیر‌های تبیینی استفاده می‌شود و می‌تواند شامل اثرات غیرخطی مانند روند‌های زمانی، اثرات فصلی، شیب‌ها و عرض از مبدأ‌های تصادفی در مدل باشند و در مدل‌هایی که وابستگی زمانی و یا فضایی دارند کاربرد دارد.
\begin{figure}\label{fig200} 
\begin{center}
\includegraphics[scale=.6]{fig200.png}
\end{center}
\tiny{\caption{خواص مارکوفی مختلف رانشان می‌دهد. $(a)$ خاصیت مارکوفی دوبه‌دو: دو رأس سیاه رنگ به شرط رأس خاکستری داده شده مستقل شرطی هستند. $(b)$ خاصیت مارکوفی موضعی: رأُس‌های سیاه و سفید به شرط رأس خاکستری داده شده مستقل شرطی هستند. $(c)$ خاصیت مارکوفی فراموضعی: رأس‌های سیاه رنگ و هاشور خورده به شرط رأس‌های خاکستری داده شده مستقل شرطی هستند. \label{fig200}}}.
\end{figure}
\begin{theorem}\label{th3}
(رو و هلد\cite{Rue})،
اگر $\bm x$ یک $GMRF$ نسبت به گراف $\Im =(\nu ,\varepsilon)$ باشد. روابط زیر معادل هم هستند،\\
در خاصیت مارکوفی دوبه‌دو: اگر $\{i,j\}\notin\varepsilon$ و $i\neq j$ باشد آن‌گاه $x_i\perp x_j\mid \bm x_{-ij}$.\\
در خاصیت مارکوفی موضعی برای هر $i\in \nu$ خواهیم داشت: $x_i\perp \bm x_{-\{i,ne(i)\}}\mid \bm x_{ne(i)}$.\\
و خاصیت مارکوفی فراموضعی برای تمام مجموعه‌های $A$ و $B$ و $C$ مجزا به‌طوری که $C$ مجزا از $A$ و $B$ باشد و $A$ و $B$ ناتهی هستند، به‌صورت $\bm x_A\perp \bm x_B\mid \bm x_C$ می‌باشد.
\end{theorem}
شکل \ref{fig200}، قضیۀ \ref{th3} را شرح می‌دهد. اثبات به‌صورت نتیجه‌ای از یک نتیجۀ کلی‌تر می‌باشد. 
 بنابراین میدان تصادفی پنهان یک میدان تصادفی مارکوفی گاوسی $(GMRF)$ با یک ماتریس دقت تُنُک $\bm Q (\bm \theta)$ است. به‌همین دلیل می‌توان از روش‌های عددی مربوط به ماتریس‌های تُنُک که بسیار سریع‌تر از روش‌های محاسباتی ماتریس‌های چگال هستند استفاده کرد. ویژگی دوم که باید این مدل‌ها داشته باشند این است که تعداد ابر پارامتر‌ها باید کم و معمولاً $ m\leq 6 $ باشد که در این‌صورت شرایط لازم برای استنباط‌های سریع فراهم خواهد بود (ایدسویک و همکاران\cite{Eidsvik}).
\section{میدان تصادفی مارکوفی گاوسی}
میدان تصادفی مارکوفی گاوسی براساس گراف قابل تعریف است. گراف $\Im$ مجموعه‌ای از رأس‌ها است که توسط خانواده‌ای از یال‌ها به‌هم وصل شده‌اند و به‌صورت زوج مرتب $(\nu,\varepsilon)$ نشان داده می‌شود، که در آن $\nu$ مجموعه‌ایی متناهی و غیر تهی از رئوس و $\varepsilon$ یال‌های آن است. اگر رأس‌ها به‌صورت $\nu=\lbrace 1,\cdots,n \rbrace$ باشند گراف را نشان‌دار\LTRfootnote{Labeled Graph} گویند. فرض کنید بردار تصادفی $\bm x=(x_{1},\cdots,x_{n})^{\prime}$ دارای توزیع نرمال چند متغیره با میانگین $\bm \mu $ و ماتریس کواریانس $ \bm\Sigma $ باشد. به‌علاوه $\Im=(\nu,\varepsilon)$ گرافی نشان‌دار باشد، که در آن $\varepsilon$ شامل همه زوج‌های $\lbrace i,j \rbrace$ است به‌طوری که رأس‌های $i$ و $j$ هیچ یال مشترکی نداشته باشند، اگر و تنها اگر $ x_{i} \perp x_{j}\mid \bm x_{-ij}$ در این‌صورت $\bm x$ یک میدان تصادفی مارکوفی گاوسی نسبت به گراف $\Im$ نامیده می‌شود.
\begin{theorem}\label{th4}
(رو و هلد\cite{Rue}): اگر $\bm x$ دارای توزیع نرمال چند متغیره با میانگین $\bm \mu$ و ماتریس دقت $\bm Q>0$ باشد آن‌گاه درایه $ij$ام ماتریس $\bm Q$ یعنی $Q_{ij}$ برابر صفر است، اگر و تنها اگر $ x_{i} \perp x_{j}\mid \bm x_{-ij} $.
\end{theorem}
با توجه به قضیه \ref{th4} درایه‌های غیرصفر ماتریس $\bm Q$ گراف $\Im$ را مشخص می‌کنند و براساس آن‌ها می‌توان استقلال شرطی $x_{i}$  و $x_{j}$ را بررسی نمود، و به‌همین دلیل در این‌جا از ماتریس دقت به‌جای ماتریس کواریانس استفاده می‌شود. بنابراین، بردار تصادفی $\bm x=(x_{1},\cdots,x_{n})^{\prime} \in R^{n}$ یک میدان تصادفی مارکوفی گاوسی تحت گراف نشان‌دار $ G=(\nu,\varepsilon)$ با میانگین $\bm \mu$ و ماتریس دقت $\bm Q_{n\times n}>0$ است اگر و تنها اگر تابع چگالی آن به‌صورت
\begin{align*}
\pi(\bm x)=(2\pi)^{-\frac{n}{2}} \mid\bm Q\mid ^{\frac{1}{2}}\exp\lbrace -\frac{1}{2}(\bm x-\bm \mu )^{\tau} \bm Q(\bm x-\bm \mu)\rbrace,
\end{align*}
باشد، به‌گونه‌ای که درایه $(ij)$ام ماتریس دقت، یعنی $\bm Q_{ij}$، مخالف صفر است اگر و تنها اگر $\lbrace i,j \rbrace \in \varepsilon$. در بیشتر موارد $\bm Q$ ماتریسی تُنُک\LTRfootnote{Sparse Matrix} است و تنها $O(n)$ درایه از $n^{2}$ درایه ماتریس $\bm Q$ غیرصفر است (رو و مارتینو\cite{Ruee})، که این به‌علت ویژگی مارکوفی میدان تصادفی است. تُنُک بودن این ماتریس باعث افزایش سرعت در محاسبات و الگوریتم‌های تکراری مانند نمونه‌گیری از میدان تصادفی مارکوفی می‌شود.
\section{تقریب لاپلاس آشیانه‌ای جمع بسته}\label{rt1}
برای تحلیل بیزی تقریبی مدل‌های گاوسی پنهان، ابتدا توزیع‌های پسین کناری متغیرهای پنهان و ابر پارامترها به‌صورت 
\begin{eqnarray}
\label{post1}
\pi(\bm x_i\mid\bm{y} )&=&\int\pi(\bm x_i\mid \bf\theta,\bm y)\pi(\bf\theta\mid\bm y)d{\bm\theta}\\
\label{post2}\pi(\bm\theta_j\mid\bf y)&=&\int\pi(\theta\mid\bm y)d{\bm\theta}_{-j}
%\end{array}
\end{eqnarray}
تشکیل می‌شوند، که در آن $\bm x$ میدان پنهان و $\bm\theta$ بردار ابر پارامترها می‌باشد و $\bm\theta_{-j}$ بردار حاصل از حذف درایه $j$ام است و برای چگالی‌های پسینی حاشیه‌ای (\ref{post1}) ، (\ref{post2}) تقریب‌هایی به‌صورت $\tilde\pi(x\mid\bf y)$ و $\tilde\pi(\bf\theta_j\mid\bf y)$ ارائه می‌شود.\\
در روش \lr{INLA} ابتدا با استفاده از بسط تیلور و خطی‌سازی قسمت درستنمایی حول یک مقدار ثابت، توزیع $\pi(\bm x|\bm y,\bf\theta)$ با استفاده از یک توزیع نرمال چند متغیره $\tilde\pi_G(\bm x|\bm y,\bf\theta)$ تقریب زده می‌شود و مقدار چگالی در مد توزیع محاسبه می‌شود. سپس برای به‌دست آوردن پسین تقریبی ابر پارامترها با به‌کار بردن تقریب لاپلاس، تقریب پسین به‌صورت  
\begin{eqnarray*}
\tilde\pi(\bm\theta|\bf y)\propto\frac{\pi(\bm x,\bm\theta,\bf y)}
{\tilde\pi_G(\bm x|\bf y,\bm\theta)}|_{x=x^*(\bf\theta)}
\end{eqnarray*}
می‌شود،که در آن $x^{*}(\bf\theta)$ مد توزیع شرطی کامل $\bm x$ به‌شرط \lr{$\bm\theta$ } است که با یک الگوریتم  بهینه‌سازی مثل الگوریتم  نیوتن رافسون محاسبه می‌شود. برای تولید نقاط از توزیع پسین تقریبی ابر پارامترها ابتدا مد تابع $\tilde \pi (\bm\theta|\bf y)$ با روش‌های بهینه‌سازی به‌صورت $\bm\theta^*$ تعیین می‌شود. ماتریس هسین در مد محاسبه می‌شود و  عکس آن به‌صورت $H^{-1} = V \Lambda  V'$ تجزیه می‌شود، که در آن $V$ ماتریس بردارهای ویژه و$\Lambda $ ماتریس قطری مقادیر ویژه آن است. مبدأ مختصات را به مد $\bm\theta^{*}$ انتقال داده، فرمول مختصات در مبدأ $\bm\theta^{*}$ به‌صورت $\theta(z)=\bm\theta^* + V\Lambda^{-1}\bm z$ تعریف می‌شود که در آن $\bm z$ مقادیر استانداردشده هستند. با شروع از مبدأ مختصات جدید روی 
هر یک از محورها نقاطی به فاصله مقادیر صحیح مناسب $\delta_z$ به‌گونه‌ای اختیار می‌شوند، که شرط $\log \tilde \pi (\theta(0)|\bf y) - \log \tilde \pi (\theta(z)|\bf y) < \delta_\pi$ برقرار باشد، سپس به‌طور مشابه نقاط درون صفحات نیز تعیین می‌شوند. با به‌کار بردن این الگوریتم \lr{ $\theta_i$}هایی از توزیع $\pi(\bm\theta|\bf y)$ تولید می‌شوند که می‌توان از آن‌ها برای تقریب
توزیع $\pi (\bm\theta|\bf y)$ استفاده کرد. پس از تولید مقادیر $\bm\theta$ از روش فوق، اکنون تقریب لاپلاس 
توزیع شرطی کامل $\pi(x_i|\bm y,\bm\theta)$  به‌صورت
\begin{eqnarray}
\label{approx}
\tilde\pi_{LA}(x_i|\bm\theta,\bm y)\propto\frac{\pi(\bm x,\bm\theta,\bm y)}
{\tilde\pi_G(\bm x_{-i}|x_i,\bm\theta,\bm y)}|_{x_{-i}={x^*}_{-i}(x_i,\bf\theta)}
\end{eqnarray}
می‌باشد، که در آن $\bm x_{-i}$ بردار حاصل از حذف درایه $i$ام است.  
 $\tilde\pi_G(\bm x_{-i}|x_i,\bm\theta,\bf y)$ تقریب گاوسی از $\pi(\bm x_{-i}|x_i,\bm\theta,\bm y)$ که $x_i$ یک مقدار ثابت و ${\bm x^*}_{-i}(x_i,\bm\theta)$ مد $\pi(\bm x_{-i}|x_i,\bm\theta,\bm y)$ است. در نهایت پس از محاسبه توزیع تقریبی پسین ابرپارامترها و شرطی کامل متغیرهای پنهان اکنون می‌توان توزیع پسین کناری متغیرهای پنهان 
را از رابطه                    
\begin{equation}
\label{posterior2}\pi(x_{i}\mid\bf y)=
\int\pi(x_{i}\mid\bf y,\bm\theta)\pi(\bf\theta\mid\bf y)d_{\bf\theta}
\approx\sum_{k}\tilde{\pi}(x_{i}\mid\bm\theta_{k},\bf y)
\tilde{\pi}(\bm\theta_{k}\mid \bf y)\bigtriangleup_{k} 
\end{equation}
محاسبه کرد، که در آن جمع روی نقاط تولید شده برای $\bm\theta$ بسته شده است، با وزن‌های ناحیه‌ای $\bf\bigtriangleup_{k} $ برای  $K=1,\cdots,k$. چون تولید نقاط از روش توضیح داده شده، با فواصل یکسان در یک مشبکه منظم تولید می‌شوند، $\Delta_k$ها با هم برابر هستند. به‌دلیل مشابه بودن رابطه‌ی (\ref{approx}) به تقریب لاپلاس و تقریب توزیع پسین (\ref{posterior2}) با جمع روی مقادیر $\bm\theta$ تولید شده این روش، رهیافت تقریبی لاپلاس آشیانه‌ای جمع بسته نامگذاری شده است. در زیر بخش بعد برای روشن شدن تقریب لاپلاس مثالی آورده می‌شود. 
\subsection{تقریب لاپلاس}
تقریب لاپلاس هر تابع $g(X)$ را با توزیع نرمال تقریب می‌زند. ابتدا از سری تیلور $\log g(X)$ به‌صورت
\begin{eqnarray*}
\log g(X)\approx\log g(\hat{X})+\dfrac{\partial\log g(\hat{X})}{\partial X}(X-\hat{X})+\dfrac{1}{2}\dfrac{\partial^2\log g(\hat{X})}{\partial X^2}(X-\hat{X})^2
\end{eqnarray*}
چون $\dfrac{\partial\log g(\hat{X})}{\partial X}=0$ پس
\begin{eqnarray*}
\log g(X)\approx\log g(\hat{X})+\dfrac{1}{2}\dfrac{\partial^2\log g(\hat{X})}{\partial X^2}(x-\hat{X})^2
\end{eqnarray*}
با در نظر گرفتن $\hat{\sigma}^2=\dfrac{-1}{\dfrac{\partial^2\log g(\hat{X})}{\partial X^2}}$ می‌توان نوشت
\begin{eqnarray*}
\log g(X)\approx\log g(\hat{X})-\dfrac{1}{2\hat{\sigma}^2}(X-\hat{X})^2
\end{eqnarray*}
و به‌طور معادل $\int g(X)dx=\int e^{\log g(X)}dx\approx const\int\exp[-\dfrac{(X-\hat{X})^2}{2\hat{\sigma}^2}]dx$ و بنابراین
\begin{eqnarray*}
g(X)\approx Normal(\hat{X},\hat{\sigma}^2).
\end{eqnarray*}
\begin{example}
توزیع $\chi^2$ را در نظر بگیرید، $p(X)=\dfrac{g(X)}{c}=\dfrac{X^{\dfrac{k}{2}-1}e^{-\dfrac{X}{2}}}{c}$ بنابراین
\begin{eqnarray*}
&L(X)=\log g(X)=(\dfrac{k}{2}-1)\log x-\dfrac{X}{2}\cr
&L^{\prime}(X)=\dfrac{\partial\log g(X)}{\partial X}=(\dfrac{k}{2}-1)X^{-1}-\dfrac{1}{2}\cr
&L^{\prime\prime}(X)=\dfrac{\partial^2\log g(X)}{\partial X^2}=-(\dfrac{k}{2}-1)X^{-2}
\end{eqnarray*}
که با حل $L^{\prime}(X)=0$ مد به‌صورت $\hat{X}=k-2$ است و با جایگذاری مد در$-\dfrac{1}{L^{\prime\prime}}(X)$،    $\hat{\sigma}^2=2(K-2)$ به‌دست می‌آید. به‌طوری‌که، می‌توان $p(X)$ را به‌صورت
\begin{eqnarray*}
\tilde{p}(X)=Normal(k-2,2(k-2))
\end{eqnarray*}
تقریب زد. در شکل (\ref{fig ur}) به‌ازای مقادیر مختلف $k$ توزیع کای‌دو با تقریب لاپلاس گاوسی، تقریب زده شده است.
\begin{figure}
\begin{center}
\includegraphics[scale=.7]{fig100.png}
\end{center}
\caption{تقریب لاپلاس توزیع $\chi^2$ با نرمال.\label{fig ur}}
\end{figure}
\end{example}
از تقریب لاپلاس می‌توان هر توزیع پسین شرطی را به‌صورت
\begin{eqnarray*}
\tilde{p}(z\mid \omega)=\dfrac{p(x,z\mid\omega)}{\tilde{p}(x\mid z,\omega)},
\end{eqnarray*}
تقریب زد، که در آن $\tilde{p}(x\mid z,\omega)$ تقریب لاپلاس توزیع شرطی $x$ به شرط $z,\omega$ است.
\section{تحلیل بیزی تقریبی مدل کاکس با تابع خطر پایه نمایی تکه‌ای}
برای استفاده از تحلیل بیزی تقریبی در مدل‌های بقا، مشابه فصل دوم فرض کنید محور زمان به $K$ بازه زمانی $0<a_1<a_2<...<a_k<\infty$ افراز شده باشد که در آن $a_k$ از هر زمان بقا مشاهده شده در مطالعه بزرگ‌تر است. بنابراین $k$ بازه زمانی به‌صورت $I_k=(a_{k-1}, a_k], k=1, ...,K$ خواهیم داشت به‌طوری که $I_1=(0, t_1]$ و تابع پایه خطر در هر فاصله زمانی به‌صورت ثابت $h_0(t)=\lambda_k$ فرض می‌شود. هریک از بازه‌های تابع خطر پایه دارای توزیع نمایی با پارامتر $\lambda_k$ است. اکنون $(t_i, \delta_i)$ و $i=1, ..., m$، زمان بقای فرد $i$ام و $\delta_i$ نشان دهنده سانسورشدن باشد، آن‌گاه تابع خطر متناسب 
\begin{equation}
\label{eq1}
h_i(t)=h_0(t)\exp( \bm z'_i\bm\beta)
=\exp\{ \bm z'_i\bm\beta +b_k\}=\exp(\eta_{ik}),
\end{equation}
می‌باشد، که در آن $t\in(a_{k-1}, a_k)$ و 
 $b_k=\log(\lambda_k)$. در مدل  (\ref{eq1}) فرض می‌شود که اثرمتغیرهای کمکی روی لگاریتم تابع خطر، تابعی خطی از متغیرهای کمکی است. 
همچنین در مدل کاکس با تابع خطر نمایی تکه‌ای به‌صورت (\ref{eq1}) مستقیماً شکل غیر خطی متغیرها را وارد می‌کنیم. و مدل به‌صورت 
\begin{equation}
\label{eq2}
\eta_{ik}=\beta_0+\sum_{r=1}^{n_r} w_{ir} f_r (u_{ir})+\sum_{j=1}^{n_B}\beta_j z_{ji}+b_k,
\end{equation}
در نظر گرفته می‌شود، که در آن $n_r$ تعدادمتغیرهای کمکی با اثرات غیرخطی، $n_B$ تعداد متغیرهای کمکی با اثرات خطی،   $\beta_j$ ها ضرایب رگرسیونی و نشان دهنده اثر خطی متغیرهای کمکی و
 $f_r(\cdot)$ 
توابع ناشناخته‌ای هستند که برای در نظر گرفتن همبستگی‌های فضایی، همبستگی‌های خوشه‌بندی یا همبستگی‌های زمانی و اثرات غیر خطی متغیرهای کمکی می‌باشند، و $w_{ir}$ وزن‌های معلوم و $\varepsilon_i$ ها اثرات تصادفی یا خطاها می‌باشند.\\
 هدف تحلیل بیزی تقریبی مدل (\ref{eq2}) و به‌دست آوردن برآوردهای بیزی پارامترهای این مدل می‌باشد.  با در نظر گرفتن پیشین گاوسی مبهم برای بردار پارامترهای رگرسیونی $\bm{\beta}$ 
 و یک پیشین گاوسی برای بردار $(b_1, ..., b_k)$ و $\bm f_r(\cdot)$،  پیشگوی $\eta_{ik}$ به‌صورت گاوسی توزیع می‌شود و با تعریف بردار $\bm{x}$ به‌صورت 
 $$\bm{x}=\{\{\bm f_r(\cdot)\}, \{\beta_j\}, \{\eta_{ik}\}, \{b_k\}\},$$
 آن‌گاه می‌توان بردار پنهان $\bm{x}$ را با یک میدان پنهان گاوسی با ماتریس دقت  $Q^{-1} (\bm{\theta_1})$، به‌صورت 
%\begin{equation*}
$\bm{x}\mid\bm{\theta_1}\sim N(0, Q^{-1}(\bm{\theta_1}))$
%\end{equation*}
مدل‌بندی نمود که در آن $\bm{\theta_1}$ بردار ابر پارامترها می‌باشد و پس از در نظر گرفتن پیشین‌ها برای اعضا  $\bm{x}$ مشخص می‌شوند.\\
اکنون فرض می‌شود $y_i=((t_i, \delta_i), i=1, ..., m)$ داده‌های بقا، به‌شرط میدان پنهان $\bm{x}$ مستقل می‌باشند. بنابراین توزیع پسین 
\begin{equation*}
\pi(\bm{x},\bm{\theta}\mid \bm{y})\propto\pi(\bm{\theta})\pi(\bm{x}\mid\bm{\theta})\Pi_{i=1}^{m}\pi(y_i\mid \bm{x},\bm{\theta}),
\end{equation*}
می‌شود. که در آن $\bm{\theta}=(\bm{\theta_1}, \bm{\theta_2})$ بردار ابر پارامترهای مربوط به متغیرهای پنهانی $(\bm{\theta_1})$ و بردار ابر پارامترها مربوط به درستنمایی $(\bm{\theta_2})$ می‌باشند. این توزیع پسین معمولاً شکل بسته‌ای ندارد، بنابراین از روش‌های تکرارشونده مثل الگوریتم مونت کارلوی زنجیر مارکوفی برای به‌دست آوردن برآورد پارامترها استفاده می‌شود که بسیار زمان‌بر و گاهی اوقات همگرایی به راحتی به‌دست نمی‌آید. 
توزیع‌های پسین حاشیه‌ای متغیرهای پنهان و  ابر پارامترها به‌صورت 
\begin{eqnarray}
\label{eq3}
\pi(x_i\mid\bm y)=\int \pi(x_i\mid \bm\theta,\bm y) \pi(\bm\theta\mid\bm y)d\bm\theta,
\end{eqnarray}
\begin{eqnarray}
\label{eq4}
\pi(\theta_j \mid\bm y)=\int \pi(\bm\theta\mid\bm y)d\bm\theta_{-j},
\end{eqnarray}
می‌باشند. 
در روش بیز تقریبی ابتدا با به‌کار بردن تقریب لاپلاس برای توزیع‌های $\pi(\bm\theta\mid \bm{y})$ و $\pi(x_i\mid\bm{\theta}, \bm{y})$ تقریب‌هایی به‌صورت $\tilde{\pi}(x_i\mid \bm{y})$ و $\tilde{\pi} (\bm{\theta_j}\mid\bm{y})$ به‌دست آورده می‌شود و سپس انتگرال (\ref{eq3}) به‌صورت
\begin{equation*}
\pi(x_i\mid \bm{y})\approx\sum_{q=1}^{Q} \tilde{\pi}(x_i\mid \bm{\theta_q}, \bm{y})\tilde{\pi}(\bm{\theta_q}\mid \bm{y})\bigtriangleup_q
\end{equation*}
تقریب زده می‌شود.
رو و هلد 
\cite{Rue}
از تقریب لاپلاس استفاده کردند و نشان دادند $\pi(\bm x\mid\bm\theta_q,\bm y)$ دارای توزیع تقریبی نرمال است و طبق خاصیت بسته بودن توزیع نرمال نسبت به حاشیه‌سازی، $\tilde{\pi}(x_i \mid\bm{\theta_q}, \bm{y})$ دارای توزیع تقریبی نرمال می‌باشد و برای به‌دست آوردن تقریبی از $\pi(\bm\theta\mid\bm y)$ از الگوریتم  توضیح داده شده در بخش (\ref{rt1}) استفاده می‌شود، یعنی ابتدا مد تابع $\log \tilde{\pi}(\bm{\theta}\mid \bm{y})$ با روش بهینه سازی به‌صورت $\theta^\ast$ تعیین می‌شود. ماتریس هسیان\LTRfootnote{Hessian Matrix} $H=\dfrac{d^{2}\log\tilde{\pi}(\bm{\theta}\mid\bm{y})}{d\theta_{i}d\theta_{j}}$ محاسبه و عکس آن به‌صورت $H^{-1}=V\Lambda \acute{V}$ تجزیه می‌شود، که در آن $V$ ماتریس بردارهای ویژه و $\Lambda$ ماتریس قطری مقادیر ویژه آن است. مبدأ مختصات را به مد $\theta^\ast$ انتقال داده، فرمول مختصات در مبدأ $\theta^\ast$ به‌صورت
%\begin{equation*}
$\theta(t)=\theta^\ast + V\Lambda^{-1}t$
%\end{equation*} 
که $t$ همان مقادیر ویژه استاندارد هستند تعریف می‌شود. با شروع از مبدأ مختصات جدید روی هر یک از محورها نقاطی به فاصله مقادیر صحیح $\delta_\pi$ به گونه‌ای اختیار می‌شوند که شرط 
%\begin{equation*}
$\log \tilde{\pi}(\bm{\theta}(0)\mid\bm{y})-\log\tilde{\pi}(\bm{\theta}(t)\mid\bm{y})<\delta_t$
%\end{equation*}
برقرار باشد، سپس به‌طور مشابه نقاط درون صفحات نیز تعیین می‌شوند. با به‌کار بردن این الگوریتم $\bm{\theta}_q$ هایی از توزیع  $\pi(\bm{\theta}\mid\bm{Y})$ تولید می‌شوند که می‌توان از آن‌ها برای تقریب توزیع $\pi(\bm{\theta}\mid\bm{Y})$ استفاده کرد. رو و همکاران\cite{{Ruee}} این روش تقریبی انتگرال (\ref{eq3})  و (\ref{eq4})  را تقریب لاپلاس آشیانی جمع بسته
\lr{INLA}
نامیدند. در ادامه روش ارائه شده در یک مثال کاربردی مورد بررسی قرار می‌گیرد.
\section{مثال کاربردی} 
داده‌های مورد بررسی مشاهدات مربوط به  $m=1043$ بیمار مبتلا به سرطان خون بین سال‌های $1982$ الی $1998$ در شمال غرب انگلستان در $n_s=24$ ناحیه می‌باشند که حدود $16$ درصد موارد سانسورشده از سمت راست هستند و $t=1,\ldots,4977$.
 هندرسون\LTRfootnote{Henderson} و همکاران
 \cite{Henderson}
، نیب و فاهرمیر\cite{Kneib} به مطالعه این داده‌ها پرداخته‌اند. برای تحلیل این داده‌ها از نرم افزار $R$ استفاده شده است که کدهای آن در پیوست \ref{appendix1} موجود می‌باشد.\\
متغیرهای کمکی برای هر بیمار $i$  شامل سن بیمار $(age_i)$، تعداد سلول‌های سفید خون $(wbc_i)$، شاخص محرومیت تاونزند\LTRfootnote{Townsend}، 
$(tpi_i)$
جنسیت بیمار$(sex_i)$ و منطقه سکونت$(s_i)$ می‌باشند. نمودار داده‌ها در شکل \ref{fig1-3} مشخص شده است. 
\begin{figure}
\begin{center}
\includegraphics[scale=.5]{fig110.png}
\end{center}
\caption{موقعیت داده‌های سرطان خون در شمال غرب انگلستان \label{fig1-3}}
\end{figure}
در این‌جا زمان در $K=20$
 بازه زمانی برابر افراز می‌گردد و چهار مدل کاکس با لگاریتم ثابت تکه‌ای خطر پایه با همبستگی فضایی، بدون همبستگی فضایی، بدون اثرات غیرخطی و با افزایش اثرات غیرخطی در نظر گرفته می‌شود و سپس با معیار مدل‌گزینی 
\lr{DIC}
مورد مقایسه قرار می‌گیرند. برای افرازهای زمانی $K=10$ و $K=50$ محاسبات تکرار و تفاوت محسوسی در نتایج ایجاد نشد.
در هر چهار مدل $b_k$ فرآیند قدم زدن تصادفی مرتبه اول با توزیع نرمال به شکل $b_{k+1}- b_{k}\sim N(0,{\tau_b}),~~k=1,\ldots,19$ فرض شد.\\
 در مدل اول اثرات سن و جنس به‌صورت خطی و اثرات تعداد سلول‌های سفید و شاخص تاونزند به‌صورت غیرخطی به‌صورت قدم زدن تصادفی مرتبه دو با پارامتر دقت $\tau_{tpi}$ و $\tau_{wbc}$ بدون اثر فضایی در نظر گرفته شد، یعنی
\begin{eqnarray*}
\pi(f_2\mid\tau_{tpi})\propto\tau_{tpi}^{\dfrac{(n_{tpi}-2)}{2}}\exp\lbrace -\dfrac{1}{2}\tau_{tpi}\Sigma_{i=3}^{n_{tpi}}(f_{2i}-2f_{2(i-1)}+f_{2(i-2)})^2\rbrace
\end{eqnarray*}
\begin{eqnarray*}
\pi(f_1\mid\tau_{wbc})\propto\tau_{wbc}^{\dfrac{(n_{wbc}-2)}{2}}\exp\lbrace -\dfrac{1}{2}\tau_{wbc}\Sigma_{i=3}^{n_{wbc}}(f_{1i}-2f_{1(i-1)}+f_{1(i-2)})^2\rbrace
\end{eqnarray*}
بنابراین مدل $1$ به‌صورت
\begin{eqnarray*}
M_1:\eta_{ik}=\beta_0+\beta_1 {sex}_i+\beta_{2} {age}_i+f_1({wbc}_i)+f_2({tpi}_i)+b_k
\end{eqnarray*}
 تعریف می‌شود. نتایج تحلیل این مدل در جدول \ref{tab1} ارائه شده است.
زمان اجرای برنامه‌های کامپیوتری برای این مدل حدود 9 ثانیه می‌باشد. در مدل دوم اثرات سن، جنس و تعداد سلول‌های سفید به‌صورت اثرات خطی و اثر شاخص تاونزند به‌صورت غیرخطی قدم زدن تصادفی مرتبه دو با پارامتر دقت $\tau_{tpi}$ و اثر فضایی به‌صورت یک میدان تصادفی مارکوفی گاوسی با پارامتر دقت $\tau_s$ در نظر گرفته شد.
\begin{table}[htp] 
\vspace{-5mm}
%\label{tab1}
\scriptsize
\begin{center}
\begin{small}
\caption{\label{tab1} برآورد بیز تقریبی پارامترهای مدل 1، (میانگین پسین، انحراف معیار و چندک‌ها) }
\begin{tabular}{cccccccccc}
\hline 
پارامتر   &  میانگین  & انحراف معیار &  &\multicolumn{5}{c}{چندک‌ها}  \\ 
%&&&&   &\multicolumn{5}{c}{$n$} & &  \\
\cline{4-9}
&&&&\multicolumn{1}{c}{$0/025$} & & \multicolumn{1}{c}{$0/5$}  & & \multicolumn{1}{c}{$0/975$}\\    
%\cline{4-8} 
% \cline{7-8}   
%\vspace{2mm}
%مدل        & $\alpha$  &~~   & $std$ & $$  &~~ & $W$ &  $K$ &~~ & RSS \\   
\hline
$\beta_0$                  &$-9.487$ & $0.2572$ &&  $-10.0186$ &&$-9.4772$& &  $ -9.0074$   \\
$\beta_{1}$            &$0.0602$&$ 0.0681$&& $-0.0734$&&$ 0.0602 $&&$ 0.1940$  \\
$\beta_{2}$            &$0.0323$& $0.0021$&& $0.0282$&& $0.0323$&& $ 0.0365$     \\
$\tau_{wbc}$           &$31543.77$& $23701.63$ && $ 4924.76$&& $25672.26$&& $93006.89$   \\
$\tau_{tpi}$              &$30383.51$&$ 21758.29$&&$5347.08$&&$25144.78$&&$86427.80$   \\  
$\tau_{b}$                &$1250.75$ &$587.03$ && $ 436.22$&& $1145.53$&& $2684.23$  \\
 \hline
\end{tabular}
\end{small}
\end{center}
\vspace{-5mm}
\end{table}

بنابراین مدل $2$ به‌صورت 
\begin{eqnarray*}
M_2:\eta_{ik}=\beta_0+\beta_{1} {sex}_i+\beta_{2} {age}_i+\beta_{3} {wbc}_i+f_1(tpi)+f_{spat}(s_i)+b_k
\end{eqnarray*}
می‌باشد. نتایج در جدول \ref{tab2} ارائه شده است و زمان اجرای برنامه‌های کامپیوتری حدود 8 ثانیه می‌باشد.\\
\begin{table}[htp]
\vspace{-5mm}
%\label{tab2}
\scriptsize
\begin{center}
\begin{small}
\caption{\label{tab2} برآورد بیز تقریبی پارامترهای مدل 2، (میانگین پسین، انحراف معیار و چندک‌ها) }
\begin{tabular}{cccccccccc}
\hline 
پارامتر   &  میانگین  & انحراف معیار &  &\multicolumn{5}{c}{چندک‌ها}  \\ 
%&&&&   &\multicolumn{5}{c}{$n$} & &  \\
\cline{4-9}
&&&&\multicolumn{1}{c}{$0/025$} & & \multicolumn{1}{c}{$0/5$}  & & \multicolumn{1}{c}{$0/975$}\\    
%\cline{4-8} 
% \cline{7-8}   
%\vspace{2mm}
%مدل        & $\alpha$  &~~   & $std$ & $$  &~~ & $W$ &  $K$ &~~ & RSS \\   
\hline
$\beta_0$                     &$-10.1638$&$ 0.2516$&&$-10.6855 $ &&$-10.1543$&& $-9.6963$   \\
$\beta_{1}$               &$0.065$& $0.0686$&&    $-0.0580$&&   $0.0765$&&    $ 0.2112$ \\
$\beta_{2}$               &$0.0345$& $0.0022$ &&    $0.0302$&&  $ 0.0345$&&     $0.0389$ \\
$\beta_{3}$              &$0.0034$& $0.0005$&&    $ 0.0025$&&   $0.0034$&&     $0.0043$ \\
$\tau_{tpi}$                  &$8791.15$&  $9767.09$&&   $391.724$&&   $5734.94$&& $35389.44$   \\
$\tau_{s}$                    &$15.38$&  $11.40$&&  $3.902$&&  $12.23$&&  $ 45.48$            \\    
$\tau_{b}$                    &$1292.06$&   $605.14$&&   $443.950$ &&  $1186.82$&&  $2762.82$    \\
 \hline
\end{tabular}
\end{small}
\end{center}
\vspace{-2mm}
\end{table}

در مدل 3 تنها اثرات سن، جنس و تعداد سلول‌های سفید به‌صورت خطی در مدل در نظر گرفته شد.
\begin{eqnarray*}
M_3:\eta_{ik}=\beta_0+\beta_{1} {sex}_i+\beta_{2} {age}_i+\beta_{3}{wbc}_i+b_k
\end{eqnarray*}
نتایج تحلیل این مدل در جدول \ref{tab3} آورده شده است و کل زمان اجرای برنامه کامپیوتری حدود 2 ثانیه می‌باشد.
\begin{table}[htp]
\vspace{-5mm}
%\label{tab3}
\scriptsize
\begin{center}
\begin{small}
\caption{\label{tab3} برآورد بیز تقریبی پارامترهای مدل 3، (میانگین پسین، انحراف معیار و چندک‌ها) }
\begin{tabular}{cccccccccc}
\hline 
پارامتر   &  میانگین  & انحراف معیار &  &\multicolumn{5}{c}{چندک‌ها}  \\ 
%&&&&   &\multicolumn{5}{c}{$n$} & &  \\
\cline{4-9}
&&&&\multicolumn{1}{c}{$0/025$} & & \multicolumn{1}{c}{$0/5$}  & & \multicolumn{1}{c}{$0/975$}\\    
%\cline{4-8} 
% \cline{7-8}   
%\vspace{2mm}
%مدل        & $\alpha$  &~~   & $std$ & $$  &~~ & $W$ &  $K$ &~~ & RSS \\   
\hline
$\beta_0$                     &$-10.1029$&$0.2499 $&&$-10.6214 $ &&$-10.0932$&& $-9.6391$   \\
$\beta_{1}$               &$0.0668$& $0.0678$&&    $-0.0662$&&   $0.0668$&&    $ 0.1999$ \\
$\beta_{2}$               &$0.0321$& $0.0021$ &&    $0.0279$&&  $ 0.0320$&&     $0.0362$ \\
$\beta_{3}$              &$0.0033$& $0.0004$&&    $ 0.0024$&&   $0.0033$&&     $0.0041$ \\
$\tau_{b}$                    &$1248.70$&   $584.26$&&   $436.42$ &&  $1144.42$&&  $2674.25$    \\
 \hline
\end{tabular}
\end{small}
\end{center}
\vspace{-2mm}
\end{table}

در مدل 4، اثرات سن، جنس و تعداد سلول‌های سفید به‌صورت خطی و اثر شاخص تاونزند به‌صورت غیرخطی قدم زدن تصادفی مرتبه دو با پارامتر دقت $\tau_{tpi}$ و اثر فضایی به‌صورت مدل $ (BYM)$\LTRfootnote{Besag, York and Mollie Model} که توسط بسج\LTRfootnote{Besag} و همکاران
\cite{Besag}
معرفی شده است، در نظر گرفته شد. در مدل $BYM$ دو نوع اثر تصادفی وارد مدل می‌شود، اثر تصادفی فضایی و اثر ناهمگنی به‌طوری که اثر فضایی، همبستگی فضایی مشاهدات را در مدل بیان می‌کند و اثر ناهمگنی برای این فرض وارد مدل می‌شود که برآوردهای به‌دست آمده بین نواحی مختلف مستقل از یک‌دیگر می‌باشند. یعنی دو اثر همبستگی فضایی درونی و استقلال برونی را با مدل $BYM$ می‌توان در نظر گرفت. مدل $4$ به‌صورت 
\begin{equation*}
M_4:\eta_{ik}=\beta_0+\beta_{1} {sex}_i+\beta_{2} {age}_i+\beta_{3}{wbc}_i+f_1(tpi)+f_{iid}(H_i)+f_{spat}(s_i)+b_k
\end{equation*}
می‌باشد. نتایج مدل 4، در جدول \ref{tab4} ارائه شده است. و زمان اجرای برنامه‌های کامپیوتری حدود 15 ثانیه می‌باشد. پارامترهای دقت اثر فضایی و اثر ناهمگنی به ترتیب $\tau_{s}$ و $\tau_{H}$ می‌باشند.\\
\begin{table}[htp]
\vspace{-5mm}
%\label{tab4}
\scriptsize
\begin{center}
\begin{small}
\caption{\label{tab4} برآورد بیز تقریبی پارامترهای مدل 4، (میانگین پسین، انحراف معیار و چندک‌ها) }
\begin{tabular}{cccccccccc}
\hline 
پارامتر   &  میانگین  & انحراف معیار &  &\multicolumn{5}{c}{چندک‌ها}  \\ 
%&&&&   &\multicolumn{5}{c}{$n$} & &  \\
\cline{4-9}
&&&&\multicolumn{1}{c}{$0/025$} & & \multicolumn{1}{c}{$0/5$}  & & \multicolumn{1}{c}{$0/975$}\\    
%\cline{4-8} 
% \cline{7-8}   
%\vspace{2mm}
%مدل        & $\alpha$  &~~   & $std$ & $$  &~~ & $W$ &  $K$ &~~ & RSS \\   
\hline
$\beta_0$                     &$-10.1486$&$0.2553 $&&$-10.6772 $ &&$-10.1392$&& $-9.6732$   \\
$\beta_{1}$               &$0.0676$& $0.0687$&&    $-0.0671$&&   $0.0676$&&    $ 0.2024$ \\
$\beta_{2}$               &$0.0341$& $0.0022$ &&    $0.0298$&&  $ 0.0341$&&     $0.0385$ \\
$\beta_{3}$              &$0.0034$& $0.0005$&&    $ 0.0025$&&   $0.0034$&&     $0.0042$ \\
$\tau_{tpi}$                  &$7155.70$&  $7387.98$&&   $370.751$&&   $4910.48$&& $27554.61$   \\
$\tau_{H}$                  &$32.57$&  $22.11$&&   $9.019$&&   $26.68$&& $90.69$   \\
$\tau_{s}$                  &$1762.40$&  $1795.47$&&   $102.103$&&   $1222.28$&& $6504.98$   \\
$\tau_{b}$                    &$1283.54$&   $606.37$&&   $445.486$ &&  $1173.76$&&  $2767.69$    \\
 \hline
\end{tabular}
\end{small}
\end{center}
\vspace{-2mm}
\end{table}
لازم به‌ذکر است که در هر چهار مدل، اثر فضایی به‌صورت فرآیند اتورگرسیو شرطی ذاتی\LTRfootnote{Intrinsic Conditiond Aut-
Regressive} توسط یک میدان تصادفی مارکوفی گاوسی مدل‌بندی شده است یعنی
\begin{eqnarray*}
[(f_{spat})_i\mid(\bm f_{spat})_{-i},\tau_s]\sim N(\dfrac{1}{n_i^s}\Sigma_{j\in\partial_i}(f_{spat})_j, \dfrac{1}{n_i^s\tau_s})
\end{eqnarray*}
که در آن $\partial_i$ مجموعه همسایگی مربوط به موقعیت $i$ام و $n_i^s$ تعداد آن‌ها، برای تمام پارامترهای دقت توزیع پیشین مبهم گاما با پارامترهای $(1,0.001)$ و برای پارامترهای رگرسیونی توزیع پیشین مبهم نرمال با واریانس بزرگ در نظر گرفته شده است. برای مقایسه چهار مدل از معیار $DIC$ استفاده شد. همان‌طور که در فصل اول اشاره شد برای مقایسه بیزی مدل‌ها معرفی شده است و به‌صورت 
\begin{equation*}
DIC=\bar{D}(\theta)+P_D=D(\bar{\theta})+2P_D
\end{equation*}
تعریف می‌شود، که در آن $\bar{D}(\theta)$ میانگین پسین انحرافی می‌باشد که میزان برازش مدل را بیان می‌کند و $P_D$ میزان پیچیدگی مدل را نشان می‌دهد. $P_D$ برابر تعداد پارامترهای موثر مدل نامیده می‌شود که مقدار آن برابر با $P_D=\bar{D}(\theta)-D(\bar{\theta})$ می‌باشد. هر مدلی که کمترین $DIC$ را دارد بهترین مدل انتخاب می‌شود. نتایج در جدول \ref{tab5} ارائه شده است.\\
\begin{table}[htp] 
\vspace{-5mm}
%\label{tab5}
\scriptsize
\begin{center}
\begin{small}
\caption{\label{tab5}مقایسه مدل‌ها }
\begin{tabular}{cccccccccc}
\hline 
$Model$   &  $P_D$  & $DIC$ & 

 & \multicolumn{1}{c}{$log-likelihood$}\\      
%\vspace{2mm}
%مدل        & $\alpha$  &~~   & $std$ & $$  &~~ & $W$ &  $K$ &~~ & RSS \\   
\hline
$M_1$                  &$16.49$ & $5381.11$ &&   $ -2656.29$   \\
$M_2$            &$24.51$&$ 5354.77$&& $ -2767.39$  \\
$M_3$            &$11.73$& $5394.86$&&  $ -26748.6$     \\
$M_4$           &$26.67$& $5356.67$ &&  $-2744.10$   \\
 \hline
\end{tabular}
\end{small}
\end{center}
\vspace{-2mm}
\end{table}

با توجه به نتایج جدول \ref{tab5}، مدل دوم و چهارم از دو مدل اول و سوم مناسب‌تر به‌نظر می‌رسند و بین این دو مدل با مقایسه تمامی موارد، ابتدا با توجه به مقادیر $P_D$ مدل چهارم به‌ عنوان مدل مناسب انتخاب می‌گردد. با انتخاب مدل چهارم و نتایج جدول \ref{tab4} به تحلیل این مدل می‌پردازیم.
\section{بحث و نتیجه‌گیری}
 به‌دلیل پیچیدگی این مدل‌ها، نداشتن شکل بسته برای توزیع پسین و برخی توزیع‌های شرطی کامل، محاسبات و اجرای برنامه‌های کامپیوتری بسیار زمان‌بر می‌باشند. تاکنون روش‌های محاسباتی مختلفی جایگزین روش‌های بیزی معمولی پیشنهاد شده است. نیب و فاهرمیر
 \cite{Kneib}
 با به‌کار بردن میدان تصادفی گاوسی مارکوفی و درستنمایی جریمه شده به تحلیل مدل کاکس با تابع خطر پایه نمایی تکه‌ای پرداختند 