{"id":9875,"date":"2020-09-28T08:00:51","date_gmt":"2020-09-28T08:00:51","guid":{"rendered":"https:\/\/www.experfy.com\/blog\/?p=9875"},"modified":"2023-10-27T13:27:33","modified_gmt":"2023-10-27T13:27:33","slug":"using-regression-with-correlated-data","status":"publish","type":"post","link":"https:\/\/www.experfy.com\/blog\/bigdata-cloud\/using-regression-with-correlated-data\/","title":{"rendered":"Using Regression with Correlated Data"},"content":{"rendered":"\t\t
\n\t\t\t\t\t\t
\n\t\t\t\t\t\t
\n\t\t\t\t\t
\n\t\t\t
\n\t\t\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t\n

Tutorial (including R code) for using Generalized Estimating Equations and Multilevel Models.<\/p>\n\n\n\n

While regression models are easy to run given their short, simple syntax, this accessibility also makes it easy to use regression inappropriately. These models have several key assumptions that need to be met in order for their output to be valid, but your code will typically run whether or not these assumptions have been met.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

\n\t\t\t\t
\n\t\t\t\t\t\t\t
\n\t\t\t
<\/div>\t\t<\/div>\n\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t\n

For linear regression (used with a continuous outcome), these assumptions are as follows:<\/p>\n\n\n\n

    \n
  1. Independence: All observations are independent of each other, residuals are uncorrelated<\/strong><\/li>\n
  2. Linearity: The relationship between X and Y is linear<\/li>\n
  3. Homoscedasticity: Constant variance of residuals at different values of X<\/li>\n
  4. Normality: Data should be normally distributed around the regression line<\/li>\n<\/ol>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
    \n\t\t\t\t
    \n\t\t\t\t\t\t\t\t\t\n

    For logistic regression (used with a binary or ordinal categorical outcome), these assumptions are as follows:<\/p>\n\n\n\n

      \n
    1. Independence: All observations are independent of each other, residuals are uncorrelated<\/strong><\/li>\n
    2. Linearity in the logit: The relationship between X and the logit of Y is linear<\/li>\n
    3. Model is correctly specified, including lack of multicollinearity<\/li>\n<\/ol>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
      \n\t\t\t\t
      \n\t\t\t\t\t\t\t\t\t\n

      In both kinds of simple regression models, independent observations are absolutely necessary to fit a valid model. If your data points are correlated, this assumption of independence is violated. Fortunately, there are still ways to produce a valid regression model with correlated data.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

      \n\t\t\t\t
      \n\t\t\t\t\t

      Correlated Data<\/h2>\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
      \n\t\t\t\t
      \n\t\t\t\t\t\t\t\t\t\n

      Correlation in data occurs primarily through multiple measurements (e.g. two measurements are taken on each participant 1 week apart, and data points within individuals are not independent) or if there is clustering in the data (e.g. a survey is conducted among students attending different schools, and data points from students within a given school are not independent).<\/p>\n\n\n\n

      The result is that that the outcome has been measured on the level of an individual observation, but that there is a\u00a0second level<\/a>\u00a0of either an individual (in the case of multiple time points) or clusters on which individual data points can be correlated. Ignoring this correlation means that standard error cannot be accurately computed, and in most cases will be artificially low.<\/p>\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

      \n\t\t\t\t
      \n\t\t\t\t\t\t\t\t\t\n

      The best way to know if your data is correlated is simply through familiarity with your data and the collection process that produced it. If you know that you have repeated measures from the same individuals or have data on participants who can be grouped into families or schools, you can assume that your data points are probably not independent. You can also\u00a0investigate your data for possible correlation<\/a>\u00a0by calculating the ICC (intraclass correlation coefficient) to determine how correlated data points are within possible groups, or by looking for correlation in your residuals.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

      \n\t\t\t\t
      \n\t\t\t\t\t

      Regression Modeling with Correlated Data<\/h2>\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
      \n\t\t\t\t
      \n\t\t\t\t\t\t\t\t\t\n

      As previously mentioned, simple regression will produce inaccurate standard errors with correlated data and therefore should not be used.<\/p>\n\n\n\n

      Instead, you want to use models that can account for the correlation that is present in your data. If the correlation is due to some grouping variable (e.g. school) or repeated measures over time, then you can choose between Generalized Estimating Equations or Multilevel Models. These modeling techniques can handle either binary or continuous outcome variables, so can be used to replace either logistic or linear regression when the data are correlated.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

      \n\t\t\t\t
      \n\t\t\t\t\t

      Generalized Estimating Equations<\/h3>\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
      \n\t\t\t\t
      \n\t\t\t\t\t\t\t\t\t\n

      Generalized estimating equations (GEE) will give you beta estimates that are the same or similar to those produced by simple regression, but with appropriate standard errors. Generalized estimating equations are particularly useful when you have repeated measures for the same individuals or units. This modeling technique tends to work well when you have many small clusters, which is often the result of having a few measurements on a large number of participants. GEE also allows the user to\u00a0specify one of numerous correlation structures, which can be a useful feature depending on your data.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

      \n\t\t\t\t
      \n\t\t\t\t\t

      Multilevel Modeling<\/h3>\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
      \n\t\t\t\t
      \n\t\t\t\t\t\t\t\t\t\n

      Multilevel modeling<\/a>\u00a0(MLM) also provides appropriate standard errors when data points are not independent. It is typically the best modeling approach when the user is interested in relationships both within and between clustered groups, and is not simply looking to account for the effect of correlation in standard error estimates. MLM has the additional advantage of being able to handle more than two levels in the response variable. The primary drawback of MLM models is that they require larger sample sizes within each cluster, so may not work well when clusters are small.<\/p>\n\n\n

      \t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
      \n\t\t\t\t
      \n\t\t\t\t\t\t\t\t\t\n

      Both GEE and MLM are fairly easy to use in R. Below, I will walk through examples with the two most common kinds of correlated data: data with repeated measures from individuals and data collected from individuals with an important grouping variable (in this case, country). I will fit simple regression, GEE, and MLM models with each dataset, and will discuss which modeling technique is best for these different data types.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

      \n\t\t\t\t
      \n\t\t\t\t\t

      Example 1: Data from the Fragile Families & Child Wellbeing Study<\/h2>\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
      \n\t\t\t\t
      \n\t\t\t\t\t\t\t\t\t\n

      The data that I will be working with first comes from Years 9 and 15 of the Princeton University\u00a0Fragile Families & Child Wellbeing Study<\/a>, which follows the families of selected children born between 1998 and 2000 in major US cities. Data are publicly available, and can be accessed by submitting a brief request on the Fragile Families\u00a0Data and Documentation<\/a>\u00a0page. Since this study follows up with the same families year after year, data points from the same family units at different time points are not independent.<\/p>\n\n\n\n

      This dataset contains dozens of variables representing the health of wellbeing of participating children and their parents. Being in psychiatric epidemiology, I am primarily interested in examining the children\u2019s mental-wellbeing. Participating children are asked if they frequently feel sad, and I will be using answers to this \u201coften feeling sad\u201d question as my outcome. Since substance use is tied to poorer mental wellbeing among adolescents, I will be using variables representing alcohol and tobacco use as predictors*.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

      \n\t\t\t\t
      \n\t\t\t\t\t\t\t\t\t\n

      *Note: Models created in this article are for demonstration purposes only and should not be considered to be meaningful. I have not considered confounding, mediation, other model assumptions, or other possible data issues in the construction of these models.<\/p>\n\n\n\n

      First, let\u2019s load the packages that we\u2019ll be using. I\u2019ve loaded \u201ctidyverse\u201d to clean our data, \u201chaven\u201d because the data we\u2019ll be reading in comes in SAS format, \u201cgeepack\u201d to run our GEE model, and \u201clme4\u201d to run our multilevel model:<\/p>\n\n\n\n

      library(tidyverse)
      library(haven)
      library(geepack)
      library(lme4)<\/pre>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
      \n\t\t\t\t
      \n\t\t\t\t\t\t\t\t\t\n
      Now let\u2019s do some data cleaning to get these data ready for modeling!<\/p>\n\n\n\n
      Data from Years 9 and 15 are housed in separate SAS files (identifiable by the .sas7bdat extension), so we have one code chunk to read in and clean each file. This cleaning has to be done separately because variable names and coding differ slightly between study years (see the\u00a0Data and Documentation<\/a>\u00a0page for codebooks).<\/p>\n\n\n\n
      There are hundreds of variables included in the datasets, so we first select those that will be used in our model and assign meaningful variable names that are consistent across data frames. Next, we filter the data to only include individuals with complete data for our variables of interest (the code below excludes individuals with missing data for these variables as well as those who refused to answer).<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
      \n\t\t\t\t
      \n\t\t\t\t\t\t\t\t\t\n
      We then recode our variables in the standard 1 = \u201cyes\u201d, 0 = \u201cno\u201d format. For the \u201cfeel_sad\u201d variable, this also means dichotomizing a variable with 4 levels which represent varying degrees of sadness. We end up with a binary variable where 1 = \u201csad\u201d and 0 = \u201cnot sad.\u201d Some regression techniques can handle multiple levels in your response variable (MLM included), but I have binarized it here for simplicity. Finally, we create a \u201ctime_order\u201d variable indicating if the observation comes from the first or second round of the study.<\/p>\n\n\n\n
      year_9 = read_sas(\".\/data\/FF_wave5_2020v2_SAS.sas7bdat\") %>% 
        select(idnum, k5g2g, k5f1l, k5f1j) %>% 
        rename(\"feel_sad\" = \"k5g2g\",
               \"tobacco\" = \"k5f1l\",
               \"alcohol\" = \"k5f1j\") %>% 
        filter(
          tobacco == 1 | tobacco == 2,
          alcohol == 1 | alcohol == 2,
          feel_sad == 0 | feel_sad == 1 | feel_sad == 2 | feel_sad == 3
        ) %>% 
        mutate(
          tobacco = ifelse(tobacco == 1, 1, 0),
          alcohol = ifelse(alcohol == 1, 1, 0),
          feel_sad = ifelse(feel_sad == 0, 0, 1),
          time_order = 1
        )year_15 = read_sas(\".\/data\/FF_wave6_2020v2_SAS.sas7bdat\") %>% 
        select(idnum, k6d2n, k6d40, k6d48) %>% 
        rename(\"feel_sad\" = \"k6d2n\",
               \"tobacco\" = \"k6d40\",
               \"alcohol\" = \"k6d48\") %>% 
        filter(
          tobacco == 1 | tobacco == 2,
          alcohol == 1 | alcohol == 2,
          feel_sad == 1 | feel_sad == 2 | feel_sad == 3 | feel_sad == 4
        ) %>% 
        mutate(
          tobacco = ifelse(tobacco == 1, 1, 0),
          alcohol = ifelse(alcohol == 1, 1, 0),
          feel_sad = ifelse(feel_sad == 4, 0, 1),
          time_order = 2
        )<\/pre>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
      \n\t\t\t\t
      \n\t\t\t\t\t\t\t\t\t\n
      We then combine data from Years 9 and 15 by stacking our two cleaned data frames using rbind(). The rbind() function works well here because both data frames now share all variable names. We next transform the \u201cidnum\u201d variable (which identifies unique family units) into a numeric variable so that it can be properly used to sort the data in the final code chunk. This step is necessary because the geeglm() function that we will be using to run the GEE model assumes that the data frame is sorted first by a unique identifier (in this case, \u201cidnum\u201d), and next by the order of observations (indicated here by the new \u201ctime_order\u201d variable).<\/p>\n\n\n\n
      fragile_families = rbind(year_9, year_15) %>% 
        mutate(
          idnum = as.numeric(idnum)
        )fragile_families = 
        fragile_families[
        with(fragile_families, order(idnum)),
      ]<\/pre>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
      \n\t\t\t\t
      \n\t\t\t\t\t\t\t\t\t\n
      The above code produces the following cleaned data frame, which is now ready to be used for regression modeling:<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t\t<\/div>\n\t\t<\/div>\n\t\t\t\t\t<\/div>\n\t\t<\/section>\n\t\t\t\t
      \n\t\t\t\t\t\t
      \n\t\t\t\t\t
      \n\t\t\t
      \n\t\t\t\t\t\t
      \n\t\t\t\t
      \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"\"\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
      \n\t\t\t\t
      \n\t\t\t\t\t\t\t\t\t\n
      Let\u2019s fit our models:<\/p>\n\n\n\n
        \n
      1. Simple Logistic Regression<\/strong><\/li>\n<\/ol>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\n\n
        First, we use the glm() function to fit a simple logistic regression model using the \u201cfragile_families\u201d data. Since we have a binary outcome variable, \u201cfamily = binomial\u201d is used to specify that logistic regression should be used. We also use tidy() from the \u201cbroom\u201d package to clean up the model output. We are creating this model for comparison purposes only \u2014 as indicated before, the independence assumption has been violated and the standard errors associated with this model will not be valid!<\/p>\n\n\n\n
        glm(formula = feel_sad ~ tobacco + alcohol, 
            family = binomial, data = fragile_families) %>% 
          broom::tidy()<\/pre>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\n
        The above code produces the following output, which the subsequent modeling approaches will be compared to. Tobacco and alcohol use both appear to be significant predictors of sadness in participating children.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"\"\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\n
        2. Generalized Estimating Equations<\/strong><\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\n
        The syntax used to specify a GEE model using the geeglm() function from the \u201cgeepack\u201d package is fairly similar to that used with the standard glm() function. The \u201cformula\u201d, \u201cfamily\u201d, and \u201cdata\u201d are arguments are exactly the same for both functions. What\u2019s new are the \u201cid,\u201d \u201cwaves,\u201d and \u201ccorstr\u201d arguments (see\u00a0package documentation<\/a>\u00a0for all available arguments). The unique identifier that links observations from the same subject is specified in the \u201cid\u201d argument. In this case the ID is \u201cidnum,\u201d the unique identifier assigned to each family participating in the study. The \u201ctime_order\u201d variable created during data cleaning comes into play in the \u201cwaves\u201d argument, where it indicates the order in which observations were made. Finally, \u201ccorstr\u201d can be used to specify the within-subject correlation structure. \u201cIndependence\u201d is actually the default input for this argument, and it makes sense in this context because it is useful when clusters are small. However, \u201cexchangeable\u201d can be specified when all observations within a subject can be considered to be equally correlated, and \u201car1\u201d is best when the internal correlations change over time. Information on choosing the right correlation structure can be found\u00a0here\u00a0<\/a>and\u00a0here<\/a>.<\/p>\n\n\n\n
        geeglm(formula = feel_sad ~ tobacco + alcohol, 
               family = binomial, id = idnum, data = fragile_families, 
               waves = time_order, corstr = \"independence\") %>% 
          broom::tidy()<\/pre>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\n
        Our GEE model gives us the following output:<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"\"\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\n
        As you can see, our beta estimates are exactly the same as those produced using glm(), but standard error differs slightly now that the correlations in the data have been accounted for. While tobacco and alcohol are still significant predictors of sadness, the p-values are somewhat different**. If these p-values were closer to 0.05, having accurate standard error measurements could easily push a p-value over or under the level of significance.<\/p>\n\n\n\n
        **Note: The test statistics for GEE and logistic regression look drastically different, but this is only because the test statistic provided in the logistic regression output is a Z-statistic and the test statistic provided in the GEE output is a Wald statistic. The Z-statistic is calculated by dividing the estimate by the standard error, while the Wald statistic is calculated by squaring the result of dividing the estimate by the standard error. The two values are therefore mathematically related, and by taking the square root of the values in the GEE \u201cstatistic\u201d column you will see a much more moderate change from the initial Z-statistics.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\n
        With the geeglm() function, it is also important to verify that your clusters have been properly recognized. You can do this by running the above code without the broom::tidy() step, so:<\/p>\n\n\n\n
        geeglm(formula = feel_sad ~ tobacco + alcohol, 
               family = binomial, id = idnum, data = fragile_families, 
               waves = time_order, corstr = \"independence\")<\/pre>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\n
        This code produces the output shown below. You want to look to the last line of the output, where \u201cNumber of clusters\u201d and \u201cMaximum cluster size\u201d are described. We had 2 observations for several thousand individuals, so these values make sense in the context of our data and indicate that clusters were registered correctly by the function. If, however, the number of clusters is equal to the number of rows in your dataset, something is not working properly (most likely the sorting of your data is off).<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"\"\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t\t<\/div>\n\t\t<\/div>\n\t\t\t\t\t<\/div>\n\t\t<\/section>\n\t\t\t\t
        \n\t\t\t\t\t\t
        \n\t\t\t\t\t
        \n\t\t\t
        \n\t\t\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\n
        3. Multilevel Modeling<\/strong><\/p>\n\n\n\n
        Next, let\u2019s fit a multilevel model using glmer() from the\u00a0lme4 package<\/a>. Again, the required code is almost identical to that used for logistic regression. The only required change is specifying random slopes and intercepts in the formula argument. This is done with the \u201c(1 | idnum)\u201d bit of code, which follows the following structure: (random slopes | random intercepts). The grouping variable, in this case \u201cidnum,\u201d is specified to the right of the | as \u201crandom intercepts,\u201d and the \u201c1\u201d indicates that we don\u2019t want the predictors\u2019 effects to vary across groups. A\u00a0useful blog post<\/a>\u00a0by Rense Nieuwenhuis provides various examples of this glmer() syntax.<\/p>\n\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t
        The lme4 package is not compatible with the broom package, so instead we pull the model\u2019s coefficients after creating a list with a summary of the model\u2019s output.<\/p>\n\n\n\n
        mlm = summary(glmer(formula = 
                            feel_sad ~ tobacco + alcohol + (1 | idnum), 
                            data = fragile_families, family = binomial))mlm$coefficients<\/pre>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\n
        Again, the output is similar to that of the simple logistic regression model, and both tobacco and alcohol use are still significant predictors of sadness. Estimates vary slightly from those produced using the glm() and geeglm() functions because groupings in the data are no longer ignored or treated as an annoyance to be addressed by correcting standard error; instead, they are now incorporated as an important part of the model. Standard error estimates are higher for all estimates in comparison to those produced through logistic regression, and Z- and p-values remain similar but reflect these important changes in the estimate and standard error values.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"\"\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t
        Example 2: Data from the Global School-Based Student Health Survey (GSHS)<\/h2>\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\n
        The second dataset that we will walk through comes from the WHO\u2019s\u00a0Global School-Based Student Health Survey<\/a>\u00a0(GSHS). This survey is conducted among schoolchildren aged 13\u201317 with the goals of helping countries to determine health priorities, establishing the prevalences of health-related behaviors, and facilitating direct comparison of these prevalences across nations. We will be using data from two countries,\u00a0Indonesia<\/a>\u00a0and\u00a0Bangladesh<\/a>, which can be downloaded directly from these countries\u2019 respective descriptive pages.<\/p>\n\n\n\n
        The data are cross-sectional: an identical survey was conducted one time among schoolchildren in both nations. I am interested in using variables from this dataset to describe the relationship between whether or not a child has friends, whether or not the child is bullied (my predictors) and whether or not the child has seriously contemplated suicide (my outcome). It is likely that these relationships differ between the two countries and that children are more similar to other children from the same country. Therefore, knowing whether a child is from Indonesia or Bangladesh provides important information about that child\u2019s responses and the assumption of independent observations is violated.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\n
        Let\u2019s load packages again:<\/p>\n\n\n\n
        library(tidyverse)
        library(haven)
        library(lme4)
        library(gee)<\/pre>\n\n\n\n
        Note that the \u201cgeepack\u201d package has been replaced with the \u201cgee\u201d package. The \u201cgee\u201d package is easier to use (in my opinion) with data that is clustered by a grouping variable such as country rather than within an individual who has multiple observations.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\n
        Next, let\u2019s load in the data (which is also in SAS format, so we use the \u201chaven\u201d package again) and conduct some basic cleaning. Data cleaning here follows a similar structure to the procedure used with the Fragile Families & Child Wellbeing Study data: important variables are selected and assigned meaningful, consistent names, and a new variable is created to indicate which cluster an observation belongs to (in this case the new \u201ccountry\u201d variable).<\/p>\n\n\n\n
        indonesia = read_sas(\".\/data\/IOH2007_public_use.sas7bdat\") %>% 
          select(q21, q25, q27) %>% 
          rename(
            \"bullied\" = \"q21\",
            \"suicidal_thoughts\" = \"q25\",
            \"friends\" = \"q27\"
          ) %>% 
          mutate(
            country = 1,
          )bangladesh = read_sas(\".\/data\/bdh2014_public_use.sas7bdat\") %>% 
          select(q20, q24, q27) %>% 
          rename(
            \"bullied\" = \"q20\",
            \"suicidal_thoughts\" = \"q24\",
            \"friends\" = \"q27\"
          ) %>% 
          mutate(
            country = 2
          )<\/pre>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t
        Again, the two data frames are stacked together. Since variables were coded consistently during collection in both countries, some cleaning can be conducted only once using this combined dataset. Missing data is eliminated, and all variables are converted from string format to numeric. Finally, variables are mutated into a consistent, binarized format.<\/p>\n\n\n\n
        survey = rbind(indonesia, bangladesh) %>% 
          mutate(
            suicidal_thoughts = as.numeric(suicidal_thoughts),
            friends = as.numeric(friends),
            bullied = as.numeric(bullied),
            suicidal_thoughts = ifelse(suicidal_thoughts == 1, 1, 0),
            friends = ifelse(friends == 1, 0, 1),
            bullied = ifelse(bullied == 1, 0, 1)
          ) %>% 
          drop_na()<\/pre>\n\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\n
        Our cleaned data frame now looks like this:<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"\"\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
        \n\t\t\t\t
        \n\t\t\t\t\t\t\t\t\t
        Let\u2019s fit our models:<\/p>\n\n\n\n
          \n
        1. Simple Logistic Regression<\/strong><\/li>\n<\/ol>\n\n\n\n
          With the exception of variable names and the data specified, the glm() code remains identical to that used with the Fragile Families study data.<\/p>\n\n\n\n
          glm(formula = suicidal_thoughts ~ bullied + friends, 
              family = binomial, data = survey) %>% 
            broom::tidy()<\/pre>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
          \n\t\t\t\t
          \n\t\t\t\t\t\t\t\t\t\n
          Unsurprisingly, whether or not a child has friends and whether or not a child is bullied are both significant predictors of the presence of suicidal thoughts in this sample.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
          \n\t\t\t\t
          \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"\"\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
          \n\t\t\t\t
          \n\t\t\t\t\t
          2. Generalized Estimating Equations<\/h2>\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
          \n\t\t\t\t
          \n\t\t\t\t\t\t\t\t\tThe gee() function in the gee package allows us to easily use GEE with our survey data. This function is a better fit than the previously used geeglm() function as data are not correlated over time, but rather by a separate variable that can be indicated with the \u201cid\u201d argument (in this case, \u201ccountry\u201d). The formula and family arguments remain identical to those used with the glm() function, and the \u201ccorstr\u201d argument used with the geeglm() function is the same here as well. However, unlike the geepack package, the gee package is not compatible with the broom::tidy() function so output is viewed using the summary() function instead.\n\n\n
          gee = gee(suicidal_thoughts ~ bullied + friends, data = survey, 
                    id = country, family = binomial, 
                    corstr = \"exchangeable\")summary(gee)<\/pre>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
          \n\t\t\t\t
          \n\t\t\t\t\t\t\t\t\t\n
          One of the reasons that I particularly like the gee() function is that the naive standard error and Z-test statistics are actually included in the output (naive meaning that these values are produced by regression where clustering is not accounted for \u2014 you\u2019ll see that these are exactly the same as those produced by the glm() function above). You\u2019ll notice drastic changes in the standard errors and Z-test statistics produced using GEE (\u201cRobust\u201d), although both of our predictors remain significant. It appears that accounting for within-country correlation has allowed for much lower standard errors to be used.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
          \n\t\t\t\t
          \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"\"\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
          \n\t\t\t\t
          \n\t\t\t\t\t\t\t\t\t\n
          3. Multilevel Modeling***<\/strong><\/p>\n\n\n\n
          ***Note: As noted above, models are for demonstration purposes only and are not necessarily valid. In this case, we would want more groups than two for our MLM model (meaning data from additional countries). If you are really only using two groups with MLM models, you should consider a\u00a0small sample size correction<\/a>.<\/p>\n\n\n\n
          Finally, we try MLM with the survey dataset. The code is exactly the same as that used with the Fragile Families study data, but with the new formula, grouping variable, and dataset specified.<\/p>\n\n\n\n
          mlm = summary(glmer(formula = 
                              suicidal_thoughts ~ bullied + friends + 
                              (1 | country), 
                              data = survey, family = binomial))mlm$coefficients<\/pre>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
          \n\t\t\t\t
          \n\t\t\t\t\t\t\t\t\t\n
          Again, beta estimates and standard error estimates are now adjusted slightly from those produced using glm(). Z- and p-values associated with the \u201cbullied\u201d and \u201cfriends\u201d variables are slightly smaller, although bullying and having friends remain significant predictors of suicidal thoughts.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
          \n\t\t\t\t
          \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"\"\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
          \n\t\t\t\t
          \n\t\t\t\t\t
          Which model is best for these examples?<\/h1>\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t\t<\/div>\n\t\t<\/div>\n\t\t\t\t\t<\/div>\n\t\t<\/section>\n\t\t\t\t
          \n\t\t\t\t\t\t
          \n\t\t\t\t\t
          \n\t\t\t
          \n\t\t\t\t\t\t
          \n\t\t\t\t
          \n\t\t\t\t\t\t\t\t\t\n
          Data from Princeton University\u2019s Fragile Families & Child Wellbeing Study would be best represented using GEE. This is due to the maximum cluster size of 2 observations, the fact that individual families have multiple data points over time, and the fact that we were more interested in accounting for grouping in the standard error estimates than actually assessing differences between families.<\/p>\n\n\n\n
          Multilevel modeling is most appropriate for data from the Global School-Based Student Health Survey (GSHS) because the data were collected cross-sectionally and can be divided into two large clusters. Additionally, the output could be further explored to determine both within- and between-group variances, and we might be interested in relationships both within and across countries.<\/p>\n\n\n\n
          How you account for violations of the independent observations assumption will depend on the structure of your data and your general knowledge of the data collection process, as well as whether or not you consider the correlation to be an annoyance to adjust for or something meaningful to explore.<\/p>\n\n\n\n
          In conclusion, regression is flexible and certain regression models can handle correlated data.\u00a0<\/strong>However, it is always important to check the assumptions of a given technique and to make sure that your analytic strategy is appropriate for your data.<\/p>\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t\t<\/div>\n\t\t<\/div>\n\t\t\t\t\t<\/div>\n\t\t<\/section>\n\t\t\t\t<\/div>\n\t\t","protected":false},"excerpt":{"rendered":"
          Regression is flexible and certain regression models can handle correlated data. However, it is always important to check the assumptions of a given technique and to make sure that your analytic strategy is appropriate for your data. <\/p>\n","protected":false},"author":924,"featured_media":9876,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"content-type":"","footnotes":""},"categories":[187],"tags":[671,670,669],"ppma_author":[3685],"class_list":["post-9875","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-bigdata-cloud","tag-correlated-data","tag-models","tag-regression"],"authors":[{"term_id":3685,"user_id":924,"is_guest":0,"slug":"emily-halford","display_name":"Emily Halford","avatar_url":"https:\/\/www.experfy.com\/blog\/wp-content\/uploads\/2020\/09\/Emily-Halford-150x150.jpg","user_url":"http:\/\/nysi.org","last_name":"Halford","first_name":"Emily","job_title":"","description":"Emily Halford is Data Analyst at New York State Psychiatric Institute."}],"_links":{"self":[{"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/posts\/9875","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/users\/924"}],"replies":[{"embeddable":true,"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/comments?post=9875"}],"version-history":[{"count":9,"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/posts\/9875\/revisions"}],"predecessor-version":[{"id":33808,"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/posts\/9875\/revisions\/33808"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/media\/9876"}],"wp:attachment":[{"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/media?parent=9875"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/categories?post=9875"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/tags?post=9875"},{"taxonomy":"author","embeddable":true,"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/ppma_author?post=9875"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}