Impact of rural Family Physician programme on maternal and child health indicators in Iran: an interrupted time series analysis
Following the Alma-Ata Declaration of 1978, the Islamic Republic of Iran (Iran) began to expand its primary healthcare (PHC) network. This publicly funded and extensive PHC network offered health houses staffed by trained health workers, known as behvarz in Persian, to rural residents. Behvarzes were the main service providers and were responsible for provision of antenatal and postnatal maternal care for uncomplicated pregnancies, as well as other PHC components. A health house was a peripheral unit of the PHC network that served about 1500 inhabitants in a rural area. Every four to seven health houses were linked to a health centre, led by one family physician, located in the closest biggest village. As such, the health centre was providing a service to approximately 7000–10 000 people. This health centre was responsible for the provision of a defined health package for a population in their associated catchment rural area, which included outpatient visits, medical laboratory tests, medication, midwifery and other services. Overall, this well-organised PHC network substantially contributed towards reduction in infant, maternal and newborn mortality rates during the last two decades of the 20th century. Since 2000, the promotion of maternal and child health (MCH) care through the provision of related services and recording of MCH-related indicators became a top priority for the Iranian PHC network. As of 2005, with the start of the second PHC-based reform known as Family Physician (FP) programme, the way each health centre (now referred to as FP team) functions within the PHC network has changed. The FP team structure was expanded to include additional staff, such as a general physician, midwives and a pharmaceutical technician. Now, one family doctor was settled in each village, so the size of the rural population serviced by the FP team was limited to 3000–4000 people. A new insurance plan (Rural Health Insurance) was also scaled up to cover the whole Iranian rural population. The previously limited duration of the service provided by the FP team (from 08:00 to 14:00) also changed to become a round-the-clock service, with triple the number of physicians available in the rural areas. In the new settings, the family physicians were now team coordinators and service providers. There was also a significant increase in the salary of general physicians involved in the FP programme. The FP programme also considerably improved accessibility and affordability of medical and health services, which among others included special antenatal and postnatal care, laboratory tests (including measurement of haemoglobin (Hb) levels to monitor iron supplementation plan) and routine ultrasounds to detect fetus abnormality. The FP teams also started providing maternity care for all women in rural areas. Like Iran, many countries (eg, Pakistan, India, Rwanda, China, African and Latin American countries) are making serious sector reforms to provide effective and high-quality MCH services. The main aim of these reforms is to decrease maternal, neonatal and infant mortalities and achieve Sustainable Development Goals. However, there is no extensive evidence as to whether such health reforms have led to improvements in health systems performance and indicators. Iranian health sector reforms are no exception to this. After a decade since the FP reform implementation in the rural areas, one could expect considerable changes to be observed in MCH indicators in Iran. However, any possible changes were not yet extensively studied, evaluated and/or reported. Therefore, the main goal of this study was to investigate the impact of Iran’s FP reform on selected MCH indicators, using the data from 1994 to 2013 for the northwest Iranian province East Azerbaijan. 

