How can I replace an ANOVA with dependent variables in R?

Question

I am conducting a study on wildfires, and so I have made some measures on burned trees. I have 4 variables (level of damage at the base, level of damage on the bark, species and diameter of the tree). I have converted them into numerical values (categories from 1 to 3 depending on the severity of the burn, categories from 1 to 3 depending on te diameter and categories from 1 to 5 depending on the species). My data looks like that, with each line corresponding to a tree measured.

 ALIVE SPECIES DIAM BASE BARK
1     1   1      2    3    1
2     1   3      2    2    1
3     1   1      2    1    1
4     1   2      1    2    1
5     1   2      1    2    1

I have tried to run a 4-way ANOVA in R to test which variable has the most impact on the survival of the trees, but I then realised that the variables needed to be independant, which they are not.

I then tried to run 4 one-way anova test, one for each variable and got these results for the base. I followed this method : https://www.scribbr.com/statistics/anova-in-r/

 Df Sum Sq Mean Sq F value   Pr(>F)    
BASE          2  2.154  1.0770   9.721 8.31e-05 ***
Residuals   277 30.689  0.1108    

Then, like the method said, i checked for homoscedasticity, but it does not seem like the results are as they should be (horizontal lines) :

But it seems strange to me to compare the variance of only 1 variable, which is what an homoscedasticity test does?

Do I need to use a different tool? Is there a method which could tell me what combination of variables has the most impact on the survival rate?

Thanks in advance! (sorry if unclear, english is not my first langage, and i am new to R)

Answer 1

Your English is perfectly understandable.

A couple points: First, the predictor variables in an ANOVA do not need to be completely independent. They can't be exactly colinear, or the model will "blow up". How much colinearity is acceptable is debated, but there are guidelines related to either VIF or condition indexes. I prefer the latter, and a condition index over 10 is considered a bit problematic and over 30 quite problematic.

Second, converting diameter to a factor is probably a mistake. I'm also interested in how and why you converted the other variables (how is damage measured?) And, when you run ANOVA, it will assume that the levels of the variables are nominal, but they are ordinal. Ordinal IVs are tricky. (Species is nominal, of course).

Third, some of the odd pattern in your plots is due to the predictor variable (here base) taking on only a few values.

Finally, if, after working on the above, you still have heteroscedasticity and other issues with the residuals, you can consider quantile regression or robust regression.