Anova for multiple subgroups

Question

To compare two models using anova in R, I can simply use anova(model1,model2) and conclude whether they are significantly different or not based on the p-value returned. However, currently I have data such that:

Cities     Localities     x1     x2     y1
City1           .          .      .      .
City1           .          .      .      .
City2           .          .      .      .
City2           .          .      .      .
City2           .          .      .      .
City3           .          .      .      .
City4           .          .      .      .
City4           .          .      .      .
City4           .          .      .      .

and so on. I want to show that there is variability between cities as groups without comparing the models for each subgroup pairwise. What can I do?

Will anova(models) where models=c(model1,model2,...) vector of fitted models corresponding to each city work? Or should I perform anova on the full dataset. In that case, how do I define the subdivisions?