For the evaluation of population-level interventions based on population indicators rather than individual measures, an interrupted time series analysis (ITSA) is recommended. ITSA is also a desired method for the evaluation of health policies, particularly when their targeted indicators have been repeatedly measured in a defined population over a time period. For our analyses, we developed ITSA regression models following the steps introduced by Bernal et al , using Stata V.13.1 (StataCorp, College Station, Texas, USA). Step 1: study design and its appropriateness Population and data The data were collected in the rural area of East Azerbaijan, a northeast province of Iran. The province covers an area of approximately 47 830 km² and constitutes roughly 5% of the country’s population (2013). In 2013, 26% of the population resided in rural areas, a decrease from 35% in 1994. The overall population growth rate in the rural area decreased from 1.42% (in 1994) to 1.31% in 2013. The total fertility ratio decreased from 2.88 in 1994 to 2.11 in 2013. Virtually all pregnancies in the rural areas of Iran are covered by health houses and related information is recorded in households’ health files from mid-1980s. As a result, the longitudinal data are available on the care delivered during pregnancy and childbirth, and development until school age, for more than three decades and for all women. For our study, we used these rural health houses data at the individual level. The sample frame consisted of a list of all pregnancies in different health houses in a specific year in associated geographical area. Since health houses have almost 100% coverage for their reference population in rural areas, the sampling was done in a simple random selection way. Given the average number of annual live births (roughly 2000 in 1994–2013), we randomly selected 10% of all pregnancy cases in a given year, which were then retrospectively followed up until the end of the pregnancy. This resulted in approximately 200 cases a year and a total of 4096 cases of pregnancies from 1994 to 2013. Intervention and comparison The multicomponent intervention, called the FP programme, was considered an intervention. This intervention was composed of a package of medical, maternal and childcare services. It also aimed to make healthcare services affordable through a new insurance plan. The programme was launched in 2005. The intervention time in the present study, therefore, ranged from 2005 to 2013. The comparison time, in contrast, was the time prior to the implementation of the FP programme (ie, 1994–2004). Outcome The Donabedian model was used for categorisation of MCH indicators. This model emphasises the importance of structure , process and outcome in healthcare services assessments. The mother’s age, level of education, occupation and gravidity were considered as structure variables in this study. They were used as confounders and their effects were controlled in all statistical models. The number of antenatal care visits (ACVs), history and frequency of ultrasounds, laboratory tests and type of delivery were considered as process indicators in this study. Maternal mortality rate (MMR), neonatal mortality rate (NMR), history of abortion and/or stillbirths were also considered as outcome indicators. The values for MMR and NMR were obtained from the annual national reports of health indicators for the study population. MMR, NMR, abortion and stillbirth ratios indicators were defined following WHO guidelines. Other outcome indicators included birth weight (BW), which was recorded regardless of the gestational age, and the value of Hb in the first test during the 3-month pregnancy period. Step 2: the impact model ITSA regression models were used to assess the causal links between the programme as an intervention and the outcomes of interest. To conduct the ITSA, the following three variables were added into the model: a time variable (trend), showing the time elapsed from the beginning of the study time frame (1994) until the end of the study period (2013); a dummy variable (programme) representing the time before and after the intervention; and an interaction term which represented the joint effect of the time trend of indicators and the implementation of the intervention. Step 3: descriptive analysis Descriptive analyses were illustrated as means (±SD) for numeric variables and frequencies for categorical variables. The descriptive time trend line charts were used to present the process and outcome indicators over the time. The scatter plots were used to show the changes of process and outcome indicators over the study time points (years). Step 4: regression analysis We applied segmented regression models, where each segmented line was fitted through ordinary least squares. The following model was used to show the relationship between the variables: y t = a t + b 1 I t + b 2 X I t + b 4 e d u t + b 5 a g e t + b 6 j o b t + b 7 n o . p r e g n a n c y t + e t where Y t is the outcome variable in time t and I is the indicator of the elapsed time. Moreover, X t is a dummy variable that denotes the ‘ Before the programme ’ (coded as 0) and ‘ After the programme ’ (coded as 1) periods. The interaction variable ( XI t ) shows the changes in the slope of response variable time trend during the study period, by considering the intervention time point. a t is the baseline level of the outcome ( Y ) when all independent variables are zero, b 1 is the unique effect of a time unit increase in pre-intervention time, b 2 is the mean change in the outcome following the programme if time remains constant and b 3 is the joint effect of time and the FP programme that represents the slope change following the intervention. The model was adjusted for potential confounders including age, mother’s education level and occupation and gravidity. To account for autocorrelation, a random error e t was added to the model. Step 5: addressing methodological issues Since data for outcome and process were reported annually, seasonality was not indicated in the analysis. We assumed no other regular fluctuation or long-term trend over the study time period. A time lag of 1 year was used for mortality outcomes (MMR and NMR) to adjust for the time of actual change and the time when intervention was received by participants. Since we did not use a concurrent control group, unmeasured time-varying confounders, such as socioeconomic change during the study period, might distort the results. Nonetheless, we assumed no time-varying confounders other than those that we adjusted for in the model. We also assumed no overdispersion in the data. The Durbin-Watson transformed statistic was used as an indicator for exploring positive or negative first-order autocorrelation. Patient and public involvement No patients or members of the public were involved in this study. 

