How to use all observations per site in a model? I want to model plant traits as a function of environmental variables. For example, tree height as a function of fire frequency. I'm doing this to test the effects of fires on plant traits (and not to predict traits at different scenarios).
Between 15-20 individuals were measured per site, in 8 sites - so I have 15-20 values of height in each site, and one fire frequency value per site.
I started by doing Spearman correlations as a preliminary approach, using the mean height per site, but I would like to use an approach where I can use all the height values, so that all information is used.
I have seen this previous question where it is advised to try hierarchical models, using the response variable as a group-level predictor.
In this other question it is mentioned the use of environmental variables as fixed effects, and site as a random effect - I suppose here site is the grouping variable.
So are both ways correct? Or should I always introduce "site" in my model to group observations?
 A: As pointed out by Robert, Site is a grouping factor in your study.  To formulate the appropriate model, you will however need to determine whether you can treat Site as a fixed or random grouping factor in your modelling.
Site as a fixed grouping factor
If you were to repeat your study again, would you select the exact same 8 sites as before because these sites are the only ones you are interested in? If yes, you should treat Site as a fixed grouping factor.  That means that you could formulate your models as linear regression models using the lm() function of R:
# effect of fire_frequency on tree_height is assumed to be 
# the same across all 8 sites
m1 <- lm(tree_height ~ fire_frequency + Site, data = yourdata)


# effect of fire_frequency on tree_height is assumed to be 
# different across sites
m2 <- lm(tree_height ~ fire_frequency*Site, data = yourdata)

Site as a random grouping factor
If you were to repeat your study again, would you select the exact same 8 sites as before because these sites are the only ones you are interested in? If no, you should treat Site as a random grouping factor, since the 8 sites were selected to be representative of a larger set of sites you are truly interested in (ideally, they would have been selected at random from that larger set of sites). That means that you could formulate your model as a linear mixed effects regression model using the lmer() function of R:
library(lme4)

m <- lmer(tree_height ~ fire_frequency + (1|Site), data = yourdata)

A third possibility would be to use a GEE model - especially since your fire_frequency variable is a site-level predictor and GEE models can offer a more natural interpretation of its effect than mixed effects models.
A: You have trees ("individuals") clustered within sites, so you are correct that site is a grouping variable. So you have repeated measures within site, so random intercepts for site would be a good approach:
tree_height ~ fire_frequency + (1 | treeID)

As far as I can tell, the approach in both links are the same. The altenative would be to fit fixed effects for site.
