5
$\begingroup$

Upfront

I am not a statistician but a medical doctor. I have a working knowledge of statistical methods in my field, but this is my first time with pairwise comparisons and due to a lack of formal math training I am having a hard time to get into this new field. This is one of the questions with a pledge for simple worded answers.

Object of research

I have 8 dietary supplements and want to know, if they differ in taste and if so, which is best. So I asked 37 people to do pairwise comparisons of random pairs. I got 224 observations/votes.

Methods

Internet research suggests, that this should be evaluated via Bradley-Terry analysis and I found the BradleyTerry2 package for R and followed their vignette. The result of a simple call to the BTm modelling function is this:

Call:
BTm(outcome = adm$winner == "A", player1 = adm$A, player2 = adm$B)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.3883  -0.6068   0.2195   0.7277   2.5585  

Coefficients:
    Estimate Std. Error z value Pr(>|z|)    
..2  -2.2017     0.5024  -4.382 1.17e-05 ***
..3  -1.3414     0.4659  -2.879 0.003991 ** 
..4   1.4892     0.5862   2.540 0.011070 *  
..5   1.1937     0.5293   2.255 0.024112 *  
..6  -1.5988     0.4932  -3.241 0.001189 ** 
..7  -1.7452     0.4896  -3.564 0.000365 ***
..8  -2.5201     0.5231  -4.817 1.46e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 310.53  on 224  degrees of freedom
Residual deviance: 185.07  on 217  degrees of freedom
AIC: 199.07

Number of Fisher Scoring iterations: 5

So the function chose "1" to be the reference category and fixed the ability of that to 0 and positioned all other supplements with estimates and standard errors of estimates around that. As all of these estimates are significantly different from 0, I can conclude, that these supplements vary significantly in taste.