study design and its appropriateness Population and data The data were collected in the rural area of East Azerbaijan, a northeast province of Iran. The province covers an area of approximately 47 830 km² and constitutes roughly 5% of the country’s population (2013). In 2013, 26% of the population resided in rural areas, a decrease from 35% in 1994. The overall population growth rate in the rural area decreased from 1.42% (in 1994) to 1.31% in 2013. The total fertility ratio decreased from 2.88 in 1994 to 2.11 in 2013. Virtually all pregnancies in the rural areas of Iran are covered by health houses and related information is recorded in households’ health files from mid-1980s. As a result, the longitudinal data are available on the care delivered during pregnancy and childbirth, and development until school age, for more than three decades and for all women. For our study, we used these rural health houses data at the individual level. The sample frame consisted of a list of all pregnancies in different health houses in a specific year in associated geographical area. Since health houses have almost 100% coverage for their reference population in rural areas, the sampling was done in a simple random selection way. Given the average number of annual live births (roughly 2000 in 1994–2013), we randomly selected 10% of all pregnancy cases in a given year, which were then retrospectively followed up until the end of the pregnancy. This resulted in approximately 200 cases a year and a total of 4096 cases of pregnancies from 1994 to 2013. Intervention and comparison The multicomponent intervention, called the FP programme, was considered an intervention. This intervention was composed of a package of medical, maternal and childcare services. It also aimed to make healthcare services affordable through a new insurance plan. The programme was launched in 2005. The intervention time in the present study, therefore, ranged from 2005 to 2013. The comparison time, in contrast, was the time prior to the implementation of the FP programme (ie, 1994–2004). Outcome The Donabedian model was used for categorisation of MCH indicators. This model emphasises the importance of structure , process and outcome in healthcare services assessments. The mother’s age, level of education, occupation and gravidity were considered as structure variables in this study. They were used as confounders and their effects were controlled in all statistical models. The number of antenatal care visits (ACVs), history and frequency of ultrasounds, laboratory tests and type of delivery were considered as process indicators in this study. Maternal mortality rate (MMR), neonatal mortality rate (NMR), history of abortion and/or stillbirths were also considered as outcome indicators. The values for MMR and NMR were obtained from the annual national reports of health indicators for the study population. MMR, NMR, abortion and stillbirth ratios indicators were defined following WHO guidelines. Other outcome indicators included birth weight (BW), which was recorded regardless of the gestational age, and the value of Hb in the first test during the 3-month pregnancy period. 

The data were collected in the rural area of East Azerbaijan, a northeast province of Iran. The province covers an area of approximately 47 830 km² and constitutes roughly 5% of the country’s population (2013). In 2013, 26% of the population resided in rural areas, a decrease from 35% in 1994. The overall population growth rate in the rural area decreased from 1.42% (in 1994) to 1.31% in 2013. The total fertility ratio decreased from 2.88 in 1994 to 2.11 in 2013. Virtually all pregnancies in the rural areas of Iran are covered by health houses and related information is recorded in households’ health files from mid-1980s. As a result, the longitudinal data are available on the care delivered during pregnancy and childbirth, and development until school age, for more than three decades and for all women. For our study, we used these rural health houses data at the individual level. The sample frame consisted of a list of all pregnancies in different health houses in a specific year in associated geographical area. Since health houses have almost 100% coverage for their reference population in rural areas, the sampling was done in a simple random selection way. Given the average number of annual live births (roughly 2000 in 1994–2013), we randomly selected 10% of all pregnancy cases in a given year, which were then retrospectively followed up until the end of the pregnancy. This resulted in approximately 200 cases a year and a total of 4096 cases of pregnancies from 1994 to 2013. 

