I have a big dataset with different kinds of human interventions in natural freshwater environments (categorical variable, one column) and the responses as control/ with intervention. It is a reunion of studies published data.

This dataset also has other columns with relevant characterization such as kind of environment (lake, river...), other associated environmental conditions, and year of measurement, for example (categorical variables) and I would also like to explore them. My question is: Is there a way to compare a control/ with intervention set of data for the interference of more than one categorical variable? The same "pairs" of control/ with intervention response variables are subjected to more than one categorical variable. Can I analyze all of it at once?

Dataset example:

structure(list(ID = c("Bayram_&_Kenanoglu_2016", "Bayram_&_Kenanoglu_2016", 
"Bayram_&_Kenanoglu_2016", "Zheng_2021", "Zheng_2021", "Chen_2018", 
"Chen_2018", "Baborowsky_et_al_2004", "Rosado_Berrios_Bouldin_2016", 
"Rosado_Berrios_Bouldin_2016", "Rosado_Berrios_Bouldin_2016", 
"Rosado_Berrios_Bouldin_2016", "Rosado_Berrios_Bouldin_2016", 
"Rosado_Berrios_Bouldin_2016", "Hamers_et_al_2015", "Hamers_et_al_2015", 
"Li_et_al_2021", "Li_et_al_2021", "Cao_et_al_2017", "Cao_et_al_2017", 
"Cao_et_al_2017", "Cao_et_al_2017", "Cao_et_al_2017", "Trentman_et_al_ 2021", 
"Hasenmueller_et_al_2017", "Hasenmueller_et_al_2017", "Rodrigues_et_al_2017", 
"Dewey_et_al_2020", "Dewey_et_al_2020", "Dewey_et_al_2020", "Shrestha_et_al_2017", 
"Shrestha_et_al_2017", "Bärlocher_et_al_2010", "Bärlocher_et_al_2010", 
"Bärlocher_et_al_2010", "Bärlocher_et_al_2010", "Bärlocher_et_al_2010", 
"Bärlocher_et_al_2010", "Vergilio_et_al_2021", "Vergilio_et_al_2021", 
"Pizarro_et_al_2014", "Pizarro_et_al_2014", "Pizarro_et_al_2014", 
"Pizarro_et_al_2014", "Pizarro_et_al_2014", "Pizarro_et_al_2014", 
"Fan_2011", "Chang_et_al_2017", "Chang_et_al_2017"), Year = c(2016L, 
2016L, 2016L, 2021L, 2021L, 2018L, 2018L, 2004L, 2016L, 2016L, 
2016L, 2016L, 2016L, 2016L, 2015L, 2015L, 2021L, 2021L, 2017L, 
2017L, 2017L, 2017L, 2017L, 2021L, 2017L, 2017L, 2017L, 2020L, 
2020L, 2020L, 2017L, 2017L, 2010L, 2010L, 2010L, 2010L, 2010L, 
2010L, 2021L, 2021L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 
2011L, 2017L, 2017L), Lat = c("41,36530848", "41,36530848", "41°19´18 .20''", 
"23°0'0'' N", "23°0'0'' N", "24,49893013", "24,49893013", "52,16436256", 
"34,93392781", "34,93392781", "34,93392781", "34,93392781", "34,93392781", 
"34,93392781", "50,78239805", "51,45862307", "37,79047935", "37,79047935", 
"33,27060865", "33,27060865", "33,27060865", "33,27060865", "33,27060865", 
"41,21300785", "38,58109023", "38,62050268", "-22,51962467", 
"33°30' 49.00'' ", "33°30' 49.00'' ", "33°30' 49.00'' ", "33,4784627", 
"33,4784627", "09.18.790", "09.18.790", "09.18.790", "09.16.728 ", 
"09.16.728 ", "09.16.728 ", "-20,23527715", "-20,23527715", "-35,08998953", 
"-35,08998953", "-35,08998953", "-35,08998953", "-35,08998953", 
"-35,08998953", "24,90596076", "24°22'20.8'N", "24°22'20.8'N"
), Long = c("41,67889289", "41,67889289", "41°21'02.40''", "90°30' 00'' E", 
"90°30' 00'' E", "117,7916689", "117,7916689", "11,68275488", 
"-91,34544184", "-91,34544184", "-91,34544184", "-91,34544184", 
"-91,34544184", "-91,34544184", "5,698366157", "6,723573688", 
"119,2299644", "119,2299644", "118,6929425", "118,6929425", "118,6929425", 
"118,6929425", "118,6929425", "-86,04450487", "-90,31187387", 
"-90,34045645", "-48,53586955", "89°25'40.00''", "89°25'40.00''", 
"89°25'40.00''", "-91,1356069", "-91,1356069", " 079.37.772", 
" 079.37.772", " 079.37.772", "079.41.667", "079.41.667", "079.41.667", 
"-43,41951321", "-43,41951321", "-71,68723385", "-71,68723385", 
"-71,68723385", "-71,68723385", "-71,68723385", "-71,68723385", 
"121,5532293", "121°18'40.0''E", "121°18'40.0''E"), Country = c("", 
"", "", "", "", "", "", "", "U.S.A.", "U.S.A.", "U.S.A.", "U.S.A.", 
"U.S.A.", "U.S.A.", "Belgium", "Germany", "China", "China", "China", 
"China", "China", "China", "China", "U.S.A.", "U.S.A.", "U.S.A.", 
"Brazil", "U.S.A.", "U.S.A.", "U.S.A.", "U.S.A.", "U.S.A.", "Panama", 
"Panama", "Panama", "Panama", "Panama", "Panama", "Brazil", "Brazil", 
"Chile", "Chile", "Chile", "Chile", "Chile", "Chile", "Taiwan", 
"Taiwan", "Taiwan"), year_intervention = c(2011L, 2011L, 2011L, 
2011L, 2012L, 2014L, 2014L, 2002L, 2012L, 2013L, 2014L, 2012L, 
2013L, 2014L, 1993L, 1993L, 2019L, 2019L, 2012L, 2012L, 2012L, 
2012L, 2012L, 2014L, 2008L, 2008L, 2014L, 2011L, 2011L, 2011L, 
2013L, 2013L, 2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 2015L, 
2015L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2005L, 2011L, 
2011L), place = c("Borçka Dam Reservoir", "Borçka Dam Reservoir", 
"Borçka Dam Reservoir", "Padma River", "Padma River", "Jiulong River", 
"Jiulong River", "Elbe River", "Maloy Bayou Cache River", "Maloy Bayou Cache River", 
"Maloy Bayou Cache River", "Maloy Bayou Cache River", "Maloy Bayou Cache River", 
"Maloy Bayou Cache River", "River Meuse", "River Rhine", "Yellow River", 
"Yellow River", "Hongze Lake", "Hongze Lake", "Hongze Lake", 
"Hongze Lake", "Hongze Lake", "Shatto ditch", "River des Perez", 
"Black Creek", "Barra Bonita", "Upper Gulf Coastal Plain", "Upper Gulf Coastal Plain", 
"Upper Gulf Coastal Plain", "River Mississipi", "River Mississipi", 
"Panama Canal", "Panama Canal", "Panama Canal", "Panama Canal", 
"Panama Canal", "Panama Canal", "Doce River Basin", "Doce River Basin", 
"South Chile river catchments", "South Chile river catchments", 
"South Chile river catchments", "South Chile river catchments", 
"South Chile river catchments", "South Chile river catchments", 
"Feitsui Reservoir", "Chichiawan Stream", "Chichiawan Stream"
), nature_pollution = c("inorganic", "inorganic", "inorganic", 
"inorganic", "inorganic", "inorganic", "inorganic", "inorganic", 
"inorganic_organic", "inorganic_organic", "inorganic_organic", 
"inorganic_organic", "inorganic_organic", "inorganic_organic", 
"inorganic_organic", "inorganic_organic", "inorganic_organic", 
"inorganic_organic", "inorganic", "inorganic", "inorganic", "inorganic", 
"inorganic", "inorganic_organic", "inorganic_organic", "inorganic_organic", 
"inorganic_organic", "inorganic_organic", "inorganic_organic", 
"inorganic_organic", "inorganic_organic", "inorganic_organic", 
"inorganic_organic", "inorganic_organic", "inorganic_organic", 
"inorganic_organic", "inorganic_organic", "inorganic_organic", 
"inorganic", "inorganic", "inorganic_organic", "inorganic_organic", 
"inorganic_organic", "inorganic_organic", "inorganic_organic", 
"inorganic_organic", "inorganic", "inorganic_organic", "inorganic_organic"
), kinf_of_environment = c("river", "river", "river", "river", 
"river", "estuary", "estuary", "river", "river", "river", "river", 
"river", "river", "river", "river", "river", "estuary", "estuary", 
"lake", "lake", "lake", "lake", "lake", "ditch", "river", "river", 
"dam", "creek", "creek", "creek", "river", "river", "stream", 
"stream", "stream", "stream", "stream", "stream", "river", "river", 
"river", "river", "river", "river", "river", "river", "reservoir", 
"stream", "stream"), trophic_state = c(NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA), kind_of_measurem = c("dif_places", 
"dif_places", "dif_places", "repeat_time", "repeat_time", "repeat_time", 
"repeat_time", "repeat_time", "dif_places", "dif_places", "dif_places", 
"dif_places", "dif_places", "dif_places", "repeat_time", "repeat_time", 
"repeat_time", "repeat_time", "dif_places", "dif_places", "dif_places", 
"dif_places", "dif_places", "repeat_time", "dif_places", "dif_places", 
"dif_places", "dif_places", "dif_places", "dif_places", "dif_places", 
"dif_places", "dif_places", "dif_places", "dif_places", "dif_places", 
"dif_places", "dif_places", "repeat_time", "repeat_time", "dif_places", 
"dif_places", "dif_places", "dif_places", "dif_places", "dif_places", 
"repeat_time", "repeat_time", "repeat_time"), cause_variation = c("dam", 
"dam", "mining", "climate_change", "climate_change", "typhoon", 
"storm", "flood", "erosion", "erosion", "erosion", "erosion", 
"erosion", "erosion", "flood", "flood", "artificial_regulation_scheme", 
"artificial_regulation_scheme", "dredging", "dredging", "dredging", 
"dredging", "dredging", "erosion_reversion_cover", "urbanization", 
"urbanization", "pollution", "deforestation", "deforestation", 
"deforestation", "agriculture", "agriculture", "agriculture", 
"agriculture", "agriculture", "urbanization", "urbanization", 
"urbanization", "mining_dam_break", "mining_dam_break", "agriculture", 
"agriculture", "pasture", "agriculture", "agriculture", "pasture", 
"typhoon", "dam_removal", "dam_removal"), particle = c("suspended_solids", 
"suspended_solids", "tailing", "suspended_solids", "suspended_solids", 
"suspended_solids", "suspended_solids", "suspended_solids", "suspended_solids", 
"suspended_solids", "suspended_solids", "suspended_solids", "suspended_solids", 
"suspended_solids", "suspended_solids", "suspended_solids", "suspended_solids", 
"suspended_solids", "suspended_solids", "suspended_solids", "suspended_solids", 
"suspended_solids", "suspended_solids", "suspended_solids", "suspended_solids", 
"suspended_solids", "suspended_solids", "suspended_solids", "suspended_solids", 
"suspended_solids", "suspended_solids", "suspended_solids", "suspended_solids", 
"suspended_solids", "suspended_solids", "suspended_solids", "suspended_solids", 
"suspended_solids", "tailing", "tailing", "suspended_solids", 
"suspended_solids", "suspended_solids", "suspended_solids", "suspended_solids", 
"suspended_solids", "suspended_solids", "suspended_solids", "suspended_solids"
), mean_contr = c(27, 28.2, 28.2, 95.8, 95.8, 53.5, 38.33, 13.86, 
59.23, 77.98, 98.01, 34.97, 35.11, 40.19, 6.75, 30.83, 63.14, 
62.69, 9.44, 9.44, 9.44, 9.44, 9.44, 97.83, 10, 10, 0.81, 342.6, 
342.6, 342.6, 156, 75, 2.3, 9.68, 70, 2.3, 9.68, 70, 56, 20, 
7.2, 7.2, 7.2, 16.8, 16.8, 16.8, 1.5, 0.39, 0.58), mean_intervention = c(46, 
40.9, 16991, 156, 169.5, 291.71, 198.54, 42.44, 57.59, 104.67, 
68.29, 30.49, 42.79, 29.78, 132, 123.25, 876.42, 198.97, 36.7, 
60.5, 55.61, 44, 34.08, 16.6, 14, 656, 21.91, 1713, 54.1, 154.7, 
303, 103, 5.94, 11.86, 86, 6.83, 15.68, 130, 868, 2956, 19.8, 
12, 8.7, 19.8, 12, 8.7, 9.2, 0.63, 0.75)), class = "data.frame", row.names = c(NA, 

I thought of something like an NMDS plot with polygons separating the groups of categorical variables, but since my response variables are control/ with intervention it seems to me that the points are necessarily going to form two groups and the polygons would be biased. I really have no clue on how to achieve this. Any suggestions?

  • $\begingroup$ This is hard to answer well without more detail about your data. But it sounds like a multiple regression would be possible, and this simpler and easier to understand & evaluate than multivariate analysis. $\endgroup$
    – mkt
    Jul 30, 2022 at 6:05
  • $\begingroup$ I am adding a link to the dataset to help comprehend what I need. $\endgroup$ Jul 30, 2022 at 13:05
  • $\begingroup$ Instead of a link, please dput() your data, or enough of it so the question can be understood: stats.meta.stackexchange.com/a/5931/121522 $\endgroup$
    – mkt
    Jul 30, 2022 at 14:04
  • $\begingroup$ I just added the dput() output. Hope it is correct now. $\endgroup$ Jul 31, 2022 at 11:07
  • $\begingroup$ Thanks. Is this your entire dataset, or a subset? $\endgroup$
    – mkt
    Aug 1, 2022 at 12:39

1 Answer 1


A regression (a.k.a. an ANOVA in this case) would be a simple and appropriate way to analyse this data. You would need to reorganise this data frame a little to analyse it appropriately, by making the measured outcome a single column (not separate columns for control and treatment), and then a new column that indicates whether the outcome corresponds to control or treatment.

After that, the simplest analysis you can run is just a regression with the measured outcome as response and treatment type (control/treatment) as the predictor. In R, it would be as simple as summary(lm(outcome_vale ~ treatment_type)). If you have much more data than you have presented here, you can try multiple regression with additional predictors (such as kind_of_environment) in the regression formula.

You will need to consider both domain expertise and the limitations of the data to choose which predictors to include in the model, whether to include interactions, and other complexities such as the functional form of the response if a predictor is continuous. The more data you have, the more power you have to detect weak or complex effects.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.