From the question, it's a little difficult for me to tell whether you're interested in comparing group means, or instead in estimate the effect of some other covariate for each group. I'll assume the latter.
The right choice is going to depend a bit on the features of your dataset, but my first thought would be to try estimating a model with an indicator for each group, interacted with the covariate of interest.
Here's an example using only two groups -- follow this same idea even for more groups:
dat = data.frame(group = factor(sample(1:2,200, replace = TRUE)), x = runif(200,0,3))
# creating a main effect of x for both groups
dat$y = 2*dat$x + rnorm(200,0,1)
# increasing the effect of x for group 1:
dat$y[which(dat$group == 1)] = dat$y[which(dat$group == "1")] +
lm.1 = lm(y ~ group*x, data = dat)
Exectuting this, you should find a coefficient of x near 3 (this is the coefficient for those in group 1), and a coefficient on the interaction term (group2:x) of approximately -1. The way to interpret this is that the effect of x in group 2 is 2 (3 + (-1) ).
With more groups, you'll see more interaction terms -- these relate the effect of x for each group. To determine each group effect, add the interaction to the coefficient on x.
It's important that this model assumes that the residuals are normally distributed. If this property doesn't hold in your data, you may need a different model.