The multicomponent intervention, called the FP programme, was considered an intervention. This intervention was composed of a package of medical, maternal and childcare services. It also aimed to make healthcare services affordable through a new insurance plan. The programme was launched in 2005. The intervention time in the present study, therefore, ranged from 2005 to 2013. The comparison time, in contrast, was the time prior to the implementation of the FP programme (ie, 1994–2004). 

The Donabedian model was used for categorisation of MCH indicators. This model emphasises the importance of structure , process and outcome in healthcare services assessments. The mother’s age, level of education, occupation and gravidity were considered as structure variables in this study. They were used as confounders and their effects were controlled in all statistical models. The number of antenatal care visits (ACVs), history and frequency of ultrasounds, laboratory tests and type of delivery were considered as process indicators in this study. Maternal mortality rate (MMR), neonatal mortality rate (NMR), history of abortion and/or stillbirths were also considered as outcome indicators. The values for MMR and NMR were obtained from the annual national reports of health indicators for the study population. MMR, NMR, abortion and stillbirth ratios indicators were defined following WHO guidelines. Other outcome indicators included birth weight (BW), which was recorded regardless of the gestational age, and the value of Hb in the first test during the 3-month pregnancy period. 

ITSA regression models were used to assess the causal links between the programme as an intervention and the outcomes of interest. To conduct the ITSA, the following three variables were added into the model: a time variable (trend), showing the time elapsed from the beginning of the study time frame (1994) until the end of the study period (2013); a dummy variable (programme) representing the time before and after the intervention; and an interaction term which represented the joint effect of the time trend of indicators and the implementation of the intervention. 

Descriptive analyses were illustrated as means (±SD) for numeric variables and frequencies for categorical variables. The descriptive time trend line charts were used to present the process and outcome indicators over the time. The scatter plots were used to show the changes of process and outcome indicators over the study time points (years). 

We applied segmented regression models, where each segmented line was fitted through ordinary least squares. The following model was used to show the relationship between the variables: y t = a t + b 1 I t + b 2 X I t + b 4 e d u t + b 5 a g e t + b 6 j o b t + b 7 n o . p r e g n a n c y t + e t where Y t is the outcome variable in time t and I is the indicator of the elapsed time. Moreover, X t is a dummy variable that denotes the ‘ Before the programme ’ (coded as 0) and ‘ After the programme ’ (coded as 1) periods. The interaction variable ( XI t ) shows the changes in the slope of response variable time trend during the study period, by considering the intervention time point. a t is the baseline level of the outcome ( Y ) when all independent variables are zero, b 1 is the unique effect of a time unit increase in pre-intervention time, b 2 is the mean change in the outcome following the programme if time remains constant and b 3 is the joint effect of time and the FP programme that represents the slope change following the intervention. The model was adjusted for potential confounders including age, mother’s education level and occupation and gravidity. To account for autocorrelation, a random error e t was added to the model. 

Since data for outcome and process were reported annually, seasonality was not indicated in the analysis. We assumed no other regular fluctuation or long-term trend over the study time period. A time lag of 1 year was used for mortality outcomes (MMR and NMR) to adjust for the time of actual change and the time when intervention was received by participants. Since we did not use a concurrent control group, unmeasured time-varying confounders, such as socioeconomic change during the study period, might distort the results. Nonetheless, we assumed no time-varying confounders other than those that we adjusted for in the model. We also assumed no overdispersion in the data. The Durbin-Watson transformed statistic was used as an indicator for exploring positive or negative first-order autocorrelation. 

No patients or members of the public were involved in this study. 

