I would suggest a multilevel model where you let the intercept vary by brand. This is similar to your standard regression model but you permit the intercept to be different depending on the brand. It seems you are interested in seeing if the brands account for variation in the response variable. A simple approach would be to present the variance accounted for by the brands as a proportion of the total variance in the response variable. I demo with an example:
# Data prep
dat <- clubSandwich::MortalityRates[, c("year", "state", "mrate", "cause")]
dat <- dat[dat$cause == "Motor Vehicle" & dat$year %in% 1970:1975, ]
dat <- dat[, c("year", "state", "mrate")]
'data.frame': 306 obs. of 3 variables:
$ year : int 1970 1971 1972 1973 1974 1975 1970 1971 1972 1973 ...
$ state: int 1 1 1 1 1 1 2 2 2 2 ...
$ mrate: num 62.2 62.4 73.6 63.4 60.3 ...
Now we have a dataset of a response variable, mrate, and two predictors, each state (brand in OP's case) measured yearly for six years.
dat$year.f <- factor(dat$year) # Make time dummies
dat$state.f <- factor(dat$state) # Make state dummies
Next I run a multilevel model. I regress the response on the year, and let the intercept vary by state:
(m0 <- lme(mrate ~ 1 + year.f, dat, ~ 1 | state.f))
Linear mixed-effects model fit by REML
Fixed: mrate ~ year.f
(Intercept) year.f1971 year.f1972 year.f1973 year.f1974 year.f1975
66.3498221 -3.4835741 0.2247485 0.3086386 -8.1646117 -10.6781779
Formula: ~1 | state.f
StdDev: 22.29717 10.7256
Number of Observations: 306
Number of Groups: 51
Total variance in outcome comes from three sources: variance accounted for by time dummies,
var(predict(m0, level = 0)); variance accounted for by the states,
as.numeric(VarCorr(m0)); and residual variance,
sigma(m0) ** 2.
So the proportion of variance accounted for the states is then:
as.numeric(VarCorr(m0)) / (
as.numeric(VarCorr(m0)) + var(predict(m0, level = 0)) + sigma(m0) ** 2)
78.8% of the variance in the response is accounted for by the states. You can also conduct a statistical test of this variation. That is: does one improve the model by including it in the model. We run this same model using generalized least squares. Assuming a compound symmetry correlation structure is equivalent to our previous multilevel model.
mg0 <- gls(mrate ~ year.f, dat, corCompSymm(form = ~ 1 | state.f))
mg1 <- gls(mrate ~ year.f, dat) # A regular regression model ignoring states
anova(mg0, mg1) # Model comparison
Model df AIC BIC logLik Test L.Ratio p-value
mg0 1 8 2479.197 2508.827 -1231.598
mg1 2 7 2814.073 2839.999 -1400.036 1 vs 2 336.8762 <.0001
The likelihood ratio test $(\chi^2(1)=337,p<.001)$ (number in parentheses, 1, comes from subtracting the values in the df column), information criteria also suggest we retain the model where we account for the non-independence of observations within states. Here, the test would not be necessary given that states account for 79% of the variation in the response variable.
Now repeat this nine additional times. There are ways to perform all ten analyses in one analyses. But with 10 different response variables, it can be difficult to keep track of everything that's happening. And parameterizing the model may be challenging. I would suggest reading up on multilevel models if you take up this approach. So far, I have ignored heteroskedasticity or potential problems with normality.