Factorial design for interactions with ordinal DV I want to make a 4 question, multiple choice survey in which each question asks about an analogous range of actions for slightly different scenarios.  Each participant would only be administered the survey once.  For example:  
"In situation W how would you behave? (5[very aggressively] -- 1 [not aggressively])."
"In situation X how would you behave? (5[very aggressively] -- 1 [not aggressively])."
"In situation Y how would you behave? (5[very aggressively] -- 1 [not aggressively])."
"In situation Z how would you behave? (5[very aggressively] -- 1 [not aggressively])."
I then wanted to sum the responses to the 4 questions, and use that total score as my DV.
I plan to have 2 IVs.  
The first will be relationship status, and consist of 7 levels (e.g. Married, single, divorced, dating, widowed, etc).
The second will be, say, partner hair-color preference, and will consist of 3 levels (e.g. blond, brunet, redhead).
So, my understanding is that this will be a 3x7 factorial design, with a total number of 21 individual conditions, and an ordinal DV.  
Ideally, I'd like to do a 3x7 factorial ANOVA to test for main effects and interaction effects, but I don't think my summed, ordinal survey response DV is amenable to this.  Is there some other analysis i could use in this situation?  Would the Kruskal-Wallis be appropriate?  Thanks!
 A: 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.
set.seed(123)

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.
library(reshape2)

Wide.Data <- dcast(Data, ID ~ Situation, value.var="Value")[, c("W", "X", "Y", "Z")]

# Get the correlation matrix.
cor(Wide.Data)

# 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)
summary(pc)
biplot(pc)
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.)
library(nlme)

fit <- gls(
  Value ~ Relationship * Hair * Situation,
  correlation = corSymm(form = ~ 1 | ID),
  weights = varIdent(form = ~ 1 | Situation),
  data=Data
)

anova(fit, type="marginal")

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,
  data=Data
)

anova(fit, type="marginal")

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.
library(plyr)

Mean.Data <- ddply(
  Data,
  .(ID, Relationship, Hair),
  summarise,
  Mean=mean(Value)
)

# Fit using the usual machinery.
library(car)

fit <- lm(Mean ~ Relationship * Hair, data=Mean.Data)

Anova(fit)
plot(fit)
summary(fit)

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