Structural indicators of education level, occupation and the age of a pregnant women had changed from 39.3% (n=92), 82.1% (n=234) and 23.63±4.73 years (n=235) in 1994 to 58.60% (n=113), 68.1% (n=191) and 26.44±6.43 years (n=193) in 2013, respectively. The gravidity for each woman decreased from 2.15±1.70 to 1.5±1.02 in the same years as well. shows the changes in structure , process and outcome indicators before and after the FP implementation. In terms of process indicators, the ratio of having at least one ultrasound during pregnancy had increased from 15.7% (n=37) in 1994 to 100% (n=185) in 2013. In contrast, the mean (±SD) number of ACV (9.13±3.67, n=232), the first measured Hb (13.22±1.08, n=139) and the ratio of normal deliveries (84.6%, n=198) had decreased to 6.43±1.80 (n=187), 12.83±0.99 (n=187) and 55.6% (n=105), respectively. In addition, the number of routine laboratory tests and the frequency of ultrasounds for each woman had increased from 25% and 0.14±0.37 in 1994 to 99.5% and 2.47±1.06 in 2013, respectively. In terms of outcome indicators, MMR, NMR and mean of BW decreased from 32.51 (deaths per 100 000 live births), 19.8 (deaths per 1000 live births) and 3429.70±495.33 g in 1994 to 21.09 (per 100 000 live births), 6.32 (per 1000 live births) and 3207.88±458.88 g in 2013, respectively. In contrast, the proportion of the history of abortion and stillbirth increased from 2.8% (n=2) and 1.4% (n=3) in 1993 to 14.5% (n=39) and 8.6% (n=7) in 2013, respectively. show the changes in process and outcome indicators. represent scatter plots of values of process and outcome indicators over the study time points. shows the results of ITSA (non-adjusted and adjusted for structure variables: mother’s education, occupation, age and gravidity) on selected indicators. For each specified response indicator, if the ‘intervention effect’ is significant, it shows that the average values of the indicator before and after the intervention times (2005) are statistically different. If the ‘time’ coefficient is significant, it means that there were changes in the response indicator over the study time regardless of the intervention time. If the ‘time–intervention interaction’ is significant, it shows that the slope of changes in the response indicator trend before and after the intervention time is different. Process indicators Antenatal care visits Over two decades covered by the study, the trend of ACVs had changed from a relatively stable trend to a sharp declining one. The independent impact of time and the programme on ACVs was not significant, but the interaction term (the slope of the decreasing trend) increased sharply after the intervention time both in unadjusted (b=−0.36, p<0.001) and adjusted models (b=−0.36, p<0.001). Routine laboratory tests and ultrasound The rate of routine laboratory tests, ultrasound and their repetitions during the pregnancy had increased over the years. There was a significant (in both unadjusted and adjusted models) independent association between time, the intervention and ultrasound (unadjusted b=3.28, p<0.001 and adjusted b=2.40, p=0.046 for time; unadjusted b=16.78, p=0.043 and adjusted b=19.31, p=0.049 for intervention). Their interaction term did not have a significant association with these indicators. In contrast, repetition of ultrasound had significant association with interaction term (unadjusted b=0.21, p<0.001 and adjusted b=0.20, p=0.008) and the unadjusted and adjusted independent effects of time (b=3.28, p<0.001). The independent impact of intervention on routine tests was not significant, but adjusted independent effect of time was significant (b=7.70, p=0.015). The slope of trend of routine test did not vary before and after the intervention. Normal delivery The ratio of normal vaginal deliveries showed a decreasing trend during the study years (unadjusted b=−0.94, p<0.001; adjusted b=−0.81, p=0.026), as well as before and after the intervention (unadjusted b=−7.41, p=0.008; adjusted b=−7.07, p=0.006). After controlling for the structural confounders, the significant difference in the slope of normal deliveries trend that was observed for the unadjusted model before and after the intervention (b=−1.91, p=0.035) disappeared in the adjusted model (b=1.33, p=0.266). Outcome indicators Mortality indicators According to the unadjusted model, NMR and MMR had declined over the study years (b=−0.77, p=0.002 and b=−0.95, p=0.004). After adjustment of structural confounders, this association, however, was true only for NMR (b=−0.56, p=0.021), but not for MMR (b=0.69, p=0.33). The independent effect of FP programme and its interaction with time did not show any significant association with these two indicators. Moreover, there was a significant increasing trend for the history of abortion in both unadjusted (b=1.54, p<0.001) and adjusted (b=1.15, p=0.009) models, regardless of the intervention. However, only the adjusted effect of intervention on the history of abortion was significant (b=4.28, p=0.047). No significant slope change was observed for the history of abortion in both adjusted and unadjusted models. Finally, the independent effect of time and intervention on the history of stillbirths was not significant over the study years. Nonetheless, the slope of the history of stillbirths trend showed a significant variation before and after the intervention in both unadjusted (b=0.84, p=0.006) and adjusted (b=1.72, p<0.001) models. Other outcome indicators For outcome indicators other than mortality indicators, BW and first measured Hb were significant associated with effects of time and the interaction term in both adjusted and unadjusted models. The coefficient for the interaction term in adjusted analysis was 33.2 (p<0.001), showing of a dramatic change in the slop of the BW trend before and after the FP programme. The maternal Hb level showed a decreasing trend during the study years in the unadjusted model (b=−0.04, p=0.007). There was no significant relation between Hb level and the intervention even after controlling for structural confounders. The joint effect of time and the intervention was not associated with Hb values in both adjusted and unadjusted models. 

