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:

library(nlme)

# 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")]
str(dat)
'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 ~ year.f, dat, ~ 1 | state.f))
Linear mixed-effects model fit by REML
  Data: dat 
  Log-restricted-likelihood: -1231.598
  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 

Random effects:
 Formula: ~1 | state.f
        (Intercept) Residual
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))[1]; and residual variance, sigma(m0) ** 2.

So the proportion of variance accounted for the states is then:

as.numeric(VarCorr(m0))[1] / (
  as.numeric(VarCorr(m0))[1] + var(predict(m0, level = 0)) + sigma(m0) ** 2)
[1] 0.7876077

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, 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.