Suppose I have 50 different brands. I am recording Variable1, Variable2,Variable3,Variable4,....,Variable10 for 50 different brands. I observed each brand 6 times. So my data will look like the following:

Brand name on the first column, Variable1 on the second and so on. I have total of 300 observations.

Goal: 1. My goal is to test if there is a difference in Variable1 between brands. I could simply use ANOVA but assumption of constant variance is not satisfied. What should I do?

I should perform goal#1 for all the variables (Variable2,....,Variable10).

What I did was I created a matrix of sample standard deviation by brand. Then I ran PCA on the matrix of sample standard deviation. I also created a matrix of sample mean by brand. After that I used cluster analysis (kmeans) on the matrix of sample mean and grouped the brands based on cluster. Now the brands falling in the same cluster will hopefully have common variance.

Now I need to test if there is difference in Variable1 between the clusters(groups)?

I also need to test if there is difference in Variable2 between the clusters.

I need to do this upto Variable10.

  • 1
    $\begingroup$ Why are you creating the clusters? Because of heteroskedasticity? $\endgroup$ Jul 15, 2018 at 1:05
  • $\begingroup$ Yes heteroskedasticity is the reason that I am creating clusters. $\endgroup$
    – Hello
    Jul 15, 2018 at 2:57
  • 1
    $\begingroup$ To calculate the standard deviation of a brand, you used the six time points on variable 1? $\endgroup$ Jul 15, 2018 at 11:16
  • $\begingroup$ Yes I only have six observations per brand. $\endgroup$
    – Hello
    Jul 15, 2018 at 13:15

1 Answer 1


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
  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 $(\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.

  • $\begingroup$ Thanks for the demo. Is this demo similar to using proc mixed in SAS? $\endgroup$
    – Hello
    Jul 16, 2018 at 4:14
  • $\begingroup$ Never used SAS. Maybe. Mixed sounds like mixed effects which is same as multilevel. $\endgroup$ Jul 16, 2018 at 4:15
  • $\begingroup$ Okay. Would it be feasible to run multilevel model when I have 50 different brands? Which output tells me if there is difference in Variable1 between brands? $\endgroup$
    – Hello
    Jul 16, 2018 at 4:23
  • $\begingroup$ That's what I have in my example, 50 states, so yes. I reported multiple outcomes, proportion of variance explained by states. And likelihood ratio rest comparing model without states to model with states. I've now emphasized the test with p value. $\endgroup$ Jul 16, 2018 at 11:43
  • 1
    $\begingroup$ There are fifty states. You're talking comparing all of them to each other? If you have a specific hypothesis about how the states are different, you can test that hypothesis using this method. If you don't, there are methods to identify state clusters given their data, like your original cluster analysis, or latent class analysis methods. I don't know them well enough to comment. The one additional thing you can obtain from this model is descriptive measures of the states means then you can identify the most extreme states on the outcome, and the average ones. $\endgroup$ Jul 16, 2018 at 16:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.