indicators Antenatal care visits Over two decades covered by the study, the trend of ACVs had changed from a relatively stable trend to a sharp declining one. The independent impact of time and the programme on ACVs was not significant, but the interaction term (the slope of the decreasing trend) increased sharply after the intervention time both in unadjusted (b=−0.36, p<0.001) and adjusted models (b=−0.36, p<0.001). Routine laboratory tests and ultrasound The rate of routine laboratory tests, ultrasound and their repetitions during the pregnancy had increased over the years. There was a significant (in both unadjusted and adjusted models) independent association between time, the intervention and ultrasound (unadjusted b=3.28, p<0.001 and adjusted b=2.40, p=0.046 for time; unadjusted b=16.78, p=0.043 and adjusted b=19.31, p=0.049 for intervention). Their interaction term did not have a significant association with these indicators. In contrast, repetition of ultrasound had significant association with interaction term (unadjusted b=0.21, p<0.001 and adjusted b=0.20, p=0.008) and the unadjusted and adjusted independent effects of time (b=3.28, p<0.001). The independent impact of intervention on routine tests was not significant, but adjusted independent effect of time was significant (b=7.70, p=0.015). The slope of trend of routine test did not vary before and after the intervention. Normal delivery The ratio of normal vaginal deliveries showed a decreasing trend during the study years (unadjusted b=−0.94, p<0.001; adjusted b=−0.81, p=0.026), as well as before and after the intervention (unadjusted b=−7.41, p=0.008; adjusted b=−7.07, p=0.006). After controlling for the structural confounders, the significant difference in the slope of normal deliveries trend that was observed for the unadjusted model before and after the intervention (b=−1.91, p=0.035) disappeared in the adjusted model (b=1.33, p=0.266). 

Over two decades covered by the study, the trend of ACVs had changed from a relatively stable trend to a sharp declining one. The independent impact of time and the programme on ACVs was not significant, but the interaction term (the slope of the decreasing trend) increased sharply after the intervention time both in unadjusted (b=−0.36, p<0.001) and adjusted models (b=−0.36, p<0.001). 

The rate of routine laboratory tests, ultrasound and their repetitions during the pregnancy had increased over the years. There was a significant (in both unadjusted and adjusted models) independent association between time, the intervention and ultrasound (unadjusted b=3.28, p<0.001 and adjusted b=2.40, p=0.046 for time; unadjusted b=16.78, p=0.043 and adjusted b=19.31, p=0.049 for intervention). Their interaction term did not have a significant association with these indicators. In contrast, repetition of ultrasound had significant association with interaction term (unadjusted b=0.21, p<0.001 and adjusted b=0.20, p=0.008) and the unadjusted and adjusted independent effects of time (b=3.28, p<0.001). The independent impact of intervention on routine tests was not significant, but adjusted independent effect of time was significant (b=7.70, p=0.015). The slope of trend of routine test did not vary before and after the intervention. 

