Equality of means between groups when variances within the same group is not constant

Question

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.

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:

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