Modelling the Evolution of Domestic Violence Occurrences in Portuguese Municipalities - What Causes Domestic Violence?

i Modelling the Evolution of Domestic Violence Occurrences in Portuguese Municipalities Ana Clara do Carmo St. Aubyn What Causes Domestic Violence? Dissertation presented as partial requirement for obtaining the Master’s degree in Advanced Analytics ii NOVA Information Management School Instituto Superior de Estatística e Gestão de Informação Universidade Nova de Lisboa MODELING THE EVOLUTION OF DOMESTIC VIOLENCE OCCURRENCES IN PORTUGUESE MUNICIPALITIES by Ana Clara do Carmo St. Aubyn Dissertation presented as partial requirement for obtaining the Master’s degree in Advanced Analytics Advisor / Co Advisor: Mauro Castelli / Maria Jordão November 2021 iii ACKNOWLEDGMENTS Throughout the writing of this study, I faced a lot of challenges that could not have been overcome without the great deal of support and assistance I received. First of all, I would like to thank my supervisors, Professors Mauro Castelli and Maria Jordão for making the writing of this thesis possible and for all the support provided when things did not go the way I wanted them to. Secondly, I would like to thank my father for his precious support and readiness to enlighten me when I was in doubt. You were a fundamental piece of this thesis. I would also like to acknowledge my friends Sofia, Rita, David and Pedro, as well as my boyfriend for never letting me give up and for sharing the ups and downs of these two semesters. Finally, but not least important, I would like to thank my mother, who understands me like no other and always calms me down when needed. Thank you for always listening to me during the process of writing this thesis even if you do not have a slight knowledge of econometrics. iv ABSTRACT Throughout the last years, domestic violence has been a widely discussed topic and an essential concern when building healthy communities. To fight the prevalence of this problem, its causes must be addressed. These causes can come from personal indicators, but they can also come from issues in society at large. It is important to shift the debate from micro-level to macro-level analyses, looking for structural factors in societies that contribute to the evolution of the number of domestic violence occurrences as these are the factors that can be addressed by regulatory bodies. When modeling domestic violence one can envisage two types of possible explanatory variables: risk factors and protective factors. The first ones, as the term suggests, increase the risk of domestic violence, causing a high number of occurrences when very prevalent. Protective factors do the opposite, buffing the risk for domestic violence. The identification of risk factors is very important for the prevention of violence and to guide policies. However, identifying protective factors is of the utmost importance as their presence in Societies can be promoted to avert the occurrences. The present document studies the evolution of domestic violence occurrences in the municipalities of the Portuguese mainland between 2009 and 2019 resorting to panel data analysis methods. The constant coefficients and the fixed effects approaches are employed to try and understand the relation between possible causes and domestic violence occurrences. While explaining a good proportion of the total variance encapsulated in the dependent variable was revealed to be a hard task, evidence of the importance of some variables in explaining domestic violence occurrences on a macro-level was found. These variables were the average number of children born to each woman in fertile age, the number of divorces per 100 marriages, the percentage of resident population with normal age for attending high school that is actually attending high school, the number of new marriages per 100 inhabitants, the percentage of men’s monthly gain that women receive on average, the number of people enrolled in employment and vocational training centers per 100 inhabitants and the number of doctors per 100 inhabitants according to the Doctors’ Professional Order. Most of these were shown to be risk factors, increasing the number of domestic violence occurrences. KEYWORDS Domestic Violence; Econometrics; Panel Data; Constant Coefficients; Fixed Effects; Risk Assessment; Explanatory Modelling v INDEX 1. Introduction ............................................................................................................................. 1 1.1. Thesis Objective and Research Questions ........................................................................ 1 1.2. The Evolution of Domestic Violence in Portugal .............................................................. 2 2. Literature Review ..................................................................................................................... 4 3. Theoretical Background ........................................................................................................... 9 3.1. Panel Data......................................................................................................................... 9 3.2. Econometric Primer .......................................................................................................... 9 3.3. Causal Relationships and Ceteris Paribus ....................................................................... 13 3.4. Econometric Assumptions .............................................................................................. 14 3.5. Statistical Tests ............................................................................................................... 16 3.6. Heteroskedasticity .......................................................................................................... 18 3.7. Model and Variable Selection ........................................................................................ 19 3.8. Modelling Issues ............................................................................................................. 21 3.9. Models For Panel Data ................................................................................................... 22 4. Data Exploration ..................................................................................................................... 24 4.1. Dependent Variable ........................................................................................................ 25 4.2. Explanatory Variables ..................................................................................................... 28 5. Methodology .......................................................................................................................... 32 5.1. Series Breaks ................................................................................................................... 32 5.2. Correlations .................................................................................................................... 32 5.3. Missing Values ................................................................................................................ 34 5.4. Summary Statistics ......................................................................................................... 36 5.5. Outlier Analysis ............................................................................................................... 38 5.6. Stationarity ..................................................................................................................... 39 5.7. Model Estimation ........................................................................................................... 42 5.7.1. Model 1 (CC – Constant Coefficients) ..................................................................... 43 5.7.2. Model 2 (CC) ........................................................................................................... 45 5.7.3. Model 3 (CC) ........................................................................................................... 47 5.7.4. Model 4 (CC) ........................................................................................................... 48 5.7.5. Model 5 (FE – Fixed Effects) .................................................................................... 49 5.7.6. Model 6 (FE) ............................................................................................................ 50 5.7.7. Model 7 (FE) ............................................................................................................ 52 5.7.8. Model 8 (FE) ............................................................................................................ 53 5.7.9. Model 9 (FE) ............................................................................................................ 53 6. Results and Discussion ........................................................................................................... 55 vi 7. Conclusions ............................................................................................................................ 57 8. Limitations and Recommendations for Future Works ........................................................... 59 9. Bibliography ........................................................................................................................... 60 10. Annexes .................................................................................................................................. 62 vii LIST OF FIGURES Figure 1.1 - Domestic Violence Occurrences (on a national level) .......................................................... 2 Figure 1.2 - Domestic Violence Occurrences by Category (on a national level) ..................................... 3 Figure 1.3 - Domestic Violence Against Spouse or Analogous Occurrences (on a national level) .......... 3 Figure 2.1 - Age of Domestic Violence Victims (2013-2017) ................................................................... 8 Figure 3.1 - Probability Density Function for DVASA ............................................................................ 10 Figure 3.2 - Conditional Probability Density Function for DVASA (GER=70,16) .................................... 11 Figure 3.3 - Graphical Relation Between GER and DVASA .................................................................... 11 Figure 3.4 - Choosing Variables and Functional Forms for Models ....................................................... 19 Figure 3.5 - Diagnostic Residual Plots (Example) .................................................................................. 21 Figure 4.1 - Absolute Change in Total and DVASA Occurrences (on a national level) .......................... 26 Figure 4.2 - Are Marriages and DVASA Occurrences Related? .............................................................. 30 Figure 5.1 - Average Contemporaneous Pearson Correlation .............................................................. 33 Figure 5.2 - Histogram for DVASA Distribution ..................................................................................... 37 Figure 5.3 - Evolution of GER, GER_Men, and GER_Women in Barrancos ........................................... 39 Figure 5.4 - Joint Distributions of Explanatory Variables with DVASA .................................................. 42 Figure 5.5 - Heteroskedasticity Test for Model 1 .................................................................................. 44 Figure 5.6 - Residuals Against Explanatory Variables (Model 1) ........................................................... 45 Figure 5.7 - Heteroskedasticity Test for Model 2 .................................................................................. 46 viii LIST OF TABLES Table 4.1 - Dataset Variable Description ............................................................................................... 25 Table 4.2 - Descriptive Statistics for the Dependent Variable (on a national level) ............................. 26 Table 4.3 - Missing Values by Domestic Violence Category (Municipalities) ........................................ 27 Table 4.4 - Missing Values for DVASA in Municipalities and Difference Between National Total and Municipality Total ......................................................................................................................... 27 Table 5.1 – Overall Summary Statistics ................................................................................................. 36 Table 5.2 - Overall, Between and Within Standard Deviations ............................................................. 38 Table 5.3 - Parameter Estimates for Model 1 ....................................................................................... 43 Table 5.4 - Parameter Estimates for Model 2 ....................................................................................... 46 Table 5.5 - Parameter Estimates for Model 3 ....................................................................................... 47 Table 5.6 - Parameter Estimates for Model 4 ....................................................................................... 49 Table 5.7 - Parameter Estimates for Model 5 ....................................................................................... 50 Table 5.8 - Parameter Estimates for Model 6 ....................................................................................... 51 Table 5.9 - Parameter Estimates for Model 7 ....................................................................................... 52 Table 5.10 - Parameter Estimates for Model 8 ..................................................................................... 53 Table 5.11 - Parameter Estimates for Model 9 ..................................................................................... 54 Table 6.1 - Model Comparison .............................................................................................................. 55 ix LIST OF ABBREVIATIONS AND ACRONYMS AIC Akaike Information Criterion APAV Associação Portuguesa de Apoio à Vítima CC Constant Coefficients DGEEC Direção-Geral de Estatísticas da Educação e Ciência DGPJ Direção-Geral da Política de Justiça DVAM Domestic Violence Against Minors DVASA Domestic Violence Against Spouse or Analogous FE Fixed Effects GAM Generalized Additive Model GBV Gender Based Violence GER Gross Enrolment Rate IPV Intimate Partner Violence KNN K Nearest Neighbors OECD Organization for Economic Co-operation and Development OMA Observatório de Mulheres Assassinadas da UMAR RESET Regression Specification Error Test SC Schwarz Criterion SFI Synthetic Fertility Index WHO World Health Organization YDI Youth Dependency Index 7 factor, as around 34% of the victims were married. This is a relevant percentage when compared to the 20.8% that were single, 16% whose marital status was unknown, 11.6% who were in a non-marital relationship, 8.7% who were divorced, 5.6% who were separated and 3.3% who were widowed. Therefore, including a measure of the number of married people may be relevant. Still regarding the marital status, (Bowlus & Seitz, 2006) finds that women who are severely abused by their husbands are significantly more likely to divorce than women who do not face this problem. However, it is important to notice that women may be more likely to report violence in a past marriage than in a present one, causing an upward bias in this probability. This study uses data from all provinces of Canada, retrieved in 1993. The initial analysis of this data also revealed that women who experienced abuse by their partners tend to have lower levels of education and come from more violent backgrounds than women who did not face abuse. The same applies to the partners – husbands who abuse their wives tend to have lower levels of education. Another finding from this initial analysis is that women who are not working are more likely to face abuse, which meets the conclusions on (Anderberg, Rainer, Wadsworth, & Wilson, 2015). According to (WHO - World Health Organization, 2010), divorces may not only be a consequence, but also a cause of domestic violence. People who are separated or divorced tend to be more vulnerable, increasing their probability of becoming victims in a future relationship. This same report by WHO lists some of the causes for domestic and sexual violence included in the literature they reviewed. Young age appears to be a risk factor for either becoming a victim of intimate violence or a perpetrator. It is foreseeable that populations with a higher proportion of young adults have higher rates of domestic violence occurrences. Lower levels of education are also consistently associated with both sides of the crime (victim and perpetrator). A higher level of education may act as a protective factor, since people with a higher level of education show lower levels of intimate partner violence. Another factor associated with both the victim and the perpetrator is poverty. Even though domestic violence cuts across all socioeconomic groups, people with lower incomes tend to be more at risk of becoming either a victim or an aggressor. One explanation for this, besides the hopelessness, stress and frustration caused by this condition, is the fact that shortage of money is a common cause for marital arguments and makes it harder for people to leave toxic relationships, as they are more financially dependent on each other. Some characteristics of the neighborhood may also influence the number of IPV occurrences, such as a lower proportion of women with higher levels of education, higher unemployment rates, a higher proportion of illiteracy and a lower proportion of women with high levels of autonomy. (Ackerson, Kawachi, Barbeau, & Subramanian, 2008) examined the role of women’s education and proximate educational context on GBV in India. A sample of 83,627 married women aged 15 to 49 years old from the 1998 to 1999 Indian National Family Health Survey was examined. The study considered that not only does the level of education of the woman herself influence her probability of becoming a GBV victim, but also does the general level of education of the community surrounding her. The results of this paper show that women with no education are 4.5 times more likely to report having suffered from domestic violence at some point in their life than women schooled for more than 12 years. Another relevant conclusion was that the probability for a woman who is living in the middle and lowest tertiles of female literacy to suffer from domestic violence at some point in her life was 1.18 and 1.10 times greater, respectively, than those of women living in the highest tertile neighborhoods. This study shows the impact that education has on GBV. 8 Figure 2.1 - Age of Domestic Violence Victims (2013-2017) 9 3. THEORETICAL BACKGROUND The theoretical background of the present study was mostly written considering (Hill, Griffiths, & Lim, 2012) and (Wooldridge, 2013). 3.1. PANEL DATA Data can be collected in multiple formats. The most widely discussed ones are cross-sectional data, pooled cross-sectional data, time series data and panel data. Cross-sectional data is data collected for multiple units across the same period. Each observation represents a unit of the relevant population. This is the “common” dataset structure. When we combine cross-sectional data from different periods we create a pooled cross-sectional dataset. In this case, each observation represents a unit of the population in a specific period in time. It is not necessarily true that the same units are studied for the different periods. If we are studying the same unit across different periods in time, we create a time series. A time series shows the evolution of that unit through a specified time span. Finally, panel data, also called longitudinal data, is a combination of cross-section and time series data. Here, we have one time series for each included unit. The identifier of the unit and the period the data refers to are shown as variables in the dataset. The present study focuses on panel data, as there is one yearly discrete time series for each Portuguese municipality. Panel data analysis is a way of studying a subject in multiple sites periodically observed over a time frame. Panel data may be considered short or long, balanced or unbalanced and fixed or rotating. A panel data is considered short when it studies many units for a short time period and it is considered long when it studies few units for a long time period. In the case of the present study, a short panel is being examined as it has 11 time periods and 278 units. A balanced panel has a number of observations equal to the number of units times the number of periods, meaning that all units are observed for all periods. If this is not the case, we have an unbalanced panel. The present study focuses on a balanced panel as all 278 municipalities are observed for all 11 years. Finally, if the same units are observed for each period, the panel data is fixed. If the set of units varies from one period to the next the panel data is rotating. In the case of this study, the 278 municipalities are fixed. 3.2. ECONOMETRIC PRIMER Econometry starts with a theory about how some relevant variables are related to others. To express our ideas regarding these relationships we use functions. For most problems, it is not enough to know in which direction the variables are related (if they increase together, vary in opposite directions, etc.). Instead, one needs to know the intensity of that relation, which means how much a change in the value of one variable will affect the value of the other. To know these parameters of the relationship we create regressions. Before generating an econometric model, one must keep in mind that relations among variables are not exact and, for this reason, no econometric model explains the exact behavior of its object of study. Instead, it describes the average or systematic behavior. This means that when presenting results from an econometric model one must always say “it is expected that…”. Since the predictions are not exact, there is a difference between the actual value and the predicted value. This difference is the random 10 component of the model and is called the error term. It represents all factors that were not included in the model and includes the behavior uncertainty. The part that is explained by the model is the systematic component of the formula and it is decided by the investigator based on what the theory already existent states about the problem that is being studied. One can say that the error term represents all things affecting the dependent variable other than the explanatory variables included in the model. It comprises the effects of relevant variables that were excluded from the model, measurement errors in both the dependent and explanatory variables, the effects of using a linear or any other form to generate results that do not follow that form and, finally, the natural randomness of observations. To complete a model specification, the researcher must choose which variables to include and how to include them, keeping in mind the algebraic form of the relations. This form is also called the functional form and in the most basic scenario, it is assumed to be linear. An econometric model looks like the formula below where Y is the variable being studied (dependent variable), Xi are the variables that are assumed to affect the behavior of Y (explanatory variables), βi are the coefficients of the explanatory variables that indicate how much they affect the dependent variable and, finally, ei is the error term. β0 is the expected value of Y if all explanatory variables are 0 and is called the constant of the model. 𝑌𝑖= 𝛽0+𝛽1𝑋1𝑖+𝛽2𝑋2𝑖+𝛽3𝑋3𝑖+𝑒𝑖 As it can be seen in the equation above that represents a multiple regression linear model, the model relates a dependent variable Y to a set of explanatory variables (X1, X2 and X3) and to a random error term e. This remains true for other econometric models, changing the set of explanatory variables, the functional form, etc. The βs are the parameters that are estimated by the model and they can be interpreted to explain the relationships among variables. In the example above, one can say that a change in X1 of one unit will represent a change of β1 units in Y, holding everything else constant. Taking this into account, one can start to understand that the values calculated using an econometric model are, in fact, a conditional average. The predicted value of Y is the average of Y if the explanatory variables assume a set of fixed values. Let us consider the example of the present study: the probability density function for Domestic Violence Against Spouse or Analogous (DVASA) is illustrated in Figure 3.1 below, with a gray dashed line indicating the average value. Figure 3.1 - Probability Density Function for DVASA Systematic Component Random Component 11 If we have a simple linear regression model, with only one explanatory variable (let us consider Education as the only explanatory variable), the predicted value for DVASA in a given point would be the average value of DVASA when Education attains a certain value. Therefore, we can calculate a conditioned probability density function and find its average to find the predicted value. Figure 3.2 below shows the plot for the probability density function of DVASA conditioned to when the variable GER (a measure of education – further explanation on Chapter 4.2) is at its most common value. One can see that the average value of this new distribution shifted slightly to the right when compared to the general distribution of DVASA. These are the differences that will be reflected in an econometric model. Figure 3.2 - Conditional Probability Density Function for DVASA (GER=70.16) There are numerous types of econometric models. The simplest one takes only one explanatory variable and is called the simple linear regression model. If we include more than one explanatory variable, we are creating a multiple regression model. Econometrics also contemplates models for time series and panel data. To understand how an econometric model is estimated the best option is to first study the simple linear regression model. In the previous example we considered DVASA as our dependent variable and GER as the only explanatory variable. We can plot the relationship between these variables to see how good of an approach we get. Figure 3.3 below shows a scatter plot of the two variables on the left. On the right side of the same figure, we added an approximation line to predict DVASA based on GER. However, adding a line is no guarantee that the expected value of the errors is 0, so the best approach is to use the Least Squares Estimators. Figure 3.3 - Graphical Relation Between GER and DVASA 12 To estimate the intercept and the slope of the optimal line that describes the relationship between the dependent variable and the explanatory variable we want to make use of all observations available. The least squares principle says that we should add the line in a way so that the sum of squares of the vertical distances between the observations and the fitted line is minimized. It is important to square these distances so that positive distances do not cancel negative ones. These vertical distances are the residuals, hence the desire to minimize them. The sum of squares we wish to minimize is given by: 𝑆(𝛽0,𝛽1)= ∑(𝑦𝑖−𝛽0−𝛽1𝑋1𝑖)2 𝑛 𝑖=1 The function above is quadratic in terms of β0 and β1 and is shaped like a bowl. To find its minimum value we need to find the bottom of the bowl, which occurs where the slope of the bowl in the direction of each axis is 0. This is the same as saying that it occurs where the partial derivatives of S concerning β0 and β1 are 0. By calculating these partial derivatives, we obtain: 𝜕∑(𝑦𝑖−𝛽0−𝛽1𝑥𝑖)2 𝑛 𝑖=1 𝜕𝛽0=0 ⇔𝛽0= 𝐸(𝑦𝑖)−𝛽1𝐸(𝑋1𝑖) 𝜕∑(𝑦𝑖−𝛽0−𝛽1𝑥𝑖)2 𝑛 𝑖=1 𝜕𝛽1=∑𝑦𝑖𝑋1𝑖 − 𝛽0∑𝑋1𝑖 −𝛽1∑𝑋1𝑖 2 We can then replace β0 as given by the first equation in the second one and obtain: ∑𝑦𝑖𝑥𝑖−(𝐸(𝑦𝑖)−𝛽1𝐸(𝑥𝑖))∑𝑥𝑖−𝛽1∑𝑥𝑖2=0 ⇔𝛽1=𝐸(𝑦𝑖)𝐸(𝑥𝑖)−∑𝑦𝑖𝑥𝑖 𝑛 𝐸(𝑥𝑖)2−∑𝑥𝑖2 𝑛 These are the formulas for the least-squares estimators in the simple linear regression model and they can be used regardless of what the dependent or explanatory variables are. Even though the simple linear regression model is the easiest to understand, most real-life econometric models take two or more explanatory variables, creating multiple regression models. Most of the conclusions regarding the simple linear regression model can be adapted for this type of model, except for some changes in the interpretation of the coefficients and in the degrees of freedom of the T distributions. In a multiple regression model, we have several explanatory variables and we want to quantify the impact of each over the dependent variable, while controlling the effects of the remaining ones. This is the notion of ceteris paribus (see Chapter 3.3). Let us continue using the same example where DVASA is the dependent variable and GER is the explanatory variable, except that this time we are adding unemployment as another explanatory variable. We then obtain the following theoretical model: 𝐷𝑉𝐴𝑆𝐴𝑖= 𝛽0+𝛽1𝐺𝐸𝑅𝑖+𝛽2𝑈𝑛𝑒𝑚𝑝𝑙𝑜𝑦𝑚𝑒𝑛𝑡𝑖+𝑒𝑖 The βs in the above model measure the change in the dependent variable DVASA given a change of one unit in the respective explanatory variable, while holding everything else constant. The expected impact of an explanatory variable in the dependent variable is given by its partial derivate: 𝜕𝑦 𝜕𝑥 13 The objective in a multiple regression model is the same as it was in the simple regression model – to minimize the sum of squares of the vertical distances between the observations and the fitted line. The difference is that this time we need to work with matrices to find the least-squares estimators. In matrix notation, the model is given by: 𝑌=𝑋𝛽+𝑒 In the previous equation Y is a matrix with only 1 column and n (n being the number of observations) rows containing the observed values for the dependent variable. X is a matrix with n rows and k+1 columns (where k is the number of explanatory variables). The first column in the X matrix only contains ones and will be used to calculate β0. The remaining columns in X contain the observed values for the explanatory variables and will be used to calculate the respective βs. The matrix β has k+1 rows and only one column, containing the values of the least-squares estimators. Finally, the matrix e has n rows and only one column containing the residuals for each observation. In this case, minimizing the sum of squares of the residuals is the same as minimizing the following (where the T means transposed): ∑ 𝑒𝑖2 𝑛 𝑖=1 = 𝑒𝑇𝑒=(𝑌−𝑋𝛽󰆹)𝑇(𝑌−𝑋𝛽󰆹)=𝑌𝑇𝑌−2𝑌𝑇𝑋𝛽󰆹+𝛽󰆹𝑇𝑋𝑇𝑋𝛽󰆹 To determine the values for the estimators that minimize the sum of squares we need to differenciate the above expression and equal it to 0, as we did with the simple regression model. We then obtain: 𝜕(𝑌𝑇𝑌−2𝑌𝑇𝑋𝛽󰆹+𝛽󰆹𝑇𝑋𝑇𝑋𝛽󰆹) 𝜕𝛽󰆹=0⇔𝛽󰆹=(𝑋𝑇𝑋)−1𝑋𝑇𝑌 3.3. CAUSAL RELATIONSHIPS AND CETERIS PARIBUS The goal in regression analysis is to find causal relationships between variables, that is, to determine whether a change in X causes a change in Y. To express our ideas regarding the relationships between variables we use functions. For example, to express a relationship between domestic violence and education one can write: 𝐷𝑜𝑚𝑒𝑠𝑡𝑖𝑐_𝑉𝑖𝑜𝑙𝑒𝑛𝑐𝑒=𝑓(𝐸𝑑𝑢𝑐𝑎𝑡𝑖𝑜𝑛) The previous function is a possible notation to say that the domestic violence occurrences in a certain place are a function of the education level of the people residing in that same place. This is the same as saying that the places where domestic violence occurs depend on the prevalent educational level in those places. However, the occurrences may not depend only on education, they may depend on many other factors as well. Keeping this in mind, to understand the relationship between domestic violence occurrences and education it is important to set aside the impacts of the remaining factors on domestic violence. The idea of ceteris paribus (c.p.) means to hold all other factors constant and is a key point in establishing causal relationships. Without holding the remaining variables constant one does not show that the change observed in y is caused by the change in x. This is also the reason why a simple correlation study is not enough to analyze causal relationships. When we are studying the causal effect of x on y, the remaining variables that influence y are called the control variables. The reason to control for these variables is simple: we believe that x is correlated with other factors influencing y, which means that not holding these variables constant will make their 14 effects reflect on the coefficient for x. Since we feed data to the regression model, it is important to correctly determine the control variables that need to be held fixed. This is a critical part of regression analysis but may be hard, as usually not all factors influencing the dependent variable are observable or accessible. When one does not include an important control variable, its effects are reflected on the partial effects of other factors, making the latter incorrect. Stating the difference between explanatory variables and control variables may be hard. Even if some variables can be considered control variables on all occasions, all explanatory variables are eventually control variables when it comes to explaining the partial effects of another variable. For panel data, the variables that determine the difference between observations (in the case of this study: municipality and year) are always control variables as, even if they can explain part of the variance in the data, their main goal is to distinguish between observations. 3.4. ECONOMETRIC ASSUMPTIONS In every econometric study, there are at least two models: the theoretical model and the empirical one. The theoretical model describes a behavior but is an abstraction of reality. To convert the theoretical model into an empirical one, some assumptions must be made. These assumptions are very important because if they are verified our conclusions are warranted. However, if our model does not meet them, the conclusions may not be true. The first assumption is that the variance of the dependent variable is a constant for each value of X. This means that in the case of the example using DVASA as the dependent variable and GER as the explanatory variable, for each value that GER may take, the average of DVASA may be different, but the measure of how much it varies around that average must be the same. This leads to a further assumption that the error term must have a constant variance as well. When this condition is satisfied, the data is said to be homoskedastic. When this condition is violated, the data is said to be heteroskedastic. The second assumption is that the dependent variable is not only random but also statistically independent. This means that the value for one observation is not dep,endent on the value of another observation. In the DVASA example this means that the value of DVASA for one municipality is not dependent on the values for other municipalities on the dataset. Instead of assuming the independence of the variables, we often assume that the covariance of the variable is 0. Again, in the DVASA example, this means that if 𝑌𝑗 and 𝑌𝑖 are the values of DVASA for two different municipalities, 𝑐𝑜𝑣(𝑌𝑗,𝑌𝑖)=0. The third assumption is also related to variance. Since the main goal of an econometric study is to understand how changes in the explanatory variables affect the dependent variable, it is important that the explanatory variables are scattered enough. Obviously, if there is no variance in the explanatory variables it is impossible to justify changes in the dependent variable from changes in the explanatory variables. Therefore, it is assumed that the explanatory variables take at least two values. The fourth assumption is related to the error term of the regression. Since the total error of the model is the sum of the deviations between the actual value of an observation and the predicted value for it and the objective of a regression is to minimize these deviations, it is assumed that the expected value of the error ,term is 0. Furthermore, in a simple regression model it is also assumed that the expected 15 value of the error term given the explanatory variable is also 0. This is demonstrated in the equation below. It is important to remember that the expected value of the dependent variable given the explanatory variable is the predicted value so 𝐸(𝑦|𝑥)=𝛽0+𝛽1𝑋1. 𝐸(𝑒|𝑥)=𝐸(𝑦|𝑥)−𝛽0−𝛽1𝑋1=0 Finally, sometimes it is also assumed that the error term follows a normal distribution centered around 0. This is an optional assumption, but it is a strong one as the probability distribution of the parameters estimated in the regressions (βs) depends on the distribution of the error term which means that this is an important assumption for statistical analysis of the model. However, if there are enough observations one can base the conclusions on the central limit theorem that establishes that if the sample is large enough, the distribution of the sample means will be approximately normally distributed. These assumptions exist so that we can determine the quality of the least-squares estimators. These estimators are supposed to be centered and efficient. An efficient estimator is one with minimal variance. When the expected value of an estimator equals the real value of the parameter it is trying to estimate we say that we are dealing with a centered estimator. This means that the expected values of the estimators for the βs must be equal to the βs. In the simple regression model this would mean that: 𝐸(𝛽1 )=𝛽1 𝑎𝑛𝑑 𝐸(𝛽0 )=𝛽0 These equations are only true if the expected value of the error is 0, the error is not correlated with the variables and the data is coming from a random sample. The variance of an estimator is also key to evaluating its reliability. It measures how much the values that the estimator may take are spread and, because of it, ends up measuring the precision of the estimator. Keeping this in mind, the smaller the variance of an estimator, the more precise it will be. For the simple regression model, when we calculate the variance of the estimators we get: 𝑉𝑎𝑟(𝛽1 )= 1 ∑(𝑥𝑖−𝑥)2𝜎2 𝑎𝑛𝑑 𝑉𝑎𝑟(𝛽0 )=𝜎2(∑𝑥𝑖2 𝑛∑(𝑥𝑖−𝑥)2) Notice that σ2 stands for the variance of the error term and appears on both expressions. We can see that the larger the variance of the error term, the larger the variance of both estimators and, consequently, the more imprecise the estimation. We can also conclude that for the estimator to be efficient, the variance of the error term must be constant. Furthermore, we can see that the sum of the squared distances of x to its average is on both expressions on the denominator. This means that the larger the dispersion of the values in x, the smaller the variance of estimators, thus, the more precise they are. This means that choosing a sample that is diverse in terms of the values in the explanatory variables contributes to a higher precision of estimators. An important theorem regarding the assumptions in an econometric model is the Gauss-Markov theorem. It states that if the fundamental assumptions are verified, the estimators for the βs are the ones with the least variance out of all the centered estimators. This means that they are the BLUE (Best Linear Unbiased Estimators). 16 When one is applying a multiple regression model, a new assumption must be verified. That is that none of the explanatory variables is a perfect linear combination of another. When this assumption is not verified, we are facing perfect multicollinearity and the least-squares estimators cannot be calculated. 3.5. STATISTICAL TESTS When we estimate the coefficients for a regression, we are making a punctual estimate for the regression parameters. These estimates represent an inference over the regression model because after we calculate the values of the parameters for our sample, we intend to make an inference for the results to be applied to the population. To make inferences we resort to interval estimation and hypothesis tests. Both these procedures are strongly based on the assumption that the residuals in the model follow a normal distribution. If this assumption is not verified, it is necessary to guarantee that the sample is large enough to assure that the least-squares estimators follow an approximately normal distribution through the central limit theorem. When resorting to this theorem the tests and intervals can be calculated, but their results are only approximate. A hypothesis test allows us to evaluate the possibility that a parameter is equal to some value. In each hypothesis test, there must be a null hypothesis, an alternative hypothesis, a test statistic, a critical region, and a conclusion. The null hypothesis is denoted as H0 and most of the time equals the parameter being tested to a specific value. It is the hypothesis that we intend to accept or reject according to the statistical evidence. Paired to any null hypothesis there is an alternative one, denoted as H1. It is usually the contrary of the null hypothesis and is the one that is accepted if the null hypothesis is rejected. To decide whether the null hypothesis is accepted or rejected, we must take into consideration the value of the test statistic. The distribution of this statistic is known if the null hypothesis is true but is unknown otherwise. The critical region or rejection region of a hypothesis test depends on the form of the alternative hypothesis. It consists of the values that have a truly low probability of occurring if the null hypothesis is true. The logic behind this is that if the test statistic falls into the rejection region it is very unlikely that the null hypothesis is true. To define this region we must define a level of significance for the test. The level of significance of a test is the probability of rejecting the null hypothesis while it is actually true. This is also called a type I error. The level of significance is commonly designated by α and is usually either 0.01, 0.05 or 0.1. Another important concept when mentioning the significance of a test is the p-value. The p-value is the minimal level of significance with which we can reject the null hypothesis. If the p-value is less than the determined level of significance we reject the null hypothesis. It is important to always remember while performing a hypothesis test that these tests are not able to prove that a null hypothesis is true or false. We can only conclude if the data is compatible with that hypothesis or not. One of the most important hypothesis tests when it comes to econometric models is about the individual significance of the estimated parameters. For this, we test the hypothesis of a parameter being 0, which would mean that there is no significant relationship between the dependent variable 23 impossible if the number of coefficients to be estimated surpasses the number of observations for each unit. One way to avoid this problem is by using an alternative version of the fixed effects model, in which only the constant in the model differs from unit to unit. This can be written, considering three explanatory variables, as following: 𝑦𝑖𝑡 = 𝛽0𝑖 +𝛽1𝑥1𝑖𝑡 +𝛽2𝑥2𝑖𝑡 +𝛽3𝑥3𝑖𝑡 +𝑢𝑖𝑡 One can see that in the previous equation only the first β has the subscript i, meaning that all differences between individuals are assumed to be captured by the intercept in the model. One way to estimate this simple version of the fixed effects model is by crating a dummy variable for each unit that serves as an identifier for which individual the observation belongs to. By assuming the value 0 or 1 according to whether the observation refers to the unit or not, the corresponding coefficient will be ignored for observations that do not belong to the respective unit (by being multiplied by 0). In order to avoid problems related to multicollinearity, one unit must serve as the base for the others, which means that one unit does not have a corresponding dummy variable and is then represented by a constant in the model. The remaining coefficients represent the differences in the constant regarding the base unit for the remaining individuals. It is important to notice that, in short panels, the coefficients for the dummy variables are calculated using as many observations as the number of times periods contemplated for each unit, meaning that probably the central limit theorem does not apply anymore. Keeping this in mind, inferences on the coefficients for the dummy variables need normally distributed errors in order to be valid. Fixed effects models can have different coefficients for the cross-sectional units, for the time periods or for both. The base idea is always the same – to reflect the fixed effects of each of these components differently. 24 4. DATA EXPLORATION The dataset is composed by 19 explanatory variables and the dependent variable (DVASA occurrences). The methods used for calculating, treating, and standardizing all of these variables are explained in the next pages of this study. The dataset follows a panel data structure, having a column for the spatial dimension (Municipality) and one for the temporal dimension (Year). A summary of the existing variables can be found on Table 4.1, below: Variable Name Description Municipality Spatial Dimension of the Dataset. Year Temporal Dimension of the Dataset. DVASA Number of DVASA Occurrences Registered by Police Authorities by 100 Inhabitants. Divorces Number of Divorces for 100 Marriages in that Civil Year. Elderly_Dependency Number of People Aged 65 and Over for Every 100 People of Working Age, that is, Between 15 and 64 Years Old. Female_Doctors Percentage of Doctors Enrolled in the Doctor’s Order Who Are Female. Fertility Average Number of Children Born for Each Woman in Fertile Age (Between 15 and 49 Years Old). GER Percentage of the Resident Population with Normal Age for Attending High School that is Actually Attending High School. GER_Men Percentage of the Male Resident Population with Normal Age for Attending High School that is Actually Attending High School. GER_Women Percentage of the Female Resident Population with Normal Age for Attending High School that is Actually Attending High School. Marriages Number of New Marriages per 100 Inhabitants. Men65 Percentage of the Total Population of the Municipality that Represents Men With 65 Years or More. Mental_Health Percentage of Total_Doctors that are Specialized in Psychiatry according to the Doctor’s Order. Middle_Aged_Women Percentage of the Total Population of the Municipality that Represents Women Between 25 and 54 Years Old. Monthly_Gain Average Gross Amount that the Employees in the Municipality Receive Every Month Including basic remuneration and Other Remuneration Paid by the Employer (Overtime, Holiday Pay or Premiums). SS_Pensions Number of Pensioners for Each Person who Cashes for Social Security. A 25 Pension is an Amount Attributed Each Month to Someone in the Event of Disability, Old Age, Occupational Disease or Death. Total_Doctors Number of Doctors by 100 Inhabitants According to the Doctor’s Order. Unemployment_Female Number of Women Enrolled in Employment and Vocational Training Centers per 100 Inhabitants. Unemployment_Male Number of Men Enrolled in Employment and Vocational Training Centers per 100 Inhabitants. Unemployment_Total Number of people Enrolled in Employment and Vocational Training Centers per 100 Inhabitants. Youth_Dependency Number of Children Under 15 Years Old for Every 100 People of Working Age, that is, Between 15 and 64 Years Old. Wage_Gap Percentage of Men’s Monthly Gain that Women Receive on Average. Table 4.1 - Dataset Variable Description 4.1. DEPENDENT VARIABLE Three datasets containing information regarding the dependent variable were retrieved from the official statistics website by DGPJ on the 4th of March of 2021. One containing information about the number of domestic violence occurrences nationwide, one with this data split by districts and a last one with the data split by municipalities. All of them contained information regarding three categories (as explained in Chapter 1. Introduction): domestic violence against spouse or analogous (DVASA), domestic violence against minors (DVAM) and others. Finally, for all three datasets, the data was collected for the period between 2008 and 2019, due to data availability. The Portuguese Criminal Code provides for and punishes the crime of domestic violence. Domestic violence assumes the nature of a public crime, which means that the criminal procedure is not dependent on a complaint by the victim, just a complaint or knowledge of the crime is enough for the Public Ministry to promote the process. Thus, in Portugal, the registered number of occurrences of domestic violence does not depend only on self-report by the victim. However, as it is a crime that commonly takes place in the privacy of a home, many cases may depend on self-report. The dataset obtained focuses on data registered by police authorities and, according to (Ellsberg, Heise, Peña, Agurto, & Winkvist, 2001), may suffer from underreporting, as it depends on self-report to some extent. The data recorded for Portugal as a whole is a discrete time series. For each category there is a set of 12 observations recorded at uniformly spaced time values, in this case, years. This remains true for the data regarding districts and municipalities, except that, for the first case, there is one time series per category and per district, and for the second case there is one time series per category and per municipality. The evolution of the number of domestic violence occurrences in Portugal can be seen above in Figures 1.1 and 1.2. It becomes clear by the analysis of these figures and of the descriptive statistics on Table 26 4.2 that domestic violence against spouse or analogous is the most prominent category out of the three. One can see that, between 2008 and 2019, the yearly average of domestic violence occurrences was 27,394. Considering the same period, the yearly average for the DVASA category was 22,977.8, a value that clearly shows how relevant this category is for the total domestic violence occurrences. The remaining categories have less significant yearly averages. DVASA DVAM Others Total Occurrences Std 1226.32 75.22 435.55 1612.36 Minimum 20394.00 430.00 3083.00 24157.00 Mean 22977.80 537.92 3879.17 27394.90 Maximum 25129.00 680.00 4651.00 30340.00 Q3 23382.80 599.00 4039.00 27877.00 Median 22851.50 515.50 3800.00 27155.00 Q1 22457.50 484.00 3647.75 26683.50 Table 4.2 - Descriptive Statistics for the Dependent Variable (on a national level) When comparing the evolution of total occurrences in Portugal (Figure 1.1) with the evolution of occurrences for DVASA (Figure 1.3), one can detect the same patterns. By calculating the difference in the number of occurrences for subsequent years, it is possible to notice that DVASA occurrences almost always justify over half of the growth or decrease in the number of total occurrences. This is only not true for 2014, when the number of DVASA occurrences increased by 35 but the number of total domestic violence occurrences decreased by 48 due to a decrement in the other categories. Figure 4.1 shows exactly this. Figure 4.1 - Absolute Change in Total and DVASA Occurrences (on a national level) Since the datasets only have 12 years worth of data, it would not be possible to perform a time series regression for Portugal as a whole, as there would not be enough degrees of freedom to provide powerful estimates. Keeping this in mind, a panel data regression will be performed with data regarding the years and municipalities. For this purpose, the dataset containing information about domestic violence occurrences by municipality must be analyzed. Portugal is divided into 18 districts and 2 autonomous regions. Each of these is subdivided into municipalities. Currently, Portugal has 308 municipalities. The municipality data retrieved from the official statistics website by DGPJ measured the three domestic violence categories for the 308 Portuguese municipalities and for an extra N.E. one, meaning not specified (não especificado in Portuguese). Since it would not be possible to find the explanatory variables values for this special 27 case, this extra municipality was eliminated from the dataset. Furthermore, 12 of the 308 municipalities did not have values for all the categories. Corvo, the smallest island in the Autonomous Region of the Azores, only had data for the DVASA category. The remaining 11 municipalities (Pampilhosa da Serra, Golegã, Ribeira de Pena, Vila de Rei, Barrancos, Vila Viçosa, Penela, Alcoutim, Alfândega da Fé, São Roque do Pico e Aguiar da Beira) were missing data for the DVAM category. The statistical confidentiality principle is stated in Diário da República (the Portuguese official gazette) in Law nº22/2008, the act that legislates on the National Statistical System. This principle, referred to in article 6 of the mentioned law, aims to safeguard citizens' privacy and secure trust in the Statistical System. Therefore, in the retrieved datasets, in order to respect the privacy of the people involved, numbers below 3 are not presented, being symbolized as missing values. Keeping this in mind, the number of missing values was calculated for each category. As one can see from Table 4.3, there were a total of 4,356 missing values among the three categories. The majority of these can be found in the DVAM and Others categories. The missing values for DVASA represent only around 2.3% of the total missing values in the dependent variable dataset. Category Number of Missing Values DVASA 100 DVAM 2,862 Others 1,394 Total 4,356 Table 4.3 - Missing Values by Domestic Violence Category (Municipalities) As mentioned before and seen on Figure 4.1, DVASA is the most prominent category in the total domestic violence occurrences in Portugal. Adding this to the facts that it is also the category with the least missing values (Table 4.3) and that it is the only category measured for all 308 municipalities, one can conclude that this is the best dependent variable for the present study. In order to better understand the missing values and to find the best way to impute them, the difference between the national values for each year and the sum of the values for each municipality in each year was calculated (including the values for N.E.). This can be seen on Table 4.4. One can see that the number of missing values for the municipalities in each year is always very close to the number of occurrences reported on the national level but unreported on a municipal level. In order to address this, the missing values were replaced by the value 1 as it was considered better to keep information about municipalities with low occurrences than to remove them altogether from the study. 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 Missing Values 21 16 6 7 5 7 9 8 4 6 4 7 Difference 27 21 8 7 6 10 7 12 4 10 3 10 Table 4.4 - Missing Values for DVASA in Municipalities and Difference Between National Total and Municipality Total It is important to keep in mind that the retrieved data was measured as an absolute value, which means that it did not consider the differences in the number of inhabitants for each municipality. Keeping the data as it was would have biased the future model, forcing it into thinking that the higher number of domestic violence occurrences in municipalities with the most population was caused by factors other than the number of inhabitants. In order to avoid this problem, the number of DVASA occurrences was 28 standardized according to the resident population in each municipality, as shown below. This way, the dependent variable is now the number of DVASA occurrences per 100 inhabitants. 𝐷𝑉𝐴𝑆𝐴𝑠𝑡𝑎𝑛𝑑𝑎𝑟𝑑𝑖𝑧𝑒𝑑 = 𝐷𝑉𝐴𝑆𝐴𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 (𝑇𝑜𝑡𝑎𝑙 𝑃𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 100 ) The data regarding resident population by municipality used to standardize the dependent variable was retrieved from the Pordata website on the 23rd of April of 2021. The definition of resident population in this case is the group of people who, regardless of being present or absent in a particular accommodation at the time of observation, lived in their usual place of residence for a continuous period of at least 12 months prior to the time of observation, or who arrived at their usual place of residence during the period corresponding to the 12 months preceding the moment of observation, with the intention of remaining there for a minimum period of one year. 4.2. EXPLANATORY VARIABLES When modelling domestic violence one can contemplate two types of variables: risk factors and protective factors. The first ones, as the name suggests, increase the risk of domestic violence, causing a high number of occurrences when very present. Protective factors do the opposite, buffing the risk for domestic violence. The identification of risk factors is very important for the prevention of violence and to guide policies. Both types of factors can be divided into modifiable (for example education) and non-modifiable (for example gender and age) factors. The first ones are the most important when it comes to defining prevention policies, for logical reasons. According to (Ellsberg, Heise, Peña, Agurto, & Winkvist, 2001) individuals belonging to families with more children are more prone to suffer assaults. As a way to include this factor in the present study, the synthetic fertility index (SFI) was considered as an explanatory variable. This index is the average number of children born for each woman in fertile age (between 15 and 49 years). In order for the generation renewal to be assured, the synthetic fertility index must be at 2.1. The data regarding this variable was retrieved from the Pordata website on the 15th of April of 2021 and included data from 2009 to 2019. Another measure for the number of children is the youth dependency index. The data for this variable was retrieved from the Pordata website on the 6th of May of 2021. The youth dependency index is the number of children under 15 years old for every 100 people of working age, that is, between 15 and 64 years old. A value less than 100 means that there are fewer young people than people of working age. This variable had data for all the municipalities (without missing values) for the period between 2009 and 2019. Healthcare workers play an important role in uncovering domestic violence occurrences and supporting the victims. According to (Cann, Withnell, Shakespeare, Doll, & Thomas, 2001), among healthcare workers, women, nurses and mental health workers tend to respond better to domestic violence cases. Keeping this in mind, data regarding the total number of doctors and the number of female doctors for each municipality was retrieved from the Pordata website on the 10th of May of 2021 in order to calculate the percentage of female doctors as below. This data covered the period between 2009 and 2019 and had missing values for some years in Pampilhosa da Serra, Oleiros and Lajes das Flores. 29 𝐹𝑒𝑚𝑎𝑙𝑒𝐷𝑜𝑐𝑡𝑜𝑟𝑠%= 𝐹𝑒𝑚𝑎𝑙𝑒𝐷𝑜𝑐𝑡𝑜𝑟𝑠𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 ∗100 𝑇𝑜𝑡𝑎𝑙𝐷𝑜𝑐𝑡𝑜𝑟𝑠𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 Following the same logic, the number of mental health workers (using psychiatry specialists as a proxy) was retrieved from the Pordata website on the 10th of May of 2021. The percentage of mental health doctors in the total of doctors was calculated as below. Once again, this variable had missing values for some years in Pampilhosa da Serra, Oleiros and Lajes das Flores. 𝑀𝑒𝑛𝑡𝑎𝑙𝐻𝑒𝑎𝑙𝑡ℎ%= 𝑀𝑒𝑛𝑡𝑎𝑙𝐻𝑒𝑎𝑙𝑡ℎ𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 ∗100 𝑇𝑜𝑡𝑎𝑙𝐷𝑜𝑐𝑡𝑜𝑟𝑠𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 Finally, the total number of doctors in each municipality was also used to create a variable that showed the number of doctors per 100 inhabitants of the municipality. This variable had data for the period between 2009 and 2019 and had no missing values. As mentioned in the APAV report regarding male domestic violence victims (APAV - Associação Portuguesa de Apoio à Vítima, 2018), elderly men (65 years or more) tend to be more at risk. This can be seen on Figure 2.1. Taking this into account, the percentage of elderly men in the total population was included in the present study as an explanatory variable. The data regarding the absolute number of men with 65 or more years was retrieved from the Pordata website on the 23rd of April of 2021. This data was then converted to a percentage of the resident population as following: 𝑀𝑒𝑛65%= 𝑀𝑒𝑛65𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 ∗100 𝑇𝑜𝑡𝑎𝑙 𝑃𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 Still focusing on the elderly population, but this time without distinguishing between genders, the elderly dependency index was retrieved from the Pordata website on the 6th of May of 2021. The elderly dependency index is the number of people aged 65 and over for every 100 people of working age, that is, between 15 and 64 years old. A value less than 100 means that there are fewer elderly people than people of working age. This variable had data for all the municipalities (without missing values) for the period between 2009 and 2019. Another way of measuring the level of dependency in a population is to consider the number of Social Security pensioners. A pension is an amount attributed each month to someone in the event of disability, old age, occupational disease or death. Data regarding the number of pensioners for each person who cashes for Social Security was retrieved from the Pordata website on the 11th of May of 2021. This data contemplated the period between 2009 and 2019 and had some missing values for the municipalities of Alenquer, Lagoa (Azores), Lajes das Flores, Santa Cruz das Flores and Corvo. It is also mentioned in another APAV report regarding domestic violence victims in general (APAV - Associação Portuguesa de Apoio à Vítima, 2018) that most of the victims tend to be women with ages comprehended between 26 and 55 years. Keeping this in mind and with the same rationale as for the elderly men variable, the percentage of the population represented by women in these ages was included. The data regarding the absolute number of women between 25 and 54 years was retrieved from the Pordata website on the 27th of April of 2021. The boundaries of the age gap were as close as possible to the ones mentioned in the APAV report. However, they are not exactly the same as this data was not available. The percentage of middle-aged women was calculated as following: 𝑀𝑖𝑑𝑑𝑙𝑒𝐴𝑔𝑒𝑑𝑊𝑜𝑚𝑒𝑛%= 𝑀𝑖𝑑𝑑𝑙𝑒𝐴𝑔𝑒𝑑𝑊𝑜𝑚𝑒𝑛𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 ∗100 𝑇𝑜𝑡𝑎𝑙 𝑃𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 30 The same APAV report mentioned that a high percentage of the victims was married, showing that it might be relevant to include a measure of marriages as an explanatory variable for the present study. However, the number of marriages in a given year does not directly affect the number of domestic violence occurrences in that same year, as the marriage of the victims happens in years before the occurrence. Nevertheless, data regarding the number of marriages national-wide was retrieved from the Pordata website on the 28th of April of 2021 to test for correlations with the number of DVASA occurrences national-wide. When testing for the correlation between absolute values of DVASA occurrences and absolute values for the number of marriages the result was 0.48 which is neither a weak nor a high correlation. However, these values should be standardized according to the population. This standardization was made resulting in DVASA occurrences per 100 inhabitants and the same for the number of marriages. The correlation was then 0.31. Also, the evolution of both variables was plotted to check for common patterns (Figure 4.2) which were mostly not found. One can conclude that it might not be relevant to include this variable in the study. However, it might be interesting to check how this variable behaves when included in a regression and for that purpose data regarding the number of marriages by municipality was also retrieved from the Pordata website on the same date. This data was then standardized to reflect the number of new marriages per 100 inhabitants. Figure 4.2 - Are Marriages and DVASA Occurrences Related? Divorces are yet another controversial variable to include. However, they might be important as, if this variable works as expected according to (Bowlus & Seitz, 2006), it can be a good drive for action. That is, even though divorces are a consequence of domestic violence occurrences, they might be able to explain some of the expected underreporting in domestic violence occurrences. Also, if an increase in divorces is connected to an increase in domestic violence occurrences, the responsible entities can look for a rise in the number of divorces and, in that case, pay closer attention to domestic violence. (WHO - World Health Organization, 2010) mentions divorces as a cause for domestic violence and not only a consequence, as separated or divorced people tend to be more vulnerable, thus becoming more prone to being victims in a following relationship. Data regarding divorces was retrieved from the Pordata website on the 7th of May of 2021 in the form of divorces per 100 marriages. The data included the period from 2009 to 2019 and had some missing values, namely there were no values at all for Odivelas, no data for Castanheira de Pêra in 2013, no data for Barrancos in 2019, no data for Porto Moniz in 2016 and 2018, no data for Corvo in 2013, 2015, 2016, 2018 and 2019. This variable is calculated the following way: 𝐷𝑖𝑣𝑜𝑟𝑐𝑒𝑠 𝑝𝑒𝑟 100 𝑀𝑎𝑟𝑟𝑖𝑎𝑔𝑒𝑠= 𝐷𝑖𝑣𝑜𝑟𝑐𝑒𝑠 𝑖𝑛 𝐶𝑖𝑣𝑖𝑙 𝑌𝑒𝑎𝑟 𝑀𝑎𝑟𝑟𝑖𝑎𝑔𝑒𝑠 𝑖𝑛 𝐶𝑖𝑣𝑖𝑙 𝑌𝑒𝑎𝑟∗100 Another relevant variable according to (Campbell, 2002) is a measure of income, as the poorest strata of the population tend to witness more cases of domestic violence. Considering this, the monthly gain 31 of employees was included in the study. This refers to the amount that the employee receives every month. In addition to the basic remuneration, it includes other remuneration paid by the employer, such as overtime, holiday pay or premiums. It is calculated as a gross amount (before deducting any discounts). This data was retrieved from the Pordata website on the 26th of April of 2021, and it contemplates the period between 2009 and 2018. There was no information for this variable when it comes to all the 19 municipalities in the Autonomous Region of the Azores for the period between 2010 and 2013, which causes a total of 76 missing values. Using the same dataset used for the monthly gain of employees, a measure of the wage gap between men and women was calculated. The data from Pordata, retrieved on the 26th of April of 2021, included the average monthly gain for all employees in a municipality as well as the average monthly gain for women only and for men. Once again, this data refers to the period between 2009 and 2018 and has no values for the municipalities in Azores for the period between 2010 and 2013. The variable here referred to as wage gap is the percentage of the men’s monthly gain that women receive on average and was calculated as following: 𝑊𝑎𝑔𝑒𝐺𝑎𝑝= 𝑀𝑜𝑛𝑡ℎ𝑙𝑦𝐺𝑎𝑖𝑛𝑤𝑜𝑚𝑒𝑛 ∗100 𝑀𝑜𝑛𝑡ℎ𝑙𝑦𝐺𝑎𝑖𝑛𝑚𝑒𝑛 According to (Anderberg, Rainer, Wadsworth, & Wilson, 2015) unemployment also influences domestic violence occurrences. Since the unemployment rate by gender was only available by regions and not municipalities, the number of people enrolled in employment and vocational training centers was used as a proxy. The values were calculated from a simple arithmetic average of the unemployed registered monthly in the employment and vocational training centers, so they are not always whole numbers. This data was retrieved from the Pordata website on the 29th of April of 2021 and had values for the period between 2009 and 2019. To test different possibilities, three variables were created from this data – female unemployment, male unemployment and total unemployment. All of them came in absolute values and had to be standardized by the number of inhabitants in the municipality. This standardization was done in the same way as the standardization of the dependent variable, resulting in the number of people enrolled in employment and vocational training centers by 100 inhabitants. It is also important to notice that there were no values regarding unemployment for the Autonomous Regions of the Azores (19 municipalities) and Madeira (11 municipalities), making it a total of 30 municipalities with no information. Education may also play an important role in explaining the evolution of domestic violence occurrences. According to (Bowlus & Seitz, 2006), victims of violence tend to have lower levels of education and so do the perpetrators. Data regarding the gross enrolment rate (GER) was retrieved from the DGEEC – Direção-Geral de Estatísticas da Educação e Ciência – on the 23rd of June of 2021. The data was available for the period between 2003 and 2019 and did not have values for the municipalities in neither Azores nor Madeira. This indicator is calculated by DGEEC based on DGEEC enrollment data and INE (Instituto Nacional de Estatística) resident population data. It is calculated the following way: 𝐺𝐸𝑅= 𝑆𝑡𝑢𝑑𝑒𝑛𝑡𝑠 𝐸𝑛𝑟𝑜𝑙𝑙𝑒𝑑 𝑖𝑛 𝐻𝑖𝑔ℎ 𝑆𝑐ℎ𝑜𝑜𝑙 𝑅𝑒𝑠𝑖𝑑𝑒𝑛𝑡 𝑃𝑜𝑝𝑢𝑙𝑎𝑡𝑖𝑜𝑛 𝑊𝑖𝑡ℎ 𝑁𝑜𝑟𝑚𝑎𝑙 𝐴𝑔𝑒 𝑓𝑜𝑟 𝐴𝑡𝑡𝑒𝑛𝑑𝑖𝑛𝑔 𝐻𝑖𝑔ℎ 𝑆𝑐ℎ𝑜𝑜𝑙∗100 Regarding the GER, three variables were included: the total GER and GER by gender. 32 5. METHODOLOGY The second step of this study, right after the data collection, was the data treatment and exploration, followed by modelling. Some of the data treatment was already described in the previous chapters, but the remaining part will be described in the present chapter. All calculations and plots were made using Python. The code used in the scope of this study can be found on a GitHub repository 11 . 5.1. SERIES BREAKS A lot of the explanatory variables had breaks caused by changes in the standards for defining and observing the indicator over time. According to the OECD (Organization for Economic Co-operation and Development) Glossary of Statistical Terms, “the specific causes of breaks in a statistical time series include changes in: classifications used, definitions of the variable, coverage, etc.”. The variables GER, GER_Women and GER_Men did not have series breaks. Fertility, Youth_Dependency, Female_Doctors, Mental_Health, Men65, Elderly_Dependency, SS_Pensions, Middle_Aged_Women, Monthly_Gain, Wage_Gap, Unemployment_Total, Unemployment_Female, Unemployment_Male and Total_Doctors all had four breaks, all in 2013. Those were in Lisbon, Loures, Santarém and Golegã. The breaks in Santarém and Golegã are caused by the fact that the parish of Pombalinho was considered, from 2013 on, a parish belonging to the municipality of Golegã, no longer being a part of Santarém. Pombalinho is a small parish, with 7.7km2 and 448 inhabitants, making this break neglectable. It was also a change in the parish configuration that caused the breaks for Lisbon and Loures. In 2013 a new parish called Parque das Nações was created, which included areas from both the municipalities of Loures and Lisbon. This parish is, from 2013 on, part of the municipality of Lisbon and it has 5.44km2 and 21,025 inhabitants. Since only 34.2% of this area and 23.7% of this population belonged to Loures before the change, the effect of this break is neglectable. The variable Marriages also had the four breaks for 2013 described in the last paragraph. However, it also had a break for each municipality in 2010 due to the fact that as of this year (inclusive), with the implementation of Law 9/2010 on the 31st of May, civil marriage between persons of the same gender became allowed. This last break cannot be considered neglectable and, so, data for the year of 2009 had to be removed from the variable in order to cancel the effects of the break. Once again, the variable Divorces had the four breaks in 2013 related to the redistribution of parishes. However, it also had a break for all municipalities for the year of 2010 for the same reason the variable Marriages had a break in 2010 (Law 9/2010). From 2010 on (2010 included), divorces were allowed for persons of the same gender. Once again, this break cannot be considered neglectable and, so, values for the year of 2009 had to be removed for the sake of the coherence of the variable. 5.2. CORRELATIONS It is important to know what the relation between variables within the dataset is. For this purpose, the correlation matrix was calculated using the Pearson Correlation Coefficient. The Pearson correlation is 11 Text is underlined as it represents a link to https://github.com/anastaubyn/MasterThesis 39 Figure 5.3 - Evolution of GER, GER_Men and GER_Women in Barrancos The second scenario, sampling problems, is not applied to the present study, as it is studying a population rather than a sample. Outliers caused by sampling mistakes might occur when data is collected about an individual who does not belong to the target of the study. In this case, the outlier is also considered an error and should be removed from the dataset. Finally, the last scenario is natural variation. These are the outliers that might be important for the study. Data distributions are centered around some point and spread from that point. This means that extreme values might occur but have a lower probability of happening. Even though these data points are unusual, they represent a natural part of the distribution and might add value to the dataset. For example, extreme values in one explanatory variable might justify extreme values in the dependent variable. For this reason, outliers apparently caused by natural variation were not removed from the dataset. Besides the usual summary statistics, histograms and boxplots are common ways of exploring the presence of outliers. The histograms for both the dependent variable and the explanatory ones were presented in Figure 5.2. and Annex II before. However, for a better understanding of the distributions of all variables, a grid of boxplots is presented in Annex III. A boxplot is a summary of the data distribution where the sides of the central rectangle represent the first and third quartiles, the line in the middle of the rectangle represents the median and the “whiskers” represent a measure of 1.5 times the inter-quartile range. Theoretically, all data points that fall outside the “whiskers” of the plot are considered outliers. One can see by the analysis of the boxplots that most of the supposed outliers fall very close to the distribution, as they are not far away from the end of the “whisker”. However, some variables have more extreme values: Divorces, GER_Women, Marriages, Mental_Health and Monthly_Gain. These extreme values may be important in justifying the variance in DVASA. 5.6. STATIONARITY When dealing with time series data (panel data has multiple time series encapsulated inside of it) one can find time-dependent structures such as trend or seasonality. When these structures are present a time series is considered to be non-stationary, as the summary statistics do not remain fixed for all time periods. This causes variations in data that are caused by natural evolution of the numbers but that the model may try to capture anyway, adding bias to the results. Time series that are free of timedependent structures are considered stationary. If the time series is stationary, the covariance 40 between two values of the series depends only on the amount of time separating those values and the summary statistics are fixed independently of the time period. This means that a stationary time series verifies the following conditions: (a) E(yt) = μ; (b) Var(yt) = σ2; (c) Cov(yt, yt+s) = γs. Since the present study is dealing with annual data, seasonality is off the table, as it is mostly found on monthly data, for example. However, trend may still be a problem. The two most common ways to detect the presence of time-dependent structures are visualizing plots of the variables of interest or using a Dickey-Fuller test. Changes in the variables over periods of time are important for visualizing whether a time series is stationary or not. If we have a variable y that is measured over some time periods t (yt), the difference (Δyt) given by yt – yt-1 is the change in the value of variable y from period t-1 to period t and is called the first difference. The plots in Annex IV show the average (yearly mean values considering all municipalities) time series for all explanatory variables plus the dependent one side by side with the time series of the changes for the same variables. Since the mean and variance of a stationary time series are constant, its changes or first differences must fluctuate around a constant value, which does not seem to be the case for the plotted variables. However, as this is a result based on average evolution and changes it may not be the most reliable one. Keeping this in mind, it is necessary to resort to a more valid method, the Dickey-Fuller test. The Dickey-Fuller test is based on the first order autoregressive model in which there are no external explanatory variables and the dependent variable is used as an explanatory variable for itself. This means that the dependent variable is related to past values of itself which means that each value of the variable contains part of the last period’s value plus an error term. If the coefficient of the last period’s value is less than one it means that the variable is stationary. This can be demonstrated as follows: 𝑦𝑡= 𝛽0+ 𝜌𝑦𝑡−1 + 𝑢𝑡 ⇔𝑦𝑡− 𝑦𝑡−1 = 𝛽0+(𝜌−1)𝑦𝑡−1 + 𝑢𝑡 One can see from the equation above that when ρ is less than one, the effects of yt-1 on yt will decrease over time to the point in which a prediction for a period far from t will be unrelated to the value of yt. However, if ρ is equal or more than one, the effects of yt-1 on the values of y for any period will never be annulated, which means that there are time-dependent structures and that the values of y rely heavily on the values of y for past periods. Keeping this in mind, the hypotheses for a Dickey-Fuller test are as following: 𝐻0: 𝜌≥1 𝑣𝑠. 𝐻1: 𝜌<1 It is important to note that the null hypothesis is that of the series not being stationary, meaning that if we fail to reject it, we are assuming non-stationarity. This causes a slight difference in the test statistic as non-stationary series have different properties altering the distribution of the usual t-statistic. To recognize this fact the statistic is called τ (tau) and its values are compared to specifically generated critical values. A problem may arise when performing a simple Dickey-Fuller test as the error term may be autocorrelated. To avoid this, we should add as many lagged differences as needed. This number of 41 differences is determined by examining the autocorrelation function (ACF) of the residuals. This variant of the Dickey-Fuller test is called the Augmented Dickey-Fuller test and its equation is as following: 𝑦𝑡− 𝑦𝑡−1 = 𝛽0+(𝜌−1)𝑦𝑡−1 +∑𝛽𝑠(𝑦𝑡−𝑠 −𝑦𝑡−𝑠−1) 𝑚 𝑠=1 +𝑢𝑡 The dataset that is the object of this study contains 5,560 time series (278 municipalities times 20 variables). The Augmented Dickey-Fuller test was applied to all of these series and it was possible to conclude that for 1,238 of them the null hypothesis was rejected with a significance level of 5%. If the significance level were to be pushed to 10%, the number of stationary series would be 1,502. Using either significance level it becomes clear that most of the series in the dataset are non-stationary. One way to avoid the problems caused by this condition is by removing the trend from the series. However, fixed effect models contemplate the possibility that the intercept may change over different time periods which can also be a solution. Finally, it is also possible to use differences as the variables. Nevertheless, since there is a much higher prevalence of the cross-sectional component in the present dataset than of the time series one this is not a severe problem. The present study is also working with short time series, making the power of Dickey-Fuller tests dubious, as it becomes hard to reject the null hypothesis even if it is false. One of the main problems of non-stationary series is that there is a danger of obtaining apparently significant results from unrelated data, which is called a spurious regression. However, the higher prevalence of the cross-sectional component avoids this problem. 42 5.7. MODEL ESTIMATION Before estimating any models, it is important to understand how the explanatory variables and the dependent variable are related in terms of the functional form to be used. This can be analyzed visually with the aid of scatter plots. Keeping this in mind, a set of scatter plots showing the joint distribution of each explanatory variable with DVASA was plotted, as shown in Figure 5.4 below. Figure 5.4 - Joint Distributions of Explanatory Variables with DVASA From the analysis of the plots above, one can conclude that it might help to include some of the variables in different functional forms. Some explanatory variables, such as the ones related to unemployment and the ones related to education, exhibit a slowly increasing pattern, in a curvilinear form. This means that for higher values of the explanatory variable, the values of the dependent variable increase slowly. This type of relationship can be summed up through a log-log model (with the respective parameter being greater than one) or through a log-linear model (with the respective parameter being greater than 0). Either way, it might be beneficial to turn the dependent variable into a logarithmic form. 43 5.7.1. Model 1 (CC – Constant Coefficients) To estimate the first model a constant coefficients approach was used. All 19 explanatory variables were used for this model in a simple linear form. From now on this model will be called Model 1. The theoretical model is as follows: 𝐷𝑉𝐴𝑆𝐴𝑖𝑡 = 𝛽0+𝛽1𝐹𝑒𝑟𝑡𝑖𝑙𝑖𝑡𝑦𝑖𝑡 +𝛽2𝑀𝑒𝑛65𝑖𝑡 +𝛽3𝑀𝑜𝑛𝑡ℎ𝑙𝑦_𝐺𝑎𝑖𝑛𝑖𝑡 +𝛽4𝑊𝑎𝑔𝑒_𝐺𝑎𝑝𝑖𝑡 +𝛽5𝑀𝑖𝑑𝑑𝑙𝑒_𝐴𝑔𝑒𝑑_𝑊𝑜𝑚𝑒𝑛𝑖𝑡 +𝛽6𝑈𝑛𝑒𝑚𝑝𝑙𝑜𝑦𝑚𝑒𝑛𝑡_𝑇𝑜𝑡𝑎𝑙𝑖𝑡 +𝛽7𝑈𝑛𝑒𝑚𝑝𝑙𝑜𝑦𝑚𝑒𝑛𝑡_𝑀𝑎𝑙𝑒𝑖𝑡 +𝛽8𝑈𝑛𝑒𝑚𝑝𝑙𝑜𝑦𝑚𝑒𝑛𝑡_𝐹𝑒𝑚𝑎𝑙𝑒𝑖𝑡 + 𝛽9𝑀𝑎𝑟𝑟𝑖𝑎𝑔𝑒𝑠𝑖𝑡 + 𝛽10𝐸𝑙𝑑𝑒𝑟𝑙𝑦_𝐷𝑒𝑝𝑒𝑛𝑑𝑒𝑛𝑐𝑦𝑖𝑡 + 𝛽11𝑌𝑜𝑢𝑡ℎ_𝐷𝑒𝑝𝑒𝑛𝑑𝑒𝑛𝑐𝑦𝑖𝑡 +𝛽12𝐹𝑒𝑚𝑎𝑙𝑒_𝐷𝑜𝑐𝑡𝑜𝑟𝑠𝑖𝑡 +𝛽13𝑇𝑜𝑡𝑎𝑙_𝐷𝑜𝑐𝑡𝑜𝑟𝑠𝑖𝑡 +𝛽14𝑀𝑒𝑛𝑡𝑎𝑙_𝐻𝑒𝑎𝑙𝑡ℎ𝑖𝑡 + 𝛽15𝑆𝑆_𝑃𝑒𝑛𝑠𝑖𝑜𝑛𝑠𝑖𝑡 +𝛽16𝐺𝐸𝑅𝑖𝑡 +𝛽17𝐺𝐸𝑅_𝑀𝑒𝑛𝑖𝑡 +𝛽18𝐺𝐸𝑅_𝑊𝑜𝑚𝑒𝑛𝑖𝑡 +𝛽19𝐷𝑖𝑣𝑜𝑟𝑐𝑒𝑠𝑖𝑡 +𝑢𝑖𝑡 Model 1 was estimated using cluster-robust standard errors that allow for a relaxation of the assumption regarding the covariance of the error term. This means that the error for observations of the same municipality may be correlated but the covariance of error terms for different municipalities must be 0. The R-squared for this model was 0.1342 meaning that 13.42% of the variability in DVASA is represented by the set of explanatory variables chosen for this model. The result of the overall significance test for Model 1 was very positive, with a p-value of 0. Table 5.3 sums up the parameter estimates and their significance for this model. Parameter Standard Error T-stat P-value Constant 0.4115 0.1620 2.5409 0.0111 Fertility 0.0282 0.0105 2.6893 0.0072 Men65 -0.0103 0.0065 -1.5837 0.1134 Monthly_Gain 0.0000413 0.000027 1.5278 0.1267 Wage_Gap -0.0003 0.0004 -0.6158 0.5381 Middle_Aged_Women -0.0088 0.0044 -1.9877 0.0469 Unemployment_Total 4.5520 2.2377 2.0342 0.0420 Unemployment_Male -4.5526 2.2369 -2.0352 0.0419 Unemployment_Female -4.5356 2.2386 -2.0261 0.0428 Marriages 0.0368 0.0186 1.9731 0.0486 Elderly_Dependency 0.0010 0.0012 0.8918 0.3726 Youth_Dependency -0.0031 0.0013 -2.3648 0.0181 Female_Doctors -0.0000685 0.0002 -0.3551 0.7226 Total_Doctors 0.0393 0.0178 2.2053 0.0275 Mental_Health -0.0006 0.0011 -0.6031 0.5465 SS_Pensions -0.0460 0.0225 -2.0474 0.0407 GER -0.0008 0.0005 -1.7929 0.0731 GER_Men 0.0006 0.0003 2.0282 0.0426 GER_Women 0.0004 0.0002 1.5352 0.1248 Divorces 0.0002 0.0000647 2.6329 0.0085 Table 5.3 - Parameter Estimates for Model 1 The last two columns of Table 5.3 show the results of the individual significance test for each parameter. P-values that exceed the significance level of 5% are highlighted in red. In this test, if the null hypothesis is rejected, one can say that there is statistical evidence that the parameter is not 0, 44 which means that the corresponding variable has an actual impact on the dependent variable. Considering a level of significance of 5%, one can see from Table 5.3 that there is no statistical evidence of some of the parameters being relevant to the model. It is the case of the parameters corresponding to GER_Women, GER, Mental_Health, Female_Doctors, Elderly_Dependency, Wage_Gap, Monthly_Gain and Men65. If one were to consider a significance level of 10%, only GER would be removed from the previous list. In future models this conclusion will be considered. In order to test for heteroskedasticity in Model 1, a graphical approach was used first, plotting the residuals against the predicted values. If any pattern is detected in the plot (increasing or decreasing variance), there is heteroskedasticity in the model. The plot is shown in Figure 5.5. Figure 5.5 - Heteroskedasticity Test for Model 1 According to Figure 5.5, there does not seem to be any evidence of heteroskedasticity in the model. However, to guarantee this result, a formal test must be made. The Breusch-Pagan test is meant to detect heteroskedasticity. It takes the squares of the residuals created by the model and makes a regression on it with the same dependent variables as the original model. Then, an overall significance test is performed to see if all the coefficients can be 0. The results for this test show a really low pvalue (1.27825047691774e-16), indicating that there is statistical evidence that at least one of the coefficients is not equal to 0, which means that there is heteroskedasticity in the model. Consequently, the estimators are not efficient. However, since cluster-robust standard errors were used, the results of statistical tests can be evaluated. No other tests on heteroskedasticity were performed for Model 1 due to the extremely low p-value for the Breusch-Pagan test. In order to test for autocorrelation among residuals, the Durbin-Watson test was applied. The DurbinWatson statistic is calculated using the residual sum of squares and the sum of squares of differences between consecutive residuals. It takes a value between 0 and 4 where the middle value (2) indicates no autocorrelation. Values between 0 and 2 indicate positive autocorrelation and values between 2 and 4 indicate negative autocorrelation. The statistic for Model 1 was 1.8941, which indicates a small level of positive autocorrelation. Considering the results of the two tests, it might be a better option to use a fixed-effects or a randomeffects model. Nevertheless, the plots of the residuals versus the explanatory variables for which the individual significance test showed statistical evidence of importance were plotted to check for wrong functional forms. These plots can be seen on Figure 5.6. 45 Figure 5.6 - Residuals Against Explanatory Variables (Model 1) After analyzing the plots in Figure 5.6 there does not seem to be any evidence for the existence of a better functional form for any of the variables. To allow better future model comparison, the adjusted R-squared was calculated for Model 1, using the following formula where N is the total number of observations and k is the total number of explanatory variables: 𝑅 2=1−(1−𝑅2)(𝑁−1) 𝑁−𝑘−1 Using the formula above, the adjusted R-squared for Model 1 is 0.1288. 5.7.2. Model 2 (CC) Taking into account the results from Model 1, and still using a constant coefficients approach, a new model was estimated keeping all variables in a simple linear form and removing variables that do not seem to be relevant. This model will from now on be referred to as Model 2. Model 2 was also estimated using cluster-robust standard errors. Its R-squared was 0.1205, lower than Model 1. This is expected since 8 variables were removed. The p-value for the overall significance test was still 0, meaning that there is really strong statistical evidence of the importance of the set of explanatory variables used in justifying the values for DVASA. Table 5.4 sums up the parameter estimates and their significancy for this model. From Table 5.4 below, one can see that, considering a level of significance of 5%, the null hypothesis for the individual significance test cannot be rejected for some parameters. These are the ones corresponding to Middle_Aged_Women, Marriages, Youth_Dependency and GER_Men. However, if a level of significance of 10% was to be considered, the only parameters for which the null hypothesis 46 could not be rejected would be the ones corresponding to Middle_Aged_Women and Marriages. These values are highlighted in red on the previous table. Parameter Standard Error T-stat P-value Constant 0.2479 0.0734 3.3766 0.0007 Fertility 0.0307 0.0108 2.8491 0.0044 Middle_Aged_Women -0.0042 0.0028 -1.5342 0.1251 Unemployment_Total 4.8803 2.2706 2.1493 0.0317 Unemployment_Male -4.8805 2.2698 -2.1502 0.0316 Unemployment_Female -4.8643 2.2716 -2.1414 0.0323 Marriages 0.0306 0.0194 1.5812 0.1139 Youth_Dependency -0.0021 0.0011 -1.8705 0.0615 Total_Doctors 0.0451 0.0170 2.6548 0.0080 SS_Pensions -0.0536 0.0184 -2.9112 0.0036 GER_Men 0.0001 0.000058 1.7544 0.0795 Divorces 0.0002 0.000067 2.6164 0.0089 Table 5.4 - Parameter Estimates for Model 2 In a similar way to what was done for Model 1, a graphical approach was used to test Model 2 for heteroskedasticity. The plot is shown in Figure 5.7 and, once again, no pattern of increasing or decreasing variance was detected. Figure 5.7 - Heteroskedasticity Test for Model 2 To have a more formal result regarding the presence of heteroskedasticity in Model 2, a Breusch-Pagan test was applied. The results for this test for this test show a really low p-value for the F-test (1.6751524272085239e-16), indicating that there is statistical evidence that at least one of the coefficients is not equal to 0, which means that there is heteroskedasticity present in the model. Consequently, the estimators are not efficient. However, since cluster-robust standard errors were used, the results of statistical tests can be evaluated. No other tests on heteroskedasticity were performed for Model 2 due to the extremely low p-value for the Breusch-Pagan test. Similarly to what happened regarding Model 1, to test for autocorrelation among residuals the DurbinWatson test was applied. The statistic of this test for Model 2 was 1.9072, which indicates a very small level of positive autocorrelation. Since this is a close enough value to 2, one can assume that there is no significant autocorrelation. 47 The value of the adjusted R-squared for Model 2 is 0.1173, lower than the value for Model 1. 5.7.3. Model 3 (CC) According to the conclusions taken from Figure 5.4, some transformations were applied to the variables, namely converting DVASA into a natural logarithmic form. A model with all variables used for Model 1 was estimated, but this time the variables were included as log-linear relationships. This model will from now on be addressed to as Model 3. Once again, Model 3 used a constant coefficients approach with cluster-robust standard errors. The R-squared for Model 3 was 0.1396. This means that 13.96% of the variability in ln(DVASA) is represented by this model. We can also get information on the between and within R-squared, which, similarly to the within and between standard deviations, represent the part of the differences between municipalities that is explained by the model and the part of the differences within municipalities that is explained by the model. For Model 3 the between R-squared was 0.2711 and the within R-squared was 0.0320. From these numbers, one can easily understand that a much larger part of the variance being addressed by the model is justified by differences between municipalities. The result of the overall significance test for Model 3 was very positive, with a p-value of 0. However, according to the T-tests for individual significance, considering a significance level of 5%, there was no statistical evidence of the relevance of many variables, including the constant. The results for Model 3 can be seen on Table 5.5 below, where the p-values greater than the intended significance level (5%) are highlighted in light red. Parameter Standard Error T-stat P-value Constant -0.2051 0.9792 -0.2094 0.8341 Fertility 0.1422 0.0682 2.0849 0.0372 Men65 -0.0653 0.0441 -1.4818 0.1385 Monthly_Gain 0.0002 0.0002 1.5164 0.1295 Wage_Gap -0.0020 0.0028 -0.7318 0.4643 Middle_Aged_Women -0.0539 0.0275 -1.9602 0.0501 Unemployment_Total 42.762 18.404 2.3235 0.0202 Unemployment_Male -42.749 18.398 -2.3236 0.0202 Unemployment_Female -42.685 18.410 -2.3185 0.0205 Marriages 0.1801 0.1138 1.5828 0.1136 Elderly_Dependency 0.0055 0.0080 0.6889 0.4909 Youth_Dependency -0.0191 0.0080 -2.3778 0.0175 Female_Doctors 0.0006 0.0013 0.4505 0.6524 Total_Doctors 0.1949 0.0875 2.2265 0.0261 Mental_Health -0.0028 0.0068 -0.4125 0.6800 SS_Pensions -0.3451 0.1494 -2.3097 0.0210 GER -0.0037 0.0027 -1.3367 0.1814 GER_Men 0.0029 0.0015 1.8992 0.0576 GER_Women 0.0013 0.0015 0.8358 0.4033 Divorces 0.0010 0.0004 2.4794 0.0132 Table 5.5 - Parameter Estimates for Model 3 48 The variables that were not statistically significant were analyzed one by one, in order to check if there were possible improvements to the functional form by converting the variables in a natural logarithmic form as well, creating log-log relationships. This was done by analyzing the plots in Figure 5.4. Men65 is included in Model 3 as a log-linear relationship in which the correspondent parameter is smaller than 0. However, it seems to resemble more the shape of a log-log relationship with the correspondent parameter being equal to -1 or below it. Thus, Men65 will be included in the next model as a log-log relationship. Monthly_Gain can also benefit from a log-log relationship with the corresponding parameter being smaller than 0. Wage_Gap does not exhibit any clear pattern and has too high of a pvalue, meaning that it will be excluded from the next model. Middle_Aged_Women was included in Model 3 as a log-linear relationship with the corresponding coefficient lower than 0. However, it seems to be best described as a log-log relationship with the correspondent coefficient greater than 0. It will then be included in the next model as a log-log relationship. Marriages does not exhibit a clear pattern and will be excluded from the next model. Elderly_Dependency also seems to benefit from a log-log relationship and will be included in this way in the next model. Both Female_Doctors and Mental_Health do not exhibit clear patterns and have very high p-values, being excluded from the next model. GER and GER_Men will also be included as log-log relationships in the next model. The adjusted R-squared for Model 3 was 0.1342. 5.7.4. Model 4 (CC) Model 4 was, once again, estimated using a constant coefficients approach with cluster-robust standard errors. The variables used were described in the final paragraph of Chapter 5.7.3. This model presented a R-squared of 0.1395, with a within R-squared of 0.0336 and a between R-squared of 0.2687. Once again, differences between municipalities are much better explained than differences within municipalities. The p-value of the overall significance F-test was 0. The parameter estimates are summed up in Table 5.6 below. The dependent variable for this model is the natural logarithm of DVASA. Parameter Standard Error T-stat P-value Constant -2.2822 2.0980 -1.0878 0.2768 Fertility 0.1445 0.0675 2.1413 0.0323 Ln_Men65 -0.6080 0.4721 -1.2878 0.1979 Ln_Monthly_Gain 0.3882 0.1351 2.8737 0.0041 Ln_Middle_Aged_Women -0.6790 0.5012 -1.3547 0.1756 Unemployment_Total 42.802 18.588 2.3026 0.0214 Unemployment_Male -42.794 18.581 -2.3031 0.0213 Unemployment_Female -42.718 18.596 -2.2972 0.0217 Ln_Elderly_Dependency 0.3086 0.3816 0.80088 0.4187 Youth_Dependency -0.00168 0.0081 -2.0704 0.0385 Total_Doctors 0.1848 0.0805 2.2963 0.0217 SS_Pensions -0.3240 0.1416 -2.2876 0.0222 Ln_GER 0.0064 0.0659 0.0974 0.9224 Ln_GER_Men 0.0792 0.0538 1.4707 0.1415 GER_Women -0.0003 0.0004 -0.775 0.4384 55 6. RESULTS AND DISCUSSION Nine models were created in the scope of this study. Four of them used the constant coefficients approach and the remaining five used the fixed effects approach. The objective was to build a model that is substantiated by statistical evidence. For this purpose, results of overall significance tests and individual significance tests must by analyzed. The R2 or coefficient of determination of a regression is a measure of the proportion of variability in the dependent variable that is justified by changes in the explanatory variables. However, the R2 may not always be the best measure to compare models as it always becomes larger when adding new variables, even if the variables added are completely irrelevant. Keeping this in mind, the R2 can only be used to compare models with the exact same number of explanatory variables. When the models have a different number of explanatory variables, the adjusted R2 can be used as an alternative. Table 6.1 below shows a summary of the measures calculated for the models created to facilitate comparison between them. In this table, one can see the dependent variable used for each model, the number of explanatory variables included in each model, the values for the overall R-squared, between R-squared, within R-squared and adjusted R-squared for each model, the number of the variables included as explanatory variables that had a p-value lower than 0.1 and were, therefore, considered relevant based on a 10% significance level, the p-value for the overall significance F-test, the p-value for the Breusch-Pagan test and the value for the Durbin-Watson statistic. Model 1 Model 2 Model 3 Model 4 Model 5 Model 6 Model 7 Model 8 Model 9 Dependent Variable DVASA DVASA Ln(DVASA) Ln(DVASA) DVASA DVASA DVASA DVASA DVASA Nº Explanatory Variables 19 11 19 15 29 24 19 17 15 Nº Observations 3,058 3,058 3,058 3,058 3,058 3,058 3,058 3,058 3,058 R-squared 0.1342 0.1205 0.1396 0.1395 0.1664 0.1657 0.1607 0.1577 0.1566 Between Rsquared 0.2592 0.2460 0.2711 0.2687 0.2957 0.2941 0.2918 0.2882 0.2877 Within Rsquared 0.0181 0.0040 0.0320 0.0336 0.0462 0.0463 0.0389 0.0365 0.0349 Adjusted Rsquared 0.1288 0.1173 0.1342 0.1353 0.1584 0.1591 0.1556 0.1530 0.1524 Nº Relevant Variables (10%) 10 9 10 9 19 16 19 14 15 Overall Significance 0 0 0 0 0 0 0 0 0 BreuschPagan ≈0 ≈0 ≈0 ≈0 ≈0 ≈0 ≈0 ≈0 0.00196 DurbinWatson 1.8941 1.9072 1.9151 1.9068 1.9245 1.9259 1.9443 1.9396 1.9416 Table 6.1 - Model Comparison Model 1, the first one created, included all retrieved possible explanatory variables in a simple linear form. It showed promising results concerning the value for the between R-squared. However, almost 56 half of the variables could not be considered relevant and the value for the within R-squared was extremely low, exposing from the start the problem of low variance of the data within municipalities. To try and achieve a more statistically relevant model, Model 2 was created, removing eight of the non-relevant variables exposed in Model 1 while keeping all relationships in a simple linear form. The remaining two non-relevant variables from Model 1 were kept in Model 2 to see if their behavior changed when included in a more restricted set of explanatory variables. However, their corresponding individual significance tests continued to show them non-relevant. As the value for the R-squared was very low, Model 3 was created using different functional forms (log-linear) to check for improvements in explainability. This model showed promising results but many variables had high p-values for their individual significance test, causing the model to be not very reliable. A further attempt with the natural logarithm of DVASA as the dependent variable was made in Model 4, but this time including some of the variables in a log-log form and removing the ones that did not show apparent patterns when plotted against DVASA. Once again, this model showed promising results but many variables had high p-values for their individual significance tests. As all the previous models were showing clear evidence of heteroskedasticity, a change in the model approach was considered, resulting in Model 5, estimated using fixed effects for the time component. In this model, differences between constants for different years were included using dummy variables for all years except 2009, to avoid multicollinearity. From Model 5 to Model 9, different combinations of variables were tried in order to optimize interpretability and reliance of the model. All of the estimated models ended up showing evidence of the presence of heteroskedasticity. This means that the estimators are not efficient and that, if cluster-robust standard errors had not been used in all estimations, statistical inference would not be valid. As already mentioned, valid statistical evidence is the main criteria when it comes to choosing the final model. Thus, all variables of that model must be significant and the model itself must pass the overall significance test. The only two models that meet these criteria are Model 7 and Model 9. Since the number of variables of these two models is different, the R-squared must not be used as a comparison measure. Instead, the comparison measure to be used is the adjusted R-squared. However, even if Model 7 shows a higher value for this measure, Model 9 has better interpretability, making it a serious contestant for final model. A choice was made to value interpretability in this case as the difference between R-squared values is not much significant. That being said, Model 9 was chosen as the final model. 57 7. CONCLUSIONS As stated in the Introduction, the present dissertation proposed to answer the following questions: ▪ How did the number of domestic violence occurrences in Portugal evolve between 2009 and 2019? ▪ How well can panel data regression explain this evolution? ▪ What are the main causes of domestic violence? ▪ How does each explanatory variable affect the number of domestic violence occurrences? In order to achieve the proposed goal, several models were estimated using the number of domestic violence occurrences against spouse or analogous or variations of this as the dependent variable. First of all, it is possible to conclude that explaining the evolution of domestic violence occurrences is a difficult process as variables on the individual level cannot be included to explain occurrences in a municipality. Due to this, the models estimated ended up justifying approximately 15 or 16% of the total variance in the dependent variable. The model selected as final can be written as following: 𝐷𝑉𝐴𝑆𝐴𝑖=0.1188+0.0180𝐷2010𝑖+0.0106𝐷2011𝑖+0.0133𝐷2014𝑖+0.0175𝐷2015𝑖+0.0222𝐷2016𝑖 +0.0345𝐷2017𝑖+0.0477𝐷2018𝑖+0.0697𝐷2019𝑖+0.0303𝐹𝑒𝑟𝑡𝑖𝑙𝑖𝑡𝑦𝑖 +0.0002𝐷𝑖𝑣𝑜𝑟𝑐𝑒𝑠𝑖+0.0001𝐺𝐸𝑅𝑖+0.0414𝑀𝑎𝑟𝑟𝑖𝑎𝑔𝑒𝑠𝑖−0.0011𝑊𝑎𝑔𝑒_𝐺𝑎𝑝𝑖 +0.0125𝑈𝑛𝑒𝑚𝑝𝑙𝑜𝑦𝑚𝑒𝑛𝑡_𝑇𝑜𝑡𝑎𝑙𝑖+0.0364𝑇𝑜𝑡𝑎𝑙_𝐷𝑜𝑐𝑡𝑜𝑟𝑠𝑖+𝑒𝑖 This means that in 2009, 2012 and 2013, if all the other variables are set to 0 and holding everything else constant, on average, the number of DVASA occurrences registered by police authorities per 100 inhabitants was 0.1188. In 2010, it is estimated that municipalities had, on average and holding all other factors constant an extra 0,018 DVASA occurrences registered by police authorities when compared to 2009, 2012 and 2013. According to the same rationale, for 2011, 2014, 2015, 2016, 2017, 2018 and 2019, it is expected that municipalities, holding all other factors constant, have an extra 0.0106, 0.0133, 0.0175, 0.0222, 0.0345, 0.0477 and 0.0697, respectively, DVASA occurrences registered by police authorities for 100 inhabitants when compared to 2009, 2012 and 2013. Fertility represents the average number of children born to each woman in fertile age. The coefficient for this variable is 0.0303, which means that, holding all other factors constant, it is expected that an increase of one child per woman in fertile age causes an increase of 0.0303 DVASA occurrences per 100 inhabitants in that same municipality. This corroborates the idea stated in (Ellsberg, Heise, Peña, Agurto, & Winkvist, 2001) that claims that women with more children are more likely to become victims of domestic violence. The variable Divorces represents the number of divorces per 100 marriages. Similarly to Fertility, this variable has also proven to be a risk factor. As the number of divorces increases, the number of DVASA occurrences is also expected to increase, which can be seen in the positive coefficient associated to this variable in the model. One can say that, if the number of divorces by 100 marriages increases by 10, holding all other factors constant, it is expected that the number of DVASA occurrences per 100 inhabitants increases by 0.002. This corroborates the theory stated in (WHO - World Health Organization, 2010), affirming that divorces cannot only be a consequence of domestic violence, but also a cause. People who are separated or divorced tend to be more vulnerable, increasing their probability of becoming victims of domestic violence. 58 GER represents the percentage of resident population with normal age for attending high school that is actually attending high school. The coefficient for this variable being positive was an astonishing find. It goes against all preconceived ideas and all ideas exposed in the literature review. It shows that as the percentage represented in GER increases, it is expected that so does the number of DVASA occurrences per 100 inhabitants. This variable was retrieved due to lack of data regarding the education level of the population by municipality, so it might be interesting to try the same regression using a more appropriate variable, if made available. The variable Marriages represents the number of new marriages per 100 inhabitants and has a positive coefficient associated, which makes it a risk factor. This was expected as most of domestic violence victims are married. However, this was a controversial variable as it should be added in a way that reflects past marriages. Nevertheless, municipalities with more marriages may also be more conservative, which adds a new hypothesis on the table. From the model one can see that it is expected that for each new marriage by 100 inhabitants, the number of DVASA occurrences by 100 inhabitants increases by 0.0414, holding all other factors constant. Wage_Gap represents the percentage of men’s monthly pay that women receive on average. The coefficient for this variable is negative, which meets the theory of exposure reduction, explained in (Aizer, 2010) which states that, as the wage gap decreases, the labor force participation of women increases and, consequently, domestic violence against them declines because women spend less time with violent partners. According to the model built, for each percentual point that women salary gains when faced with men monthly gain, it is expected that the number of DVASA occurrences by 100 inhabitants decreases by 0.0011, holding all other factors constant. Unemployment_Total represents the number of people enrolled in employment and vocational training centers per 100 inhabitants. As expected, the coefficient for this variable is positive, making unemployment a risk factor. According to the model, it is expected that for each 10 additional persons enrolled in employment and vocational training centers, the number of DVASA occurrences per 100 inhabitants increases by 0.0125. Finally, Total_Doctors was added to the dataset to try and fight underreporting, as doctors are the professionals who deal directly with the consequences of domestic violence and can report occurrences to the authorities. This variable represents the number of doctors per 100 inhabitants according to the Doctors’ Professional Order. It has a positive associated coefficient, as expected. According to the model built, it is expected that each new doctor per 100 inhabitants leads to an increase of reported DVASA occurrences by 100 inhabitants of 0.0364, holding all other factors constant. 59 8. LIMITATIONS AND RECOMMENDATIONS FOR FUTURE WORKS Explaining a good proportion of the total variance encapsulated in the dependent variable was revealed to be a hard task. When choosing to study domestic violence occurrences using a social structural theory, it is almost inevitable that explanatory variables are most of the times very static. A lot of macro-level variables tend to evolve very slowly, causing a problem of low variance in short time series. This is what happened with data gathered for each municipality and this is the reason why, for every model created, much higher values were presented for the between R-squared than for the within R-squared. Nevertheless, the main goal of a model like the ones built in the present document is to allow authorities to act on the problem and, for that matter, differences on a spatial level are more important than on a temporal level. Nevertheless, it is my opinion that studies focused on the individual-level should be done to find more prominent causes of domestic violence. Furthermore, it is impossible to know for sure whether the variance in the dependent variable reflects reality. As mentioned before, domestic violence is a crime that commonly takes place in the privacy of a home and, for that reason, many cases may depend on self-report. This means that according to (Ellsberg, Heise, Peña, Agurto, & Winkvist, 2001), the number of occurrences registered may suffer from underreporting. Explanatory variables also have their own issues. Due to unavailability of data, some explanatory variables are not exactly those originally thought for the theoretical model. It is the case of all variables related to unemployment, variables related to education, among others. The variables that ended up being used are as close as possible to the original idea, but they might have led to a not so ideal set of independent variables in the final model. Finally, the only two panel data analysis approaches used were the constant coefficients and the fixed effects. There was an attempt to apply a random effects approach, but it did not show very promising results and, considering the scope of this work, was not pursued. However, I believe that it may be possible to build interesting models on this subject using random effects. Another option would be to try different slopes for different years by creating interactions between the dummy variables already created and other explanatory variables. 60 9. BIBLIOGRAPHY Ackerson, L., Kawachi, I., Barbeau, E., & Subramanian, S. (2008). Effects of Individual and Proximate Educational Context on Intimate Partner Violence: A Population-Based Study of Women in India. American Journal of Public Health, 507-514. Aizer, A. (2010). The Gender Wage Gap and Domestic Violence. American Economic Review, 18471859. Amirthalingam, K. (2005). Women's Rights, International Norms, and Domestic Violence:. Human Rights Quarterly,, 683-708. Anderberg, D., Rainer, H., Wadsworth, J., & Wilson, T. (2015). Unemployment and Domestic Violence: Theory and Evidence. The Economic Journal, 1947-1979. APAV - Associação Portuguesa de Apoio à Vítima. (2018). Homens Vítimas de Violência Doméstica 2013 - 2017. Lisboa: Associação Portuguesa de Apoio à Vítima. APAV - Associação Portuguesa de Apoio à Vítima. (2018). Vítimas de Violência Doméstica 2013 - 2017. Lisboa: Associação Portuguesa de Apoio à Vítima. APAV - Associação Portuguesa de Apoio à Vítima. (2020). Vítimas de Homicído Relatório APAV 2019. Lisboa: Associação Portuguesa de Apoio à Vítima. Barber, C. (2008). Domestic Violence Against Men. Nursing Standard, 35-39. Bowlus, A., & Seitz, S. (2006). Domestic Violence, Employment and Divorce. International Economic Review, 1113-1149. Brasil, E., Alves, F., & Soares, S. (2018). Dados 2017. Almada: OMA - Observatório de Mulheres Assassinadas da UMAR. Brasil, E., Alves, F., & Soares, S. (2019). Dados 2018. Almada: OMA - Observatório de Mulheres Assassinadas da UMAR. Campbell, J. (2002). Health Consequences of Intimate Partner Violence. The Lancet, 1331-1336. Cann, K., Withnell, S., Shakespeare, J., Doll, H., & Thomas, J. (2001). Domestic Violence: A Comparative Survey of Levels of Detection, Knowledge and Attitudes in Healthcare Workers. Public Health, 89-95. Devries, M., Mak, T., Garcia-Moreno, C., Petzold, M., Child, C., Falder, G., . . . Watts, H. (2013). The Global Prevalence of Intimate Partner Violence Against Women. Science, 1527-1528. Ellsberg, M., Heise, L., Peña, R., Agurto, S., & Winkvist, A. (2001). Researching Domestic Violence Against Women: Methodological and Ethical Considerations. Studies in Family Planning. Hill, R. C., Griffiths, W. E., & Lim, G. C. (2012). Principles of Econometrics. John Wiley and Sons. 61 Soares, S., Branco, E., & Alves, F. (2020). Relatório Anual 2019. Almada: OMA - Observatório de Mulheres Assassinadas da UMAR. WHO - World Health Organization. (2010). Preventing Intimate Partner and Sexual Violence Against Women: Taking Action and Generating Evidence. Geneva: World Health Organization. Wooldridge, J. M. (2013). Introductory Econometrics: A Modern Approach. Michigan: South-Western. 62 10. ANNEXES Annex I. Contemporaneous Correlation Matrixes 63 64 Annex II – Histograms for Explanatory Variables