Now the next step (chapter 4, page 11, in https://cran.r-project.org/web/packages/BradleyTerry2/vignettes/BradleyTerry.pdf ), according to the package's vignette, is to calculate "quasi variances", which leads to the following table:

   estimate        SE   quasiSE   quasiVar
1  0.000000 0.0000000 0.3657419 0.13376715
2 -2.201721 0.5023954 0.3183952 0.10137549
3 -1.341350 0.4659285 0.2914081 0.08491867
4  1.489221 0.5861939 0.5199500 0.27034804
5  1.193664 0.5292600 0.4656436 0.21682395
6 -1.598792 0.4932263 0.3113142 0.09691653
7 -1.745216 0.4896299 0.3092022 0.09560602
8 -2.520092 0.5231470 0.3418877 0.11688718

I am having a hard time trying to understand, what is happening here and why we do this. Apparently, "1" now has a standard error, too, and all other standard errors have become a little smaller. This allows me to draw a plot of the estimates with 1.96 times their quasi standard error:

enter image description here

In the original model, the estimate of "4" was significantly different from 0, but now the error bars of "1" and "4" do overlap a lot. This appears to be a contradiction: Are the estimates of "1" and "4" significantly different, or are they not?

Specific questions

  1. How do I properly do an omnibus test on the Bradley-Terry model to proove, that it performs significantly better then a null model, i. e. that my supplements differ in taste?
  2. Is the original model or the "quasi-variances-thing" the best way, to present my model and how can I do inference on whether the differences between two arbitrary supplements are significant?
  3. Please give or point to an explanation of why (not how) we compute quasi variances in this context in as simple terms as possible.
Improve this question
$\endgroup$
1
+50
$\begingroup$

I think what you have done so far is fine; your decision to look into the quasi-variances is correct. We use the quasi-variances as that we care for is the difference between the estimated coefficients. To quote Firth (Sociol. Methodol., 2003): "Quasi standard errors are potentially useful in any situation where contrasts among parameters $\beta_1, \dots, \beta_K$ are of primary or exclusive interest." The use of a Bradley-Terry model falls squarely in this category as the coefficients $\beta$ are effectively useless without a reference level.

The problem: While we can readily estimate the magnitude of the difference $\hat{\beta}_j - \hat{\beta}_k$, getting the associated standard error of this difference is often tricky. Just using the standard idea: $Var(\hat{\beta}_j - \hat{\beta}_k) = (Var(\hat{\beta}_j) + Var(\hat{\beta}_K))^{\frac{1}{2}}$ might be wildly off because these are not uncorrelated estimates; the proper term is: $(Var(\hat{\beta}_j) + Var(\hat{\beta}_K) -2Cov(\hat{\beta}_j, \hat{\beta}_k))^{\frac{1}{2}}$. If we ignore the covariance term, our visualisation would also be problematic as the associated "comparison intervals" would assume uncorrelated samples. Finally even if we assume uncorrelated samples, based on our available information in the lm output we would not be able to estimate the variance associated with the "reference" level.

The solution: Calculate some statistics, the quasi-variances, that have the property to behave like the standard errors of $\hat{\beta}_j$ and $\hat{\beta}_k$ are uncorrelated. To estimate the quasi-variances, we use the estimated variance-covariance matrix. In particular, the main idea is to formulate a generalised linear model $M$ with exponential link where the response $y$ is the true standard error ($Var(\hat{\beta}_j - \hat{\beta}_k)$), the covariates $X$ correspond to the variance contrasts between the levels of the factor we are interested in and MLE estimates $\theta$ (the $\beta$'s in this model $M$) are the quasi-variances. Given these $\theta$s we can immediately get the associated standard error for the difference $\hat{\beta}_j - \hat{\beta}_k$ ignoring any covariance terms by using $( \theta_j + \theta_k)^{\frac{1}{2}}$ as the approximate standard error.

OK, to put some numbers on all of these. Reminder:

   estimate        SE   quasiSE   quasiVar
1  0.000000 0.0000000 0.3657419 0.13376715
...
4  1.489221 0.5861939 0.5199500 0.27034804

The difference $\hat{\beta}_1 - \hat{\beta}_4$ is $-1.489$. If we just ignore the covariance and use the standard errors of these estimates we get ... nothing comprehensible because $\beta_1$ is known to be $0$. We could argue that the standard error is $(0 + 0.586^2)^{0.5} \approx 0.586$; to do that we would have to assume that the variance of $\beta_1 =0$ is a plausible thing and that $\beta_1$ does not covary with $\beta_4$. Both assumptions seem implausible so let's move on. Using the quasi-variances we have $(0.13376715 + 0.27034804)^{0.5}$ $\approx$ $0.636$ as our quasi-standard error. This suggests that it is plausible that $\hat{\beta}_1 - \hat{\beta}_4$ is significantly different than $0$.

David Firth has two accessible (and aptly named) papers on the matter:

  1. Overcoming the Reference Category Problem in the Presentation of Statistical Models in Sociological Methodology (2003) and
  2. Quasi-variances (with De Menezes) in Biometrika (2004).

I would suggest you look at both papers as they are generally easily to follow and cover all the relevant details. I would start with the Sociol. Methodol. paper first and use the Biometrika to fill in any gaps in my understanding. Finally, regarding your particular sub-questions:

  1. As mentioned, what you did seems fine. Based on the graph and the SE of the differences, the two supplements are probably different in a statistically significant way. Once more, notice that the parameter $\beta_1$ is 0 by convention for the purposes of estimation; try setting $\beta_8$ as your reference level and see what happens. In any case, go ahead and calculate the quasi-standard errors to convince yourself and your readers. :)
  2. Yes, the "quasi-variances-thing" is indeed the best way to infer if the difference between the two supplements is significant. Do not get discouraged immediately if the "comparison intervals" slightly overlap, there may still be a statistically significant difference between the means of the two samples (as here).
  3. Basic exposition and references provided in the main body of the post.
Improve this answer
$\endgroup$
  • $\begingroup$ Thank you very much, that is a great answer. Short additional question: As for an omnibus test: I would like to present a $p$-value, as medical audiences are traditionally focused on $p$-values. Can I use the above given result Null deviance: 310.53 on 224 df, Residual deviance: 185.07 on 217 df for a likelihood ratio test as p = 1 - pchisq(310.5-185.1, 224-217)? (Found this for logistic regression) $\endgroup$ – Bernhard Jun 24 '18 at 12:03
  • $\begingroup$ I am glad I could help! Disclaimer: I would pick bootstrap over information criteria (eg. AIC(c)) and I would pick information criteria over LRT in most cases. That said, yes, your proposal is fine. In any case, it is clear that your current model is better than the null model. If you wish, you can try defining a "pseudo" null model (e.g.BTm(outcome = adm$winner == "A", player1 = adm$A, player2 = adm$A)) and then use an LRT directly too (e.g. epiDisplay::lrtest(nullModel, usedModel)). (the "pseudo"-null model has residual deviance that equates the null deviance) $\endgroup$ – usεr11852 Jun 24 '18 at 15:31
  • $\begingroup$ Thank you very much, with your answer and the links I feel confident enough that my analysis is sound. Accept and bounty well earned. I am curious though, how Bootstrapping and how AIC lead to inference of the model. Is there a short answer or a link answer or should I open another question for that? $\endgroup$ – Bernhard Jun 24 '18 at 18:56
  • $\begingroup$ Great, than you! I read these papers years back and then... forgot about them! :D Model-based inference is a huge discussion but in general I prefer to bootstrapping instead of being based on the assumption of a parametric model (e.g. a $\chi^2$ distribution behind the statistic of interest) to get certain $p$-values. Similarly I would pick an information criterion to say: "Model A is more plausible than model B" as a second option. Unfortunately I don't have one link to point you too... Check the book Information Criteria and Statistical Modeling by Konishi & Kitagawa for starters! :) $\endgroup$ – usεr11852 Jun 24 '18 at 20:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.