Account for multiple treatments to the same subject when the treatments are a numerical variable

Question

254 individuals were asked their purchase intent of a product (on a 1-5 scaled) at different prices \$699, \$799, \$899, \$999, and \$1099. The data looks like:

> head(price.responses)
  699 799 899 999 1099
1   4   4   3   3    2
2   5   5   4   4    5
3   4   4   4   5    5
4   1   2   3   4    5
5   1   1   3   5    5
6   5   3   3   2    1

I'm transforming the data like so:

people = rep(sapply(1:254, function(x) { paste("person", x, sep="")}), 5)
s.price.responses = stack(price.responses)
colnames(s.price.responses) = c("purchase_intent", "price")
s.price.responses$price = as.numeric(as.character(s.price.responses$price))
s.price.responses = cbind(s.price.responses, people)

> head(s.price.responses)
  purchase_intent price  people
1               4   699 person1
2               5   699 person2
3               4   699 person3
4               1   699 person4
5               1   699 person5
6               5   699 person6

A simple chi-squared test shows that the distributions of purchase intent is not the same across the the different price points

chisq.test(table(s.price.responses[,c("price", "purchase_intent")]))

        Pearson's Chi-squared test

data:  table(s.price.responses[, c("price", "purchase_intent")])
X-squared = 194.2906, df = 16, p-value < 2.2e-16

I'd like to do something like a repeated measures anova test which accounts for the fact that the same person was put through multiple samples (\$699, \$799,....) and a linear regression which would give an estimate of the relationship between purchase_intent ~ price. Right now I'm doing:

price.regression = lm(purchase_intent ~ price + price:people, data=s.price.responses)

Am I accomplishing what I what with this model? I.e. am I getting an estimate of the effect of price on purchase_intent once accounting for the inherent variation among people? Are there better methods for achieving this?

Answer 1

The straightforward answer is, while your model may be giving you the proper estimate for the price coefficient, it is not giving you the proper standard errors, degrees of freedom, and tests of significance. This is because there are two levels of error in your data: person-level error (e.g., certain people may have more strong purchasing intentions for products at any price; certain people may react more to price than others) and observation-level error (e.g., a person may have been distracted at the very end of the experiment and thus may have been less thoughtful about his/her responses at that time).

The proper way to model data of this structure is a mixed model or linear mixed effects model, which is a generalization of the linear regression model that allows more complex error structures than the classic linear model. Linear mixed effects models divide effects into fixed effects and random effects; fixed effects are assumed to be the same for all units within your data, whereas random effects are allowed to vary across units.

One of the classic applications of the linear mixed effects model is to data where you have multiple observations per person, as in your case. The "units" across which some of your effects can vary are people; as I mentioned above, it could be that different people have baseline purchase intentions that are different from others (i.e., the intercepts for some people may be different than for others) and different people may react differently to price than others (i.e., the slopes for price for some people may be different than for others).

Linear mixed effects models is an expansive and rapidly evolving topic, so I highly recommend doing some reading if you wish to apply this method to your data. Fortunately, there are a wide variety of resources available (including many on this site). You can also find some good annotated resources using the lme4 package in R, which is one of the standard packages for fitting linear mixed effects models. . Just bear in mind that, as with most statistical techniques, the technique itself is useless without the knowledge of how to interpret the results of the technique.

Finally, I have written some R code that simulates some data that might be of the same structure as your own and shows how to fit a linear mixed effects model to those data.

require(lme4)
require(plyr)

set.seed(34)

id <- rep(1:254, 5)
id <- id[order(id)]
price <- rep(c(699, 799, 899, 999, 1099), 254)
d <- data.frame(id, price)

# Fixed effect of price
d$purchase_intent <- -.01 * price + 11.99

# Add in random effects at the person level
d <- ddply(d, "id", mutate,
           id_int = rnorm(1, sd = 1), # Random intercept for each person
           id_price = rnorm(1, sd = .01), # Random slope for each person
           purchase_intent = purchase_intent + id_int + id_price * price)

# Add in observation-level error
d$purchase_intent <- d$purchase_intent + rnorm(254 * 5, sd = 5)

# Rescale purchase_intent to be between 1 and 5
d$purchase_intent <- 5 * (d$purchase_intent -min(d$purchase_intent))/(max(d$purchase_intent) - min(d$purchase_intent))

# Fit the model. The parenthesis indicates the random effects for id.
# The (price|id) indicates that we are adding a random slope for each
# unit defined by id.  lmer also automatically adds a random intercept
# for each unit defined by id
mod <- lmer(purchase_intent ~ price + (price|id), data = d)
summary(mod)

Linear mixed model fit by REML ['lmerMod']
Formula: purchase_intent ~ I(price/100) + (I(price/100) | id)
   Data: d

REML criterion at convergence: 1140.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.3888 -0.5980  0.0045  0.6066  3.0172 

Random effects:
 Groups   Name         Variance Std.Dev. Corr 
 id       (Intercept)  0.046593 0.21585       
          I(price/100) 0.003676 0.06063  -0.20
 Residual              0.076685 0.27692       
Number of obs: 1270, groups:  id, 254

Fixed effects:
              Estimate Std. Error t value
(Intercept)   3.410020   0.051806   65.82
I(price/100) -0.049844   0.006683   -7.46

Correlation of Fixed Effects:
            (Intr)
I(pric/100) -0.814

---

whuber is right that you may want to consider treating purchase_intent as an ordinal variable. In psychology (my field), most reviewers would allow you to treat purchase_intent as an interval variable as long as each unit in the variable is well-anchored and as long as the instructions were clear to the participants. However, if you are not comfortable with the interval assumption, you can use a cumulative link mixed model using the ordinal package.

Using the clmm function, the function from ordinal that fits mixed models, is relatively straightforward -- the clmm function uses very similar syntax as the lmer function. The main limitation to clmm is that, at present, clmm cannot accommodate random slopes.

See the code chunk below for an example of how to use clmm. This code chunk picks up with the dataset that I created in the previous code chunk.

require(ordinal)

# Make purchase_intent an ordered factor
d$purchase_intent <- ordered(round(d$purchase_intent, digits = 0))

mod <- clmm(purchase_intent ~ I(price/100) + (1 | id), data = d)
summary(mod)

Cumulative Link Mixed Model fitted with the Laplace approximation

formula: purchase_intent ~ I(price/100) + (1 | id)
data:    d

 link  threshold nobs logLik  AIC     niter     max.grad cond.H 
 logit flexible  1270 -970.87 1955.75 335(2599) 5.60e-04 4.9e+03

Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 10.98    3.313   
Number of groups:  id 254 

Coefficients:
             Estimate Std. Error z value      Pr(>|z|)    
I(price/100) -0.30418    0.05116  -5.945 0.00000000276 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:
    Estimate Std. Error z value
0|1 -14.7632     1.3767 -10.724
1|2 -11.3372     0.7626 -14.866
2|3  -5.5937     0.5599  -9.991
3|4   0.4206     0.5042   0.834
4|5   7.9041     0.9802   8.064