{'fraction_non_alphanumeric': 0.04694556429405415, 'fraction_numerical': 0.014715925362941423, 'mean_word_length': 4.495641646489104, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamic Modelling of Health and its application to the large scale analysis of Body Mass Index, using data from consecutive set of surveys

The methods used so far for the analysis of time changes in population health suffer from the lack of causality in their design. This results in problems with their implementation and interpretation. Here the method is presented with causality directly implemented in the design. This is done by, first, building up a dynamic model of population, postulating existence of Driving Force acting at subjects, while they move along their cohort lines, causing the changes of their substantial health indicators , State Variables, at rate proportional to this Force. The correspondent rates , named Cohort Trends, or C-trends, describe health history in each birth cohort. Having initial value and C-trends , the model allows to calculate health level (the means of State Variables) in each birth cohort, and thus, in the whole population. The task for statistical method is to identify the dynamic model (evaluate C-trends and Initial values) using data from a set of consecutive independent cross-sectional surveys. This is done by an iterative algorithm, running multiple regression procedure at each step, until the specified smoothing conditions are fulfilled. The algorithm can operate with surveys having different age ranges. The illustrative example shows the results of analysis of Body Mass Index for men , using 7 surveys in period 1972-2002 with age ranges 25-64 and 25-74. Since C-Trend is proxy for Driving Force, the year - age pattern of C-Trends provides unbiased information for health authorities on efficiency of health promotion actions or negative effects of uncontrolled harmful factors.


Introduction
This paper deals with the principles and methods of the analysis of time changes in population health, based on dynamic modeling of population.
The first step of analysis is building up a dynamic model of population, postulating existence of Driving Force acting at subjects, while they move along their cohort lines, causing changes of their substantial health indicators, State Variables, at rate proportional to this Force. The correspondent rates, named Cohort Trends, or C-trends, describe health history in each birth cohort. Having initial values and C-trends for each birth cohort, the model allows to calculate health level (means of State Variables) in each birth cohort, and thus, in the whole population.
In dynamic view, all traditional health indicators fall into two main categories: Driving Force (this includes smoking, physical activities, diets habits) and State Variables (among those BMI, Blood pressure, Cholesterol levels). The third category includes all class indicators, like gender, social group, place of residence.
Driving Force and State Variable indicators play different role in dynamic model and should be analyzed using different schemes.
In this paper we focus on analysis of State variable BMI. The data for this variable are taken from a set of independent consecutive health surveys, with randomly selected samples of population. Age ranges of samples and time periods between surveys could be different. Total calendar period of analysis could be quite large, in the example presented here, it is 30 years, covering 7 health surveys.
The methods used so far, addressing time changes in population health, did not differentiate between Driving Force and State Variable indicators. The problem has been approached by assessing trends over time for means and other statistics (e.g. percentiles) of the agecategories-specific distributions of parameters of interest, such as traditional risk factor indicators (e.g. systolic and diastolic blood pressure, cholesterol, body mass index), their categorical derivatives (e.g. prevalence of high blood pressure, prevalence of high cholesterol, prevalence of obesity) and event and mortality rates, see, for example, (Tunstall-Pedoe et al., 2003;Dobson et al., 1998a;Dobson et al., 1998b;Kuulasmaa et al., 2000;Gregg et al., 2005;Kautiainen et al., 2002;Chen et al., 2003;DiLiberti and Lorenz, 2001).
Various terms are used in literature for such trends: "trends" (Tunstall-Pedoe et al., 2003;Dobson et al., 1998a;Dobson et al., 1998b;Kuulasmaa et al., 2000), "secular trends" (Gregg et al., 2005;Kautiainen et al., 2002), "time trends" (Holford, 1991). Here we will call them "secular trends". There were several modifications in methods used to assess secular trends. One of these, "trends by linear regression", was the key element of the analysis in the WHO MONICA Project (Tunstall-Pedoe et al., 2003;Dobson et al., 1998a;Kuulasmaa et al., 2000). (Tunstall-Pedoe et al., 2003) contains extended list of publications, the survey-specific tabulations could be found in (Tolonen et al., 2000). First, the age-group specific trends were assessed, using linear regression; then they were aggregated using direct age standardization with fixed weights. The aggregated trends were subject to correlation analysis in order to test MONICA hypotheses.
A different modification uses the multiple logistic regression procedure applied to the whole set of data (Gregg et al., 2005). As a result, the marginal characteristics are obtained directly from the procedure. This is equivalent to direct standardization, with weights corresponding to the analyzed population.
In previous two examples, the method was applied to the data of wide age-range (40 years and more), while spans between consecutive surveys were 3-10 years. Some studies of adolescents deal with samples of age range 5-8 years, being sampled every year or every other year. In that case the trends are first examined visually by age (Kautiainen et al., 2002;Chen et al., 2003), since, as a rule, they exhibit a wide diversity of age-specific patterns.
To cope with such a diversity, several approaches have been commonly used. The first one is the APC ( Age, Period, Cohort) adjustment (Selvin, 1996;Holford, 1991), which accounts for high variability of secular trends across age, time period and cohorts (Chen et al., 2003) The second one divides the overall time period into several segments (DiLiberti and Lorenz, 2001), for which the corresponding plots suggest linear trends. The third one first aggregates the age-specific values, using direct age standardization, then estimates trends by applying linear regression to the aggregated values (Kautiainen et al., 2002).
The above examples bring to the fore the following fact: in many cases the plots of the analyzed parameters versus time display wide diversity across age-groups, and do not suggest linear trends within age groups. Even if the analysis is applied to the data comprising only two surveys (Dobson et al., 1998a;Dobson et al., 1998b), the unmeasured data between these two surveys may exhibit high nonlinearity.
Various types of summary rates, have been criticized, for example, by Holford (Holford, 1991): "Although this approach has the advantage of simplicity, it suffers from significant limitations; important details in the trends are lost in the averaging process involved in generating a summary rate. In many instances, these details have contributed significantly to the understanding of time trends for disease". However, the APC adjustment suggested by him (Holford, 1991) and Selvin (Selvin, 1996), as an alternative to summary rates, is also not free of problems. The main one stems from the fact that parameters used for adjustment "are hopelessly entangled" (Holford, 1991).
A variety of methods have been suggested recently to solve the problem of identifiability in the APC analysis. One of them, Partial Least Square (PLS), is illustrated in (Jiang et al., 2013) in analysis of BMI. (Havulinna, 2014) claims that model identifiability is becoming less of a problem with Bayesian APC models, and suggests two extension. one involving interactions between the age, period and cohort effects. Second extension uses autoregressive integrated (ARI) models.
At the same time, in review of the latest methods of Age-Period-Cohort Analysis, it was acknowledged that "Although a variety of approaches has been proposed to solve the APC conundrum, each has limitations. Yet another challenge is a criticism often lodged against general -purpose methods of APC analysis, namely, they provide no venue for testing specific, substantive, and mechanism-based hypotheses and thus are mere accounting devices of algebraic convenience that may be misleading." (Yang and Land, 2013, p3).
According to (Holford, 1991), secular trends "are significant because they can be highly suggestive as to what might be expected in the future and they are an effective approach to understanding disease etiology". We believe, that in many cases, the real data reject the idea of such a trends. For example, the real data on population size shown in Figure 1 hardly suggest age-specific linear trends over time. Rather they show fairly smooth changing along cohort lines.
To develop dynamic approach, a certain knowledge and skills are required not only in statistics, but also in physics, theory of control, theory of optimization, differential equations and others (Luenberger, 1979). Perhaps, this is the reason why the dynamic systems approach so far has not been widely used in population health research, though some examples of such applications can be found in the literature (see, for example, (Hargrove, 1998)).
General purpose of the dynamic system applications is formulated in (Brown University, 1995): "These three categories correspond roughly to the need to predict, explain, and understand physical phenomena" (for more discussion see also (Dobbin and Gatowski, 1999).
In turn, to build up a dynamic model, first we derive some principles, which we call the Principles of Dynamic Modeling in Health Research. These Principles are independent of the target task, so they could be applied to any other task, for example, to follow-up analysis with end-points.
The aims of this paper are as follows: • to develop the Dynamic Model of Population Health for the case of one State Variable • to build up the Dynamic Regression Method and algorithm allowing for multiple surveys with different age ranges with large volume of data involved • to run the illustrative analysis Note that each of the parts above is worth of more detailed, separate presentation.
Therefore, the challenge was to provide concise and logically completed descriptions of all parts, clearly outlining logical interrelations between them.
The earlier version of the Dynamic Regression Method was developed and presented by (Moltchanov and Mik'halskii, 2008), where C -trends were suggested as an alternative to circular trends, used so far. Historically, the method developed in MONICA for checking consistency of the reported demographic data ( (Moltchanov et al., 1999)) served as the prototype for the method developed by (Moltchanov and Mik'halskii, 2008). Some general aspects, such as criteria for commonly used health indicators to serve as system State Variables, have been considered earlier by (Moltchanov, 1993).
Section 2 describes Dynamic Model paradigm being applied to health research on Individual and Population Level. In Section 3 we present the analytical form and general form of statistical model for the case of continuous, normally distributed one parameter. The data used for analysis are described in section 4. The detailed description of the DRM "fast" algorithm is given in section 5. The example of application, is given in section 6. Section 7 contains conclusion and discussion.

Dynamic model on individual level, notation, definitions
To describe population history, it is convenient to use y-a plane, where y is real-valued calendar time in years, vertical axis, a is real-valued age in years, horizontal axis, with points specified as (y, a) (vertical coordinate first). This is matrix standard for indexing elements, which will be of use later.
Consider a population defined on observational frame (real compact) C: The usual assumptions must be fulfilled for the population: it should include all the residents of a geographically outlined area with well defined administrative boundaries. And neither administrative boundaries, nor rules for residentship are changed within observational frame C.
To introduce dynamic model terminology, we consider health history of one subject of the population. Let x W (y, a) be a weight of this subject at point (y, a). After time period ∆t, measured in years, the subject will move to the point (y + ∆t, a + ∆t), located on birth cohort line for this subject. For the sake of simplicity, we assume that weight changes linearly over this time period and both points belong to the observational frame C. Under these assumptions, the rate of change of weight at point (y, a) will be expressed as According to the current state of knowledge, weight change in adult is, in fact, change in amount of body fat, which is determined by balance of calories taken with meal and burned throughout the body activity over a certain time period (see, for example, http://www.weightlossforall.com/calories-per "One pound of body fat equals roughly 3500 calories."). Or one kilogram of body fat equals 7716.2 calories. Using this ratio, we can calculate how balance of energy results in weight change.
Let function f B (y, a) be the balance of energy in kilocalories per time unit. Then the above verbal statement could be expressed as u W (y, a) = k conv · f B (y, a), where k conv = 0.1296 (3) or, introducing the weight-specific Driving Force f W (y, a), Note that alternatively we could define f W (y, a) = f B (y, a), and get the equation u W (y, a) = k conv · f W (y, a), which is the form for general case (any State Variable, hence subscript is omitted): u(y, a) = k · f (y, a), k is a fixed coefficient.
Let t be a relative time, t = 0 when subject is at point (y, a). In general case, to represent time change of a state variable, we use function v(t) -a smoothed version of x(y + t, a + t), being free of daily cycles, which are common property of biological indicators ( for details see (Moltchanov, 2012). Then (5) will be transformed into: Now we rephrase the above description of the case using terminology introduced in (Moltchanov, 2012) for dynamic modeling of population health: An Object ( our subject) is moving over time along cohort line carrying, as a system, its State Variable x W (y, a) (weight in our case) and being all the time affected by weight-specific Driving Force f W (y, a). The Law of Motion postulates that change of State Variable occurs due to Driving Force at rate proportional to the value of this force, which is reflected in (6).
For the rate of change of a State Variable u(y, a), while the object moves along cohort line, we will use the term "Cohort trend" or "C-trend", and, for linguistic convenience, we will use term "Modifier" interchangeably with term "Driving Force" as proposed in (Moltchanov, 2012).
Note, that Driving Force, as a function of time, may be a non-continuous function. While State Variable, according to (6), is a continuous function of time, moreover, it is rightdifferentiable one. That property labels health parameters, which may serve as State Variables. Along with weight, the measurements of blood pressure, cholesterol, hight, as well as schooling years, are State Variables. While current smoking status, physical activity, dietary habits, along with other behavioral characteristics, are not State Variables, rather they are Modifiers.
Observe that continuous function of any number of State Variables is itself a State Variable. The expression (6) remains valid for this variable with the resulting Modifier being a linear combination of the contributing Modifiers.
In our future consideration and example we will deal with such a variable, The Body Mass Index, defined in pseudo-code notation as Let x BMI be State Variable for BM I. H -height of the subject in meters, which we assume to be constant for a subject in analysis frame C. To get equation for Law of Motion for BM I, we have to divide both sides of (6) by H 2 : Here we have introduced notation f BaH (y, a) for the Energy Balance adjusted for height.

Dynamic Model on Population Level: Axiomatic Setup
We use the term "Population Level" rather than "Population", since, our target population, in fact, might be a gender-specific subpopulation. Each subject may enter this population due to birth (if a min = 0), or crossing left-low boundaries, or migration in. Each subject may leave this population due to death, migration out or crossing the right-upper boundaries. If a subject with coordinates (y 0 , a 0 ) is within the population during time t , at that time it has coordinates (y 0 + t, a 0 + t). Thus, we may say that it is moving along cohort line.
Consider all subjects having coordinates on half-open interval ((y 0 , a 0 −∆a), (y 0 , a 0 )] at time t = 0. At time ∆t all those left in population will arrive at ((y 0 + ∆t, a 0 − ∆a + ∆t), (y 0 + ∆t − ∆a, a 0 + ∆t]. In other words, the birth cohort of width ∆a moves from (y 0 , a 0 ) to (y 0 + ∆t, a 0 + ∆t). We may think of such a cohort as of a container moving on plane (y, a). The contents of each container in process of movement is changed due to migration and death. If the rate of contents update is negligible (say, less than 1% per year), we may ignore it in our analysis. If not, the analysis has to take this into account.
Each container fits the definition of the dynamic model object, if we regard the corresponding State Variable as mean of State Variables for currently available subjects. The dynamic equation then could be obtained from ones for each subject, having form (6), by taking means of both sides: Since the whole selected observational frame could be covered by collection of non-overlapping cohorts of selected width, we may conclude that, in case of population, the overall dynamic model is a collection of the dynamic models specific for each cohort.
In fact, it is more convenient and straightforward to go ahead with the theoretical abstraction for birth cohort, which deals with cohort containers of infinitesimal age range, characterized by multidimensional distribution of the parameters of interest, not by physical subjects.
In general case, we suggest that there potentially exists a set of random variables (r.v.) X i , i = 1...k representing the corresponding set of measurable indicators of interest (State Variables) defined at each point (y, a) of compact C.
In this paper we restrict ourselves to the case of one indicator, so that subscript of X will be omitted. To make the following description more illustrative let us keep in mind the Body Mass Index (BMI) as an example of the indicator in question. We introduce the following notation: v(y, a) . = E(X(y, a)).
For the sake of simplicity while describing the core dynamic model, we assume: The dynamic equations describe changes of the distribution of r.v. X for a birth cohort taken at point (y, a) over time interval dt: We introduce function u(y, a), rate of change of the State Variable along cohort line: The Low of Motion postulates existing at each point (y, a) such a Driving Force f (y, a), which causes change of the State Variable, while point (container) moves along cohort line, at rate proportional to the value of the force, u(y, a) ∝ f (y, a), or, with f (y, a) properly scaled, Note that the above formulation of Low of Motion and equation (11) are similar to those for model on individual level (equation (5) ). The principal difference is that we are dealing now with population means v(y, a) ) and its cohort trends u(y, a), and population means of Modifier f (y, a).
In this paper we consider the case, when Driving Force f (y, a) does not depend on the properties of the cohort (value of v(y, a) ). The generalization of the model for the case of multidimensional distribution and state-dependent dynamics is formulated in (Moltchanov, 2012) For the sake of convenience we will use terms "Mean levels" or "levels" for the values of function v(y, a).
Let v 0 (y, a) be the value of v(y, a) at low-left boundary of the compact C for a (birth) cohort crossing the point (y, a): Then v(y, a) can be expressed as Thus, if the values of v 0 (y, a) at low-left boundary and u(y, a) on C are known, then the function v(y, a) could be evaluated for each point on C. In other words, we can set up initial levels v 0 (y, a) and the Driving Force pattern on C and we can simulate behavior of the population health in terms of mean levels v(y, a), using (11), (12), (13). Observe At the same time, continuity on y and a is not an obligatory property of this function.
3 Statistical Model: The core formulation

General Formulation of the Task
Suppose that a set of measurements is available (x k , y k , a k ), k = 1, . . . , K, for subjects of gender-specific subpopulation (men, for certainty), selected in a set of the independent crosssectional surveys. We assume that for each survey, the stratified by gender and age group random sample scheme was used. The age group stratification could be different in different surveys, as well as age range.
To start with analysis, first, we have to define the observational frame for the analysis, by setting up parameters in (1): The general formulation of the task is to estimate the functions v 0 (y, a) and u(y, a) on C, using the available measurements (x k , y k , a k ), k = 1, . . . , K.
To solve this problem, one option would be to formulate the optimization problem in functional space: to minimize the functional I: applying some additional requirements on functions u(., .) and v 0 (., .), such as continuity (piece-wise continuity), and/or restricted variation.
However, it seems more convenient to transform the above problem into the discrete -scale analogue and to take the advantage of the simplicity of the analysis and adaptation of the numerical methods available in the standard statistical packages.

Discrete-Scale Model
Let i and j be an integer value of time in years and an integer value of age in years correspondingly. Our intention is to build up the integer-values proxies of the equations (9 -16).
Let P (i, j) be a parallelogram-shaped element (convex hull) defined by its angle points: excluding its left and upper boundaries, which could be written as We impose for function u(., .) the conditions of being constant on each P (i, j) and for functions v(., .) -being constant on a and linear on y with constant slope u(i, j).
Formally this could be expressed as follows: u(y, a) = u(i, j), ∀i, j, y, a : (y, a) ∈ P (i, j) v(y, a) = u(i, j) · (y − i) + v(i, j), ∀i, j, y, a : (y, a) ∈ P (i, j) We derive minimal and maximal values for i and j from the correspondent values for y and a using definition (17): (i min , j min ) : (y min , a min ) ∈ P (i min , j min ), (i max , j max ) : (y max , a max ) ∈ P (i max , j max ) For convenience, from now on we will use relative scale for age and time, defined by trans- Consider functions u(i, j) and v(i, j) defined on integer-valued two-dimensional domains Now the main dynamic equation (10 ) could be rewritten as Let v 0 (i, j) be the value of v(., .) at low-left boundary of the domain V corresponding to a (birth) cohort crossing the point (i, j): Combining (22) and (23), we rewrite equation (12) as: From (24), it follows that if v(i, j) is set up on the low-left boundary of V and u(i, j) is set up on the whole U then v(i, j) could be calculated for the whole V.
Finally, assembling (9), 24) and (19) for each available observation (x k , y k , a k ), k = 1, . . . , K, we obtain: Let z be a vector with components v 0 (i, j) and u(i, j) ordered in the following way: Using vector z and introducing vector of coefficients b k , we can rewrite (25) in the form This form represents a particular case of Gauss-Markov Setup for the Least Squares Linear Estimation problem (Rao, 1973).
Let B 0 be a matrix composed of row vectors b T k in (27), z and x 0 stand for column vectors of the parameters z j and the variables x k correspondingly, and S 0 be a scalar function defined as Note that if rank (B 0 ) = dim(z), then estimates obtained by unconditional minimizing of function S 0 (z) are unique ones. Such a case takes place only if the observations cover all the elements P (i, j) when surveys cover the whole analysis period without gaps.
In practical cases, minimizing of S 0 results in singular or ill-posed Inverse Problem, and so-called regularization techniques are needed to obtain meaningful solution estimates. Most of these techniques employ the idea of smoothing of some function having clear physical interpretation (Neumaier, 1999).

Smoothing
We define the following indicator of smoothness of function u(., .) (C-trends); We define the following indicator of smoothness of function v(., .) Each term in this sum represents the square for a proxy of the second derivative of function v(., .) with respect to age or with respect to calendar time at point (i, j).
Replacing v(., .) by v 0 (., .) and u(., .) using (24), and the last ones by vector z, we will transform the previous expression to the following form: Now we can add one or both constraints S k (z) ≤ α k with some selected α k ≥ 0, k = 1, 2, to the model (27). Observe that indicators S 0 , S 1 , S 2 are quadratic functions in finite vector space E n with elements (vectors) z and n = dim(z). The optimization problem for point estimation for our case, could be formulated as min x∈En S 0 (x), subject to S k (x) ≤ α k , with given α k > 0, k = 1, 2.
Let n 0 , n 1 and n 2 be numbers of rows in matrices B 0 , B 1 and B 2 correspondingly. Let λ 1 , λ 2 be some non-negative scalars. Introducing matrices and vectors where I 0 , I 1 and I 2 are identity matrices of rank n 0 , n 1 and n 2 correspondingly, we can formulate the problem of least squares estimation in the following form (a modification of Gauss-Markov setup which fits form of Aitken setup (Rao, 1973) x for which the point estimation problem is In (Moltchanov and Mik'halskii, 2008) it was shown that problems (31) and (34) are equivalent: problem (31) with given α 1 , α 2 possesses the same solution as problem (34) with some λ 1 , λ 2 , and vice versa, or both don't possess any solution. In particular, it was shown that for existence of a unique solution to problem (34) it is sufficient to have 4 data points such that the corresponding points (y, a) on plane y, a satisfy condition: no any 3 of them are located on a common straight line.
More discussion on statistical properties of problem (34) could be found in (Moltchanov, 2012). As soon as parameters λ 1 , λ 2 are given in setup (32, 33), the following could be obtained routinely:ẑ -point estimate of vector z, covariance matrix of this estimate Cov(ẑ), andσ 2estimate of σ 2 .
The task remains, to formulate criteria of smoothness of a solution above and to organize an iterative process by selecting λ 1 , λ 2 and solving Linear Regression Problem (33) at each iteration, until the predefined criteria are satisfied.
The design of the corresponding algorithm should reflect substantially the size and the structure of data to be analyzed.

Data
To illustrate the method and to demonstrate its performance, the data will be used comprising 7 cross-sectional surveys, conducted in North Karelia, Finland, during the period 1972 -2002. Details of these data are shown in Table 1: • Study population: North Karelia, Finland, men.
• Source of data: cross sectional independent surveys conducted in years 1972 -2002 every 5 years.
• in year 1972 survey was conducted in the period of 8 months, February -September, in year 1997 -period was 6 months, January-June, in other years surveys were conducted in 4 or 3 month, starting in January.
• Participation rate varies from 66% to 94% • Total number of observation with non-missing BMI, age and date of examination, is 11045.
Original measurements of interest were gender, date of birth, date of examination, weight and height. At the study sites, height and weight were measured using a standardized protocol. Height was measured to the nearest 0.1 cm. Body weight of the participants wearing usual light indoor clothing without shoes was measured with a 0.1 kg precision on a balanced beam scale.
The analysis variables included in the model are: BMI -the Body Mass Index, defined as weight(kg)/height(m) 2 . AGE -age in full years, defined as year of examination minus year of birth.
YEAR -date of examination measured in years.

Data aggregation
The original individual data might be of quite large size (11045 in our example), which equals to rows number n 0 of matrix B 0 in (32), and, along with columns number, proportional to product of year range by age range, this results in resource-expensive computation at each step of iteration. Data aggregation reduces significantly number of rows n 0 . Aggregation is applied to the original measurements (x k , y k , a k ), k = 1, . . . , K, producing summary statistics for (age · year) cells with size 1, containing arithmetic means (x c , y c , a c ), and number of original measurements in each cell (n c ), where c = 1, . . . , C -collection of non-empty cells. In general case, the cells are excluded with n ≤ n exc . In our example, n exc is set to 5, and there were no excluded cells, and number of cells is equal to sum of all age ranges.
For convenience, we rename the aggregated data to the form of the original ones, (x k , y k , a k ), k = 1, . . . , K, where K = C. After that all the above formulas remain valid.

Analysis Domain
Each data point at cell (i, j) causes inclusion into the analysis u(i, j) from the cells (i−k, j −k), k = 0, ...min(i, j). Overlapping of such cells creates collection of cohort segments starting at some points on low-left border of domain U. Some of such segments contain several data points, other only one.
For future use, we define selection options: • Domain=1: all cohort segments are selected; • Domain=2: only cohort segments with 2 and more data points are selected.
At the next step, selected collection is modified by inclusion additionally those cells which fill the internal gaps in all vertical and horizontal segments. The resulting analysis domain U a , may contain significantly less cells, compared with original U, thus, diminishing number of columns in the analysis problem. We introduce vector u ind , indicating (with 1/0) the components of vector u, corresponding to domain U a and vector u a part of vector u with components included in the analyses.
Domain U a also determines subset of components of vector v 0 included in the analysis, which is a segment described by i l and i r indexes of the first and last included components.

Between-cohorts smoothness instead of smoothness of v(., .).
We replace indicator of smoothness (29) by the following one: allowing form (30) with the corresponding matrix B 2 .

Iteration step analysis outcomes and control of iterations.
Let z ind be index vector indicating (with 1) the components of the original vector z (29) included in the analysis. Let, further, z a be subset of the components of z -vector of parameters to be estimated at each step of iteration. Correspondingly, all design matrices will be modified B i → B a,i by excluding the columns for components, not participating in the analysis. Besides, in matrix B 1 all the rows will be excluded, where not all the components belong to domain U a After point estimatesẑ a , covariance matrix of this estimate Cov(ẑ a ), and estimate of σ 2 will be found, we recalculate this solution to the original scale, gettingẑ, Cov(ẑ) and also Corr(ẑ) -Pearson correlation matrix. For this purpose we will use matrix A za2z : z = A za2z z a , in which we set all rows for non-participating components to missing values(.). So, that matrices Cov(ẑ) and Corr(ẑ) have non missing values only if both components are "participating".
Let r i,j,j+1 be coefficient of correlation of u i,j and u i,j+1 -horizontal link, and r i,j,i+1 be coefficient of correlation of u i,j and u i+1,j -vertical link. Observe that the closer to 1 the coefficients of correlations are, the more "smooth" is the plot of values of u i+,j horizontally or vertically.
We define the following measure of smoothness of surface u(i, j): where sum is taken only over non-missing values, and n ru is number of such a values. Thus, r u is an average coefficient of correlation between two adjacent horizontal or vertical u(i, j) valus in domain U a .
Similarly, we define measure of smoothness for initial values of cohorts The target condition for analysis is set up in terms of reference values r u and r v and accuracy levels δ u , δ v . At the end of each iteration the following condition is checked: If this condition is fulfilled, then iterations stop. Otherwise, the new values for weights are defined as follows: and the next iteration is started.
To measure difference in C-trends over age and calendar year, the pairwise comparison tests are performed for mean values of C-trends, evaluated for a set of age-year clusters, defined by cluster sizes, ∆ a and ∆ y .
Let U c be matrix of such mean values, u c =Shape(U c , 1) T and matrix A u2uc such that u c = A u2ucû . As soon as, matrix A u2uc is created for given ∆ a , ∆ y , U c could be calculated, as well as variance/covariance values for its elements in a format of C = Cov(û c ) = A u2uc Cov(û)A u2uc T . For each cluster, statistics and corresponding probabilities are computed for pairwise comparison of mean C-trends for current cluster and for adjacent one for older age group, and for current one and for adjacent one for the next calendar years period ( if the corresponding clusters exist). Using classical paradigm, this is done by testing linear hypotheses in form H 0 : u ci − u cj = 0. General expression for F-value (see for, example, SAS/Stat manual, (SAS Institute Inc., 2011c)) in this case takes a simple form Corresponding probability is computed using SAS function probF (see (SAS Institute Inc., 2012b)) as P r = 1 − probF (F, 1, n − r) Note, that in Bayesian view, these probabilities should be referred to as tail-area probabilities for posterior predictive distributions ( (Gelman et al., 1995) , p.169).
Results of pairwise tests are presented graphically in figure, produced by PROC GCON-TOUR, properly annotated ( see Figure 5 in example of application).
For reference, we will call this algorithm DRM3(A), with prefix DRM3 to differentiate it from those developed in (Moltchanov, 2012): DRM2(R) for "oRiginal" data, and DRM2(A)for "Aggregated" data.
6 Example of Application

Analysis Setup
The algorithm modification DRM3(A) is used, preprocessing original data, described in Table  1, into aggregated format. These data are shown in Figure 2. The analysis was set up for the (maximal) age range 25-74 and for the calendar year period 1972-2002.
To control iterations, the accuracy levels were selected δ u = 0.05, δ v = 0.05. We have run analysis for several pairs of reference levels r v and r u . Let R(r v , r u ) be set of outcomes (in terms of estimates, figures and tables) for analysis run with r v , r u ). As a main outcome, we present analysis R(0.7, 0.9). For comparison, we have produced also analyses R(0.7, 0.7), R(0.7, 0.8) and R(0.7, 0.95).

Outputs
The results of the each analysis are visualized the set of 3-dimensional figures and special plots.
First, we present results for analysis R(0.7, 0.9). Figure 2 displays the values representing means of BMI calculated for each age and year, for which the survey data are available (number of cases in each cell exceeds 9). To visualize the along-cohort changes, the columns corresponding to the same birth cohorts in different surveys are drawn using similar shades of gray.
Note that this figure is the same for different analyses R(r v , r u ). Figure 3 displays estimates for the mean levels of BMI for the whole domain, with study age range plus one year, and study period plus one year. Figure 4 displays C-trends with upper or low part of the 95% confidence intervals, shown at boundaries only. Figure 5 displays mean levels of C-trends for specified 5-year age-year clusters, with P-values for differences between clusters. Figure 6 displays mean C-trend with 95% Confidence Interval for each cluster. For each points, horizontal position indicates age range of the correspondent cluster. Figure 7 displays original data (mean ± CI), estimates of v and estimates of C-trend with 95% Confidence Interval along cohort line for a selected cohort.
These figures illustrate the principle "one figure is better than one hundred tables", though all the underlying data are available and could be presented in a set of tables. Figure 8, along with Figures 5 and 6, illustrates the choice of feasible solution.
1. Confidence limits for C-trends estimates are decreasing ( Figure 6), as well as maximal curvature of two-dimensional plot of C-trends estimates over age and year.
2. The curvature of C-trends plot along cohort is decreasing ( Figure 6 ). Relatively high curvature may occur due to a gap in the data (5 year gap between surveys) or due to odd values of some original measurements. In that case this plot exhibits non-monotonic function of age.
Thus, locally increased curvature is a side effect of low r u ; a feasible solution with reasonable high (r u = .9) is free of such side affects nad less sensitive to possible biased original measurements. At the same time, the odd original values might be not a bias in measurements, rather it could result from hight rates of migration-in or -out of population.
Summing up, we may conclude that, as a feasible solution for the data available, the outcome R(0.7, 0.9) could well be chosen. This is referred to as results in our further discussion. Figure 3 shows that mean BMI levels increase along cohort lines throughout the study period, although they are different for different birth cohort. Specific peaks and troughs follow cohort lines.

BMI dynamics: the main findings
Figures 4 6 show that C-trends are decreased with age; they differ significantly with periods, especially for the age groups 25-29 and 30-34. In particular, the alarmingly high C-trends are observed for the period 1997-2002, compared with period 1992-1997 for age groups 25-29,30-34 (significant difference is indicated in Figure 5). In addition, the high rate of increasing C-trends over year is observed for period 1992-2002 for age range 25-30.
The young generation coming to the observation frame exhibits rate of increasing BMI as high as 0.4 -0.5 units per year. Which is about two times higher than in years 1972-1982. This is clear challenge to health management.
Recall that C-trends are proportional to the external Driving Force (Modifier) which, in case of BMI, is the average hight-adjusted balance of calories. Therefore the above findings induce the task for expert in social and economical areas to find out, why difference in calories consumed with food and burned through physical activity is so alarmingly increased over time.

Conclusion and Discussion
In this paper we have presented a novel formulation of the key principles of dynamic modeling in application to health research, which justify the structure and interpretation of the core models dealing with C-trends.
In particular, according to these principles, traditional risk factors' indicators fall into two categories, State Variables and Modifiers (see section 2 ), having different dynamical nature and, hence, playing different roles in the model and analysis.
As corollary of this, circular trends for State Variables have no sense at all. At the same time, only State Variables may determine instantaneous hazard rate of failure. In dynamic models, causality is postulated: changes are due to Driving Forces (Modifiers), existing in the real world. In case of consecutive survey data, C-trends are believed to be proxies for Driving Forces, providing the tool for three main practical tasks: analysis, prediction and control of health on population level ( see section 2) We have used these principles as a framework for developing the dynamic model of simulating the temporal changes in characteristics of a real-world object -population. In the course of this process, first, we have identified two interacting objects, population and its environment, on the top aggregation level. Further system analysis has led us to breaking down the study population into a set of potentially infinitesimally narrow birth cohorts, carrying over time health state profiles expressed in terms of health related indicators (State Variables).
The model employs the health field concept, suggesting existence of an influencing factors (Modifiers), generated by environment and acting on the population, specific for each calendar year and age, and causing within-cohort changes of the health indicator with rate of change corresponding to the strength of this factors.
For illustrative purposes we have selected one-parameter case with continuous, normally distributed parameter and with strength numerically equal to rate of change. While keeping model reasonably realistic, these simplifications help to highlight the key properties of the dynamic model of population health and method of its identification -the Dynamic Regression Method.
In the illustrative example, we have shown that the Dynamic Regression Method applied to aggregated data, DRM3(A), provides a sensible view on the BMI dynamics. It reveals clear difference between dynamics the levels of the parameter and its C-trends. From practical prospectives, it is C-trends, not levels, which primarily seem to be modifiable by preventive activities or involuntary changes affecting the population. It is worth noting that outcomes from the DRM3(A) analysis serve as data for the next-level analysis, involving additional information and aiming at finding reasonable explanation of the observed dynamics (diagnostic property of DRM). One of the important complementary component for such an analysis is dynamics of the population size (we have developed a modification of the DRM for that type of data, this is a subject for one of the next publication). If there is significant migration "in" or "out" of the study population, the observed effects could be entirely or partially due to the population instability (health selective effect). The outcomes from the DRM analysis could be used straightforwardly for prediction of the age-specific profile of the State Variable, say, for 5 year period, by applying the C-trends at the last year of the study period to the estimates of the parameter's levels at that year. Such a projection will not cover the cohorts, not included in the study age range at the last study year.
Recall that this method has been developed as an alternative to the secular trends and APC approach used so far. In this respect, it is worth noting that the model presented here is characterized by local cohort trends (C-trends), which have clear interpretation: changes in the State Variable of the same physical entity per time unit. If we will formally calculate a characteristics resembling age-specific secular trend, we will obtain a difference between two different physical entities (birth cohorts), caught occasionally at the moments of measurement. Hence, it may behave quite arbitrarily. In other words, in the view of the dynamic modeling approach, secular trends do not exist in nature. In one special case only, when all the age profiles of a State Variable are the same over calendar years (stationary case), formally calculated secular trends will be equal to zero at each age within the study age range. Only in that trivial case, secular trends possess both, predictive and diagnostic power. However, even in this case, secular trends are kind of statistical fallacy, due to missing causality. As to APC approach, the main methodological drawback of it is treating Age, Period and Cohort as linked algebraically, while Cohort brings a differential component to the problem: namely, derivative of the target variable along cohort line is a function of Age and Period. This immediately removes the notorious conundrum of the APC method -linear dependency of the participating factors.
There are certain restrictions in using the current version of DRM3 methods, imposed by the size of the problem, due to using matrix operations. Switching to Bayesian framework and employing Markov Chain Monte Carlo methods (Gelman et al., 1995) may solve these problems.
The simplified dynamic equation used in the current model could be modified, accounting for the fact that rate of change may depend also on the current level of the State Variable. The more comprehensive model needs to be developed, comprising multiple State Variables, and corresponding C-trends as a linear functions of current State Variables. Also, the model is to be developed processing measurements which are Modifiers, not Sate Variables, such as smoking and physical activity. In some sense, this is an inverse problem to one presented here.
Collection of such models could be a powerful practical tool for prediction of population health for about 5 year span.