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

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.

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

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.

• Thanks for the demo. Is this demo similar to using proc mixed in SAS? Jul 16, 2018 at 4:14
• Never used SAS. Maybe. Mixed sounds like mixed effects which is same as multilevel. Jul 16, 2018 at 4:15
• 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? Jul 16, 2018 at 4:23
• 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. Jul 16, 2018 at 11:43
• 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. Jul 16, 2018 at 16:35