The ratio of normal vaginal deliveries showed a decreasing trend during the study years (unadjusted b=−0.94, p<0.001; adjusted b=−0.81, p=0.026), as well as before and after the intervention (unadjusted b=−7.41, p=0.008; adjusted b=−7.07, p=0.006). After controlling for the structural confounders, the significant difference in the slope of normal deliveries trend that was observed for the unadjusted model before and after the intervention (b=−1.91, p=0.035) disappeared in the adjusted model (b=1.33, p=0.266). 

Mortality indicators According to the unadjusted model, NMR and MMR had declined over the study years (b=−0.77, p=0.002 and b=−0.95, p=0.004). After adjustment of structural confounders, this association, however, was true only for NMR (b=−0.56, p=0.021), but not for MMR (b=0.69, p=0.33). The independent effect of FP programme and its interaction with time did not show any significant association with these two indicators. Moreover, there was a significant increasing trend for the history of abortion in both unadjusted (b=1.54, p<0.001) and adjusted (b=1.15, p=0.009) models, regardless of the intervention. However, only the adjusted effect of intervention on the history of abortion was significant (b=4.28, p=0.047). No significant slope change was observed for the history of abortion in both adjusted and unadjusted models. Finally, the independent effect of time and intervention on the history of stillbirths was not significant over the study years. Nonetheless, the slope of the history of stillbirths trend showed a significant variation before and after the intervention in both unadjusted (b=0.84, p=0.006) and adjusted (b=1.72, p<0.001) models. 

According to the unadjusted model, NMR and MMR had declined over the study years (b=−0.77, p=0.002 and b=−0.95, p=0.004). After adjustment of structural confounders, this association, however, was true only for NMR (b=−0.56, p=0.021), but not for MMR (b=0.69, p=0.33). The independent effect of FP programme and its interaction with time did not show any significant association with these two indicators. Moreover, there was a significant increasing trend for the history of abortion in both unadjusted (b=1.54, p<0.001) and adjusted (b=1.15, p=0.009) models, regardless of the intervention. However, only the adjusted effect of intervention on the history of abortion was significant (b=4.28, p=0.047). No significant slope change was observed for the history of abortion in both adjusted and unadjusted models. Finally, the independent effect of time and intervention on the history of stillbirths was not significant over the study years. Nonetheless, the slope of the history of stillbirths trend showed a significant variation before and after the intervention in both unadjusted (b=0.84, p=0.006) and adjusted (b=1.72, p<0.001) models. 

For outcome indicators other than mortality indicators, BW and first measured Hb were significant associated with effects of time and the interaction term in both adjusted and unadjusted models. The coefficient for the interaction term in adjusted analysis was 33.2 (p<0.001), showing of a dramatic change in the slop of the BW trend before and after the FP programme. The maternal Hb level showed a decreasing trend during the study years in the unadjusted model (b=−0.04, p=0.007). There was no significant relation between Hb level and the intervention even after controlling for structural confounders. The joint effect of time and the intervention was not associated with Hb values in both adjusted and unadjusted models. 

Following the improvement of rural inhabitants’ access to medical services, many studies predicted the improvement in MCH indicators as the main achievements of the FP plan in Iran. These improvements were expected mainly due to proper involvement of physicians and midwives in the antenatal care. It should be considered that the reductions in NMR and MMR after the provision of maternal and child health services in the first health system reform (implementation of the PHC network in 1985) in Iran was partly due to a higher rate of maternal and neonatal mortality in those years (high baseline rates). Therefore, the recent reform (FP programme) that upgraded the PHC network in rural area also targeted previously unachieved and unrealised MCH goals. The FP programme enriched the MCH services through provision of family healthcare services, improvement of services quality, and promotion of perinatal care and safe delivery in rural areas. This study shows that the FP programme can promote process indicators, as well as proximal outcome indicators (ie, BW), better than the mortality indicators (MMR and NMR). There was even some increase in the ratio of mothers who have experienced stillbirths and to some extent NMR after the implementation of the programme in 2005. Although this increase is likely to originate partly from a better registration of a fetal death following the intervention, there should be other determinants of maternal and neonatal health that are different from sheer availability and affordability of health services that the FP programme secured. The frequency of pregnancy ultrasound had increased over the years as well. This might be due to the free provision of these services after the programme implementation. These results are in line with other studies conducted in Cuba. At least one ultrasound session during gestation is recommended to diagnose abnormalities. However, easier access to physicians, along with less concern about financial bearings, had led to more use of healthcare services by pregnant women. Our findings also showed reductions in the number of vaginal deliveries. In fact, the FP programme failed to stop the increasing trend of caesarean sections among the study population. It should be borne in mind that the FP intervention was also likely to increase caesarean sections through the provision of better diagnostic services and over-medicalisation of pregnant women. This is in line with another study that has reported similar findings. Maternal Hb level is supposed to be influenced by nutritional and supplementary nutritional services. As some improvements were expected in the provision of micronutrients supplementation for pregnant mothers following the increase in the access to health services, an increase in the normal range of Hb levels was expected. However, our study did not show such a result. Interestingly, there was a critical change in the decreasing trend of BW, indicating that service availability and affordability is a strong determinant of this indicator. This reduction is one of the main positive effects of the FP programme in Iran. Nonetheless, there is no universal agreement on indicators of FP success. Many studies have referred to indicators such as MMR, NMR, infant mortality rate, childbirth by skilled and trained personnel, caesarean delivery rate, abortion ratio, stillbirths ratio and BW. According to these studies, a better access to health services can lead to a decrease in prenatal mortality, history of stillbirths and infant death rates during the first week delivery. It can also increase the survival chance of low BW infants. Although the reduction of MMR has been considered as a significant achievement of the FP reform by some studies, the programme did not have a significant effect on the decreasing trend of MMR in our study. It should be noted that any improvement in MMR trend is also highly dependent on social determinants of health. However, in countries with medium and low level of MMR, more targeted intervention should be considered. The findings of our study showed that the FP programme had no significant effect on NMR trend. This result differs from a previous study conducted in Iran. Another difference was also seen in the ratio of abortions, which increased after the programme in our study but was found to have a reducing trend in the previous studies. This difference can be partly explained by the differences in the methods of data collection. Previous studies were based on national and provincial statistics reports (aggregated data) and their time periods were shorter than the period in our study. They also did not adjust for possible confounders. Our study used data for a relatively longer period of time, covering almost 100% of its target population. A random sampling approach was used to select the study samples from recorded registries over the years. These characteristics make the study result less prone to selection bias. However, the effect of time-varying confounders could not be ruled out completely. In fact, having no concurrent control group to discriminate the effect of FP from changes that would have happened without FP is a major limitation of the study. Although the urban population (with more than 200 000 inhabitants) could have been considered as the comparison groups, the fundamental difference between these settings could produce a number of methodological challenges. Finally, as there is little difference between the features of the study sample and the characteristics of the representative sample reported in the latest national census in 2016, demographic health survey in 2011 and health system reports, the findings of this study can be generalised, with some level of uncertainty, to rural areas of the entire country. 

Our findings showed that introduction of the FP programme had a significant effect on process indicators and proximal outcome variables (ie, BW) in Iran. In addition, due to the accurate reporting and documentation of MCH records, the number of stillbirths and abortions significantly increased after the FP programme implementation. Interestingly, the effect of the reform on decreasing trends of MMR and NMR was not significant. It shows that along with some service availability factors that were addressed by the FP programme, there are some other effective factors that should be considered in order to improve the MCH indicators. 

Reviewer comments Author's manuscript 

