The Kruskal-Wallis test is not appropriate in this case, because it deals with a design that is similar in structure to a one-way analysis of variance (ANOVA). Also, even if it were appropriate, it does not really deal with the theoretical problems of the distribution of these data any better than ANOVA might be expected.
Linear models can be thought of as a first-order approximation to the underlying phenomenon that you are measuring. Pragmatically, a primary question is "How well does the model fit the data?" The tests for statistical significance are reasonably robust under some deviations from the underlying assumptions.
In this case, there could be an issue if answers to some of the scenarios were negatively correlated. Averaging (or summing) the answers from each subject kind of implicitly assumes that they all measure the same thing. That is very reasonable in this case, but it is worth checking out.
A linear model is going to give you more bang for the buck in terms of being able to accurately model your design, quantify things, and slice and dice the results. So, if a linear model is a reasonable approximation then I would use it. There are other modeling methods that could be used here --- it would be comforting to see that they agree with the linear model chosen.
Here are some approaches that you might use in
R. First, I am going to set up some totally random fake data.
# Make up some data.
Relationship <- c("married", "single", "divorced", "dating", "widowed", "other1", "other2")
Hair <- c("blonde", "brunette", "redhead")
Situation <- c("W", "X", "Y", "Z")
Data <- data.frame(
ID = rep(1:100, each=4),
Relationship = rep(
sample(Relationship, size=100, replace=TRUE, prob=rep(1/7, 7)), each=4),
Hair = rep(sample(Hair, size=100, replace=TRUE, prob=rep(1/3, 3)), each=4),
Situation = rep(Situation, 100),
Value = sample(1:5, size=400, replace=TRUE, prob=rep(1/5, 5))
Next, look at the pattern of results:
# Plot the data to look for patterns in the response values.
Wide.Data <- dcast(Data, ID ~ Situation, value.var="Value")[, c("W", "X", "Y", "Z")]
# Get the correlation matrix.
# Try a scatterplot matrix. Doesn't work too well. It would
# be nice to have a matrix mosaic plot. Still...
pairs(sapply(Wide.Data, jitter, amount=0.25))
# Look at the principal components.
pc <- princomp(Wide.Data)
loadings(pc) ### To look at the simple factor loadings.
pc$scores ### If you wanted to use the scores.
# Tabulate each pair. May need to stare at the tables a while.
with(Wide.Data, table(W, X))
# ... Do the rest of the pairs the same way.
Here the goal is to explore for obvious relationships among the answers. There are many other plots and summaries that would be good to generate. The principal components analysis could give you an idea of whether the dimensionality of the questions is going to be well represented by a one dimensional summary. If the loadings on the first component are roughly equal and the variance explained by the first component is high then an average is well supported.
Next, look at linear model fits. A model that accurately represents the design would have errors in responses from each subject arbitrarily correlated with differing variabilities:
# Fit an unstructured correlation matrix for measurements from the
# same subject. (Note that this depends critically on the sort order
# of the Situation items being the same for each subject. The safe way
# is to define a "Time" variable that gives correct order.)
fit <- gls(
Value ~ Relationship * Hair * Situation,
correlation = corSymm(form = ~ 1 | ID),
weights = varIdent(form = ~ 1 | Situation),
This model is pretty rich in parameters and you could expect problems fitting it. Also, it is modeling at the granularity of the individual response value, so the linear model approximation might be tenuous. But, maybe with enough data it would not be bad.
A simpler model would provide a simple random effect for each subject, implying all errors in response values from each subject are equally (and positively) correlated:
# Fit a simpler model with a random effect for each subject.
fit <- lme(
Value ~ Relationship * Hair * Situation,
random = ~ 1 | ID,
This model is pretty likely to work but again may be too granular to fit well. Finally, you can fit a model to some summary of the four response values (such as the average, or some other combination):
# Compare with results from simply averaging the values.
Mean.Data <- ddply(
.(ID, Relationship, Hair),
# Fit using the usual machinery.
fit <- lm(Mean ~ Relationship * Hair, data=Mean.Data)
# Plot effects. Pair-wise contrasts. Et cetera. Et cetera.
One of these last two approaches is likely to be pretty good, if you have a decent amount of data. Interpreting the results is going to be dependent on your population and sampling strategy.
Other options for fitting the data:
- Carry out the analysis for each question separately (better to average if appropriate)
- Model the ranks using modern approaches such as nparLD or rlme
- Dichotomize the responses and fit binomial regressions
- Fit a multinomial mixed effect model (good luck, but see this post)
These approaches tend to lose information in varying degrees, which changes the type of conclusion you can muster after all is said and done.