\documentclass[11pt,a4paper]{article}
%\usepackage{palatino}
\usepackage{graphicx}
\usepackage{rotating}
\usepackage{amsmath}
\usepackage{fancyvrb}
\setlength{\parindent}{0pt}
\setlength{\parskip}{6pt}
\setlength{\arrayrulewidth}{1pt}
\renewcommand{\textfraction}{0.05}
\setlength{\textwidth}{16cm}
\setlength{\textheight}{24.5cm}
\setlength{\oddsidemargin}{0pt}
\setlength{\topmargin}{0pt}
\setlength{\headheight}{0pt}
\setlength{\headsep}{0pt}
\newcommand{\bm}[1]{\boldsymbol{#1}}
\newcommand{\var}[1]{\mbox{var}\left({#1}\right)}
\newcommand{\covar}[1]{\mbox{cov}\left({#1}\right)}
\newcommand{\corr}[1]{\mbox{corr}\left({#1}\right)}
\newcommand{\vcite}[2][~]{#1$[#2]$}
%\newcommand{\vcite}[2][]{#1$^{\bm{#2}}$}
\linespread{1.5}
\renewcommand{\baselinestretch}{1.5}

\newcommand{\D}{\displaystyle}  %colocar no numerador e denominador das frações em tabelas para fica com o tamanho normal

\usepackage{fancyvrb}
\DefineVerbatimEnvironment{Rin}{Verbatim}{fontshape=sl}
\DefineVerbatimEnvironment{Rout}{Verbatim}{}
\newcommand{\R}{{\normalfont\textsf{R}}{}}
\newcommand{\code}[1]{\textsl{\texttt{#1}}}
\newcommand{\codeout}[1]{\texttt{#1}}


\title{Comparing modelling strategies for longitudinal data: \\the effect of Vitamin A in childhood diarrhoeal disease}
\author{ V. L. Andreozzi$^\dag$ \& M. S. Carvalho \\[3pt]
{\small \em National School of Public Health - Fiocruz, Brazil.}\\[6pt]
M. L. Barreto\\[3pt]
{\small \em Federal University of Bahia, Brazil.}
}
\date{}

\begin{document}

\maketitle
\renewcommand{\footnoterule}{}
\renewcommand{\thefootnote}{}
\footnotetext{\hspace*{-0.23in}\dag {\em Correspondence to:}
Dr V.L. Andreozzi, National School of Public Health, Fiocruz. Rua Leopoldo Bulh\~{o}es 1480/827, Manguinhos, Rio de Janeiro, RJ, Brazil, Cep: 21041-210. {\tt
E-mail:valeska@fiocruz.br}}


{\linespread{1}
\begin{abstract}

\noindent \textbf{Introduction}:The choice of modelling strategies can improve or impair the best designed studies. It is well known the role of recent developed techniques, such as generalised estimating equations and random effects models, on analysing longitudinal data. On the other hand, late updates on survival analysis models allow the inclusion of multiple time-to-event periods and random effects. In this paper we compare the use of these methods in a longitudinal study on the effects of Vitamin A in preventing childhood diarrhoeal disease. 

\textbf{Methods}: Data comes from a randomised, double-blind, placebo-controlled clinical trial where over 1000 children aged 6-48 months were assigned to receive either vitamin A or placebo, every 4 months over 1 year. Demographic, clinical, nutritional and socio-economic characteristics were collected at the baseline. The subjects were visited at home three times a week, to register the number of liquid or semi-liquid motions a day. Three response variables were considered: number of liquid or semi-liquid motions per week, numbers of days with diarrhoea per week, in the longitudinal models; numbers of days without diarrhoea in the survival models. 

% \textbf{Results}: The survival models are more comparable to the longitudinal models for the number of days with diarrhoea, since both of them model the occurrence of diarrhoea episodes. For the two ways of looking at the problem, no significant effect was found for the vitamin A supplementation. The only significant effect was found when modelling the severity of the disease, using the number of liquid or semi-liquid motions per week. The most important covariates were the age of the child, as expected, and the number of domestic appliances (refrigerator, radio, TV). 

\textbf{Conclusion}: The comparison of the models showed reasonable consistency; even the differences found can be attributed to different measures and ways of looking at the problem. The major benefit trying those various models is the ability of detecting small change in point and confidence interval estimates, including independent variable of interest and all potential risk factors. The variety of statistical models is increasing sharply with the development of statistical tools and availability of software. Unfortunately the choice among them is usually more dependent on previous knowledge then on the real question addressed by the study. Quite often data like this one is only analysed using simple summary measures, in classical models easily adjusted in common statistical packages, causing loss of information so hardly gathered. There is a strong need in the epidemiology environment of incorporation of an entire set of statistic techniques dealing with dependent data, particularly from longitudinal studies. 

\paragraph{Keywords: Longitudinal Study, Marginal models, Random effects models, Survival Models}
\end{abstract}
}

\section{Introduction}

Several statistical methods are available to estimate measures of association between risk factors and disease outcome in epidemiological studies. The choice of the model to fit the data is guided by the type of response variable, desired measure of association, possible non-linearity, and, as often, by the previous knowledge of the researcher in charge of the analysis. In longitudinal studies the main focus of interest is on changes over time, either evaluated as time-to-event or as repeated measures taken during the follow-up period. 

One of the most widespread models in the medical scientific literature is the semi-parametric proportional hazards survival model, presented by Sir David Cox in 1972 (ref - DR Cox (1972) Regression models and life tables Journal of the Royal Statistical Society Series B 34:187-220.). A brief search in the bibliographic database Medline, using the words "Cox" and "survival", limited for the period between 2002 and 2006, returned 7041 papers. However, the inclusion of the terms "multiple" or "repeated" and "events" reduced the amount to a mere 60. In spite of being avaible for years, the use of extended Cox models is still limited among epidemiologists and other health data analysts. On the other hand, the search based on the words "random", "effects" and "longitudinal" returned just 309 papers, some 33\% in statistical journals. 

In this paper we will compare different statistical approaches -- using survival and repeated measures models -- to estimate the effect of vitamin A supplementation on children diarrhoea in a longitudinal study. According to the World Health Organisation \cite{Feachen148,whovita}, supplementation of vitamin A have shown a reduction of the frequency of moderate and severe diarrhoea. Diarrhoeal disease is widely recognised as a leading cause of morbidity and mortality, with estimated 3.3 million children deaths per year in the last decade \cite{bern142}. 

The determinants of diarrhoea mortality are well known. In brief, the biological risk factors comprise maternal factors (age at birth, parity and birth intervals); environmental contamination (intensity of household crowding, water contamination, household food contamination or potential fecal contamination); nutrient deficiency (nutrient availability to the infant or to the mother during pregnancy and lactation). The socio-economic determinants, generally operating through the previous factors, are grouped into three broad categories: individual ones (such as mother's educational level, beliefs about disease causation, food preferences); household factors (income/wealth effects acting upon food availability, quality of water supply, housing conditions, fuel/energy availability, means to purchase what is necessary for the daily practice of hygienic/preventive care, access to information); and community-level factors, including the ecological setting (climate, temperature, altitude, season, rainfall) and political, economy and health system variables \cite{masuy96}.

In this frame, the role of vitamin A supplementation on diarrhoea mortality and incidence was investigated by Barreto \textit{et al} \cite{barreto116}, in a longitudinal randomised, double-blind, placebo-controlled community trial in the North-East of Brazil. The main method used in the original paper was a Poisson regression model, with the number of diarrhoea episodes per child-day as dependent variable and the intervention group (vitamin A or placebo) as variable of interest. In this paper, the vitamin A supplement presented a significant protective efect (rate ratio = 0.80 [0.65--0.98]) only for serious diarrhoea episodes. Summarize data of each individual into one statistics, such as incidence of diarrhoea episode per child, reduces the analysis of longitudinal data to the analysis of independent observations, for which classical procedures are available. Although this is a valid procedure, the main disadvantage is the impossibility to study the disease evolution within-child, not counting the amount of unused information collected on various individual and environment covariates. A latter publication \cite{Valpaper}, using random effects model, estimated a risk for diarrhoea infection 1.57 times greater for a child administered a placebo as opposed to vitamin A. However, a temporal effect, showing a decreasing daily probability of diarrhoea for both groups during the trial, was observed, suggesting the use of some further adjustment. 

To thoroughly explore this problem, we fitted five different models: three using repeated measures techniques and two modelling the time-to-event. We do not intend to state which of the five is the best one, as there is no goodness-of-fit measure to compare such different approaches. Besides, the main idea is to discuss the robustness of the main estimates, incentivating, whenever possible, the use of different modelling approaches. In section \ref{sec:longitmodels} we present the repeated measures approach, marginal models in subsections \ref{subsec:marginal} and random effects in subsection \ref{subsec:rem}). In section \ref{sec:survmodels} we describe the Cox survival models with extensions to deal with multiple independent events (subsection \ref{subsec:ag}) and with random effects (subsection \ref{subsec:frail}. The following sections present the dataset (section \ref{sec:data}) and the modelling strategies (\ref{sec:met}). Results and discussion are presented respectively in sections \ref{sec:results} and \ref{sec:discussion}.

\section{Longitudinal Models} \label{sec:longitmodels}

There are two main perspectives to analyse longitudinal data: the marginal models and the random effects models. The first one is focused on modelling the mean response over time across an observed population, considering the correlation among the repeated observations of each individual as a nuisance parameter. The random effects approach, on the other hand, treats the repeated measures correlation as an unmeasured factor characterizing individuals and condition the effect of covariates on each individual trajectory. 


\subsection{Marginal Models} \label{subsec:marginal}
Marginal models extend the standard generalized linear model (GLM) framework, appropriate for independent observations, to include adjustment for within subject correlation\vcite{3,7}. A relationship between the explanatory variables and the \textit{marginal} mean response is specified, thus averaging out any unobserved subject-specific effects. In this approach the within-subject correlation function is separately specified from the mean relationship. The parameters affecting the mean response therefore have the same population averaged interpretation as in the standard GLM. 

In the general marginal model framework the marginal mean response for the $i^{th}$ subject at time $t$, $\mu_{it}$, is related to a vector of $p$ subject-specific (possibly time-varying) covariates, $\bm{x}_{it}$, via a specification: 

\begin{equation}
g(\mu_{it})=\bm{x}_{it}^\prime\bm{\beta},
\end{equation}
and the marginal variance of the response, $\sigma^2_{it}$, is related to its mean, such that

\begin{equation}
\sigma^2_{it}=\phi v(\mu_{it})
\end{equation}
Here $\bm{\beta}$ is a vector of unknown fixed population effects, $\phi$ is an unknown overdispersion (or scale) parameter, and the link function $g(\cdot)$ and the variance function $v(\cdot)$ are assumed known. 

Where the specification of the marginal model departs from that of the standard GLM is that, while independence continues to be assumed for observations taken on \textit{different} subjects, correlations between observations taken on the \textit{same} subject are no longer assumed to be zero, but rather characterised by a further known correlation function

\begin{equation}
Corr(Y_{it},Y_{ik})=c(\mu_{it},\mu_{ik},\alpha)
\end{equation}
$c(\cdot)$ which will in turn involve a number of further unknown parameters  $\bm{\alpha}$. 

The parameters of marginal models may be estimated using a generalized estimating equation (GEE) approach. This methods decompose the variance of the response variable in three parts, named: dispersion parameter ($\phi$), a diagonal variance matrix explained by the mean ($A$) and the correlation structured ($R$) as follows:

\begin{equation}
V_i=\phi A_i^{1/2} \times R(\alpha) \times A_i^{1/2} 
\end{equation}
The correlation structured is assumed to be the same for all subjects. Various alternative specifications for the within-subject correlation structure are possible, such as, independente ($corr(Y_{ij},Y_{ik})=0$), exchangeable ($corr(Y_{ij},Y_{ik})=\alpha$), autoregressive structures ($corr(Y_{ij},Y_{ik})=\alpha^t$ for $t=0,1,2,\cdots$) and unstructured ($corr(Y_{ij},Y_{ik})=\alpha_{jk}$). By appropriate choice of the link, variance functions and correlation structures, a wide class of response variables can be modelled within this framework. 

Details are given in  Zeger \& Liang \cite{Zeger77} and Liang \& Zeger \cite{Zeger56}, who demonstrate that this method results in consistent estimates of $\bm{\beta}$ even when the correlation structure in the model is misspecified and who provide expressions for the standard errors of the $\bm{\beta}$ estimates which are robust to any such misspecification. GEE methods for marginal models can handle different number of observations on each subject\vcite{16} including values missing completely at random, and unbalanced designs;  this being relevant to the varying follow-up periods in the clinical trial considered here. A variety of software is available to implement GEE methods, here we use the public available GEEPACK package \cite{geepack} developed for marginal models within the public domain statistical package R \cite{r}.

\subsection{Random Effects Models} \label{subsec:rem}

There is a long tradition of using random effects to handle heterogeneity in the analysis of clinical trials, particularly where longitudinal and/or repeated measures data is involved (e.g. see Laird and Ware\cite{Laird60}; Cnaan \textit{et al}.\cite{cnaan57}). The basic idea is to model additional sources of variability in the response by introducing additional terms into the relationship between mean response and explanatory variables which are random quantities (i.e. subject to a probability distribution) rather than fixed parameters. For example, the inclusion of appropriately specified subject-specific random effects will allow for heterogeneity between subjects and also induce within-subject correlation into the  model, both of which  commonly result from unobserved covariates i.e. important dissimilarities between individuals not reflected in the explanatory variables which  were explicitly recorded in the study.

In the case of discrete longitudinal data, such as that considered here, the appropriate modelling framework for the inclusion of random effects is the Generalised Linear Mixed Model (GLMM) (e.g. see McCulloch and Searle \cite{mac2001}, Molenberghs and Verbeke (2005) (LIVRO DE DADOS DISCRETOS)). That model formulation is hierarchical involving firstly, a specification for the \textit{conditional} mean of the response given values for the random effects, and, secondly, a specification of a joint probability distribution for the random effects. More specifically, conditional on a vector of $q$ unobservable random effects $\bm{\gamma_{it}}$, the mean response  $\mu_{it}$ is specified as

\begin{equation}
g(\mu_{it})=\bm{x}_{it}^\prime\bm{\beta} + \bm{w}_{it} ^\prime \bm{\gamma}_{it}
\end{equation}
and the conditional response variance, $\sigma^2_{it}$, is related to its mean via

\begin{equation}
\sigma^2_{it}=\phi v(\mu_{it})
\end{equation}
Here, as usual,  $\bm{\beta}$ is a vector of unknown fixed population effects, $\bm{x}_{it}$ and  $\bm{w}_{it}$ are vectors of relevant subject-specific covariates (where some of the same  covariates may appear in both sets), $\phi$ is an unknown overdispersion (or scale) parameter, and the link function $g(\cdot)$ and the variance function $v(\cdot)$ are assumed known. The random effects, $\bm{\gamma}_{it}$, are then taken to have a common underlying multivariate distribution--in practice often a multivariate normal distribution with zero means and an unknown dispersion matrix with an appropriately specified structure.

The first part of this formulation is similar to that of the standard GLM, but the important distinction is that it involves the \textit{conditional} mean and variance of the response given values for the random effects. Here, a within-subject covariance structure is indirectly generated by assuming that each individual has unique random effects, distributed around the fixed population effects\vcite{20}.  In the random effects model the fixed parameters are interpreted as effects on an \textit{average subject} in the study. 

There are a range of likelihood, likelihood related or GEE methods available for parameter estimation in the GLMM (e.g. see McCulloch and Searle \cite{mac2001};  Diggle \textit{et al}. \cite{diggle8}). To perform maximum likelihood estimation, it is necessary to be able to compute the likelihood, i.e, to integrate out the random effects. This is feasible when the model has just a subject-specific intercept. However, for numerous random effects and discrete response variable likelihood methods can experience computational difficulties and a hierarchical Bayesian approach implemented via MCMC methods can be used instead (REFerences). This approach, in a nutshell, calculates the posterior distributions for the model parameters that is proportional to the product of their prior distributions and the data likelihood. The median and quantiles of the posterior distributions provide point-estimates and credible intervals. A model with random effects was fitted to this data set and published elsewhere \cite{(ref val)}. 

An extension of the model described above, allowing a non-linear time varying effect (observed in exploratory data analysis), is called Generalised Additive Mixed Model (GAMM). The model can now be written as
\begin{equation}\label{modtempo}
g(\mu_{it})=\bm{x}_{it}^\prime\bm{\beta} + \bm{w}_{it} ^\prime \bm{\gamma}_{it} + \delta_t
\end{equation}
where $\delta_t$ corresponds to some function of time. In the Bayesian framework employed here, $\delta_t$ is estimated using a non-parametric smoothing prior, where a second order random walk is a common choice (e.g.~Fahrmeir, 1994, 1998): 

\begin{equation}\label{rw2}
\delta_t = 2\alpha_{t-1} - \alpha_{t-2} + \epsilon_t
\end{equation}

We use an inverse gamma hyperprior for the unknown variance $\sigma_\alpha^2$. 
This hyperparameter determines the smoothness of $\bm{\alpha}$, which can also be estimated from the
data. We used parameters $a=0.001$ and $b=0.001$ and differences of the sensitivity study with respect to this
settings were minor. We adopt locally uniform priors for the regression coefficients $\bm{\beta}$.

%The model selection was based on the Deviance Information Criterion that combines model fit with model complexity. VAMOS USAR???? SE SIM, REFERENCIA DE DIC 



\section{Survival Models} \label{sec:survmodels}
The Cox survival model is defined as a semi-parametric model because a parametric form is assumed for the covariate effects while the baseline hazard function is treated non-parametrically. It is also called a proportional hazards model because the ratio of hazard rates of any two individuals is proportional and constant over time. However, when an individual can experience multiple events, the dependency between events within individual has to be taken into account. As in the longitudinal approach we also have two different ways to deal with multiple events: marginal models and frailty or random effects models. The former treats the population-averaged effects of covariates, while the later includes differences between individuals through an estimated random effect \cite{therneau2000}. 

\subsection{Anderson and Gill Model} \label{subsec:ag}
Depending on the nature of the recurrent event, different marginal models are adequate. One can postulate that diarrhoea episodes are mutually independent events in each individual, occuring as a non-homogenous Poisson process. This supposition means that the usual proportional hazard assumption is unaffected by earlier events occurred to the individual, as, conditionally on covariates, the number of events in non-overlap interval times are independent. This approach is called Anderson and Gill model (AG)  --- ref do AG!!!!). Mathematically, the definition of the AG model is equal to the classic Cox proportional hazard model:
$\lambda_{it}(t|\bm{x}_i)=\lambda_0(t)\exp(\bm{x}_i^{\prime}\bm{\beta})$,  where $\lambda_0(t)$ is an arbitrary baseline hazard rate, $\bm{\beta}$ is the vector of parameters. The difference between this two models is in the construction of the risk set and in the variance estimator.


\subsection{Frailty Models} \label{subsec:frail}

Let $\bm{z}_i$ be a random variable representing an unknown random effect, related to children, with unit mean and variance $\epsilon_i$. Large values of $\bm{\epsilon}_i$ reflect a greater degree of heterogeneity among children. For an individual with regressor variables vector $\bm{x}_i$ and frailty effect $\bm{z}_i$, the form often assumed for the hazard rate is $\lambda_{it}(t|\bm{z}_i)=\bm{z}_i\lambda_0(t)\exp(\bm{x}_i^{\prime}\bm{\beta})$. The standard assumption is that observations within each child share a common frailty effect $\bm{z}_i$.

Individuals or groups with frailty values greater than 1 tend to experience the event at a rate faster than under the basic Cox model. Conversely, when frailty values are less than 1, survival times tend to be longer. Note that, when the variance of $z$ approaches zero, the model reduces to the basic Cox model. The most common model proposed in the literature for the frailty effect is the one-parameter gamma distribution. One attractive feature of the gamma distribution is that it is mathematically tractable. Computationally, frailties can be viewed as an unobserved covariate, estimated through the Estimation-Maximization algorithm. An alternative approach, less computationally demanding is to treat the frailty model as a penalized Cox model. This algorithm is available in the public domain statistical package \R, which was used for the analysis. 



\section{Data} \label{sec:data}

The data considered throughout this paper was collected in a double-blind, placebo-controlled, community trial in northeastern Brazil in a cohort of children aged 6--48 months at baseline and followed up from December, 1990 to December, 1991. Children were randomly assigned to receive vitamin A or placebo at the start of the trial and every four months thereafter. Demographic, health, nutritional and socioeconomic characteristics were similar in the treatment and the placebo groups. Children were visited three times per week at home by local field workers who collected various morbidity data including the occurrence of diarrhoea. Full details of the clinical trial have been published elsewhere \cite{barreto116} including discussion of a number of children who were lost to follow-up (mostly through migration away from the study area) these being roughly equally distributed between the two groups. Subsequent analyses in this paper focus on the resulting set of 1114 children (562 placebo group, 552 Vitamin A supplement group).

\section{Modelling strategies} \label{sec:mod}
Diarrhoea was measured through the number of days with diarrhoea in each week, such a day being defined according to the standard  definition ($\geq 3$ liquid or semi-liquid stools in 24~hr). The longitudinal models assumed that the days with diarrhoea, named $Y_{ij}$, followed a Binomial distribution:

\begin{equation}
Y_{ij} \sim Bin(\pi_{ij} , m_{ij})
\end{equation}
where $\pi_{ij}$ is the daily probability of diarrhoea of child $i$ and week $j$ and $m_{ij}$ is the numbers of days in week $j$ that the child $i$ was followed up. The logit link function was adopted to associate the covariates and the daily probability of diarrhoea $Pr(Y_{ij})=\pi_{ij}$. The measure of association derived from this models is the odds ratio (OR). 
%ATENCAO MUDANCA ACIMA!!!!

We decided to apply both an auto-regressive (AR) and an exchangeable (EX) correlation structure. The introduction of an AR correlation structure indicates that the probability of a day with diarrhoea decreases with a time lag. Longer episodes and short time between episodes, or very short episodes (or no episode), folowed by long healthy periods would be well fitted with this supposition. On the other hand, an EX correlation structure means that the relationship between the probability of any two days with diarrhoea is a fixed parameter. 

% CONFERE O TEXTO????

% descrever estrutura dos random effects
% descrever burnin, sample, etc\\



For the survival model each episode of diarrhoea was defined as a sequence of days with diarrhoea and the episode was considered finished after three days without diarrhoea, to allow for intermitent episodes. During a diarrhoea episode the child was excluded from the group at risk. Figure \ref{fig:sobrev} shows the time line for five children. Each vertical stripe marks the beginning of a diarrhoea episode. Notice that after each one the time line is blank, as the child is not in risk of a new episode. Loss of follow-up is marked as a "\textbf{o}" (Child 5). The data set is organised as follows: variable \textit{reg} is the unique identifier of each child; \textit{beg} is the begining of each observation period without diarrhoea while \textit{end} indicates either the begining of a diarrhoea episode or a period not observed; \textit{sex} and \textit{age} are examples of covarites and \textit{group} is either \textit{pla}, for placebo, or \textit{vit} for vitamin A; \textit{status} is $1$ for the begining of a diarrhoea episode and $0$ for censoring. Observing, for example, child 2 (feminin, 48 months old) record: she was followed-up since the first day, and was counted in the population at risk until day $244$, when her first episode began, lasting for 3 days. She re-entered the risk group at day 247, and this period is registred in line 2 of the dataset. After three more episodes she was followed until the end of the study time, when she was censored. 

% colocar a figura ao lado do quadro


\begin{figure}
\begin{center}
\includegraphics[angle=-90, width=6cm]{sobrev.eps}
\caption{\small Time lines of children indicating at-risk group and events to fit an AG model.}
\end{center}
\end{figure}

\begin{Verbatim}
   reg group  sex   age beg end status
     2   pla  fem    48   1 244      1
     2   pla  fem    48 247 286      1
     2   pla  fem    48 289 291      1
     2   pla  fem    48 293 361      1
     2   pla  fem    48 363 379      0
\end{Verbatim} 





The modelling strategy was to keep all covariates that were significant in at least one of the models in order to compare the models. We have quality of the drinking water, group (vitamin A or placebo), sex of the child, age in months as a time dependent covariate, type of water supply, number of domestic appliances (radio, refrigerator, television), duration (in months) of exclusively breastfeeding, mother's level of education, as covariates.

% R, lib survival, geepack e bayesX (indicar o script no final)






\section{Results} \label{sec:results}


\section{Discussion} \label{sec:discussion}



The first part of this formulation is similar to that of the standard GLM, but the important distinction is that it involves the \textit{conditional} mean and variance of the response given values for the random effects. This may be contrasted with the GLMs in Section~\ref{subsec:marginal} which model the \textit{marginal} mean and variance of the response i.e. that  averaged over the probability distribution for any unobserved random effects (subject heterogeneity). In the random effects model a within-subject covariance structure is indirectly generated by assuming that each individual has unique random effects, distributed around the fixed population effects\vcite{20}. In the marginal model subject-specific random effects are `averaged out' of the model, within-subject correlation being handled via a function specified separately to the mean structure. The important point here is that there is a subtle difference in the interpretation of fixed effect parameters in the marginal and random effect models. In the former, they are `population averaged' effects applicable to groups of subjects exhibiting similar heterogeneity to that observed amongst the subjects in the study. In the random effects model the fixed parameters are interpreted as effects on an
`average subject' in the study. Some have argued that the latter inference is more useful in clinical trials, where assessment of treatment effect on a average patient is of more interest than its average effect on a population as a whole\vcite{21}.


All this modelling exercise allows us to make a few assertions:
\begin{itemize}
\item Age is by far the most important covariate in diarrhoeal disease determination.
\item The best family socio-economic indicator is the number of household appliances.
\item The mother's education level is not as associated with diarrhoeal disease in this community as described in the literature, in spite of a significant number of mothers in all categories (281 with no education, 588 with primary school and 326 with high school). 
\item Calendar time, or time since the beginning of the study, when modelled, presented a strong protective effect. It should be noticed that it was included in the models in a linear way, and visual inspection of the trend showed a non linearity.
\end{itemize}

A significant effect of vitamin A could be only detected preventing the total number of liquid and semi-liquid motions, and not for days with the disease neither for days free from the disease. Considering the severity models, the random effects one presents one advantage over marginal model: it is able to capture the larger variability among children, via the random intercept. When marginal model is used the interpretation of the coefficients combines the within-subject relationship with a between-subjects relationship.

The model for days with diarrhoea presented different effects for our main interest, depending on the working correlation structure: significant effect when exchangeable structure were adopt, and not significant when an auto-regressive was chosen. Unfortunately there is not a specific method to choose the correlation structure. When the responsible variable is continuous, the correlation structure of the observed data serves as a guide for its choice, not our case. Considering the literature our results are moving in the opposite direction because the GEE analysis was not robust against the wrong choice of a correlation structure for non-continuous data. So the choice of the correlation must be made based on the knowledge of the researcher on problem to be analized.

The most interesting feature of the survival model with random effects is its ability to capture and measure through the R$^2$, the unexplained variability, in spite of its limit in including the calendar time. It should be noted that 
usually at this point the question that arises is: which is the best model? How to choose among them? First of all we should discuss in theory which one is a better approach to the problem. Besides, many improvements could be made to those models using the same data set, for instance, the inclusion of the time in a non linear way, as well as to the study design, in order to treat the possible learning effect detected via the decay time trend. What we want to emphasise is that \textit{All models are wrong. Some models are useful} \cite{box79}. And what we really want to know is: should we propose the supplementation of vitamin A as an strategy for controlling diarrhoeal disease in deprived populations? Well, the answer is not as straightforward as we would like. Although we have some evidence in favour of a real protective effect of vitamin A on the severity of diarrhoea, even when controlled by socio-economic level, mother's education and so on, the time trend can be a strong confounder, and should be properly addressed. Besides, obviously a cost-effectiveness study should be done, comparing the efficacy of this strategy with alternatives such education, casually detected via the time trend. Statistical modelling is a tool in medical sciences and an art. No model itself is right, and we should consider previous knowledge, robustness, performance using different suppositions, instead of just looking at p-values and positive results.

\bibliography{tese}
\bibliographystyle{unsrt}


\end{document}
