3
$\begingroup$

I'm playing a game with a friend where we first select a random coin from a bag. This coin has a bias $q$ chance of coming up heads and my friend knows what $q$ is when he selects the coin, and tells me "q is (some number)", but may be lying to me. I don't know what $q$ is.

I then choose heads or tails and win if I am correct.

We play this game $N$ times, each time choosing a different coin with a different $q$. However, I suspect that my friend is cheating, eg. that the value of $q$ is not the value he claims it to be.

Thus I have a sequence of values:

\begin{align} Q &= [q_1, q_2, q_3, ..., q_N] \\ R &= [T, T, H, H, ..., T] \end{align}

How can I construct a statistical test (modified version of the binomial test?) to reject the null hypothesis that my friend is not cheating from my observations?

Is this a known variant of the binomial test?

$\endgroup$
  • $\begingroup$ A lot of important information seems lacking. What is the point of "selecting a random coin" when its behavior is completely determined by its "known bias"? At what point in this process does the value of $q$ become "known" and what role does your opponent's "claim" play in the process? What would be the motivation for cheating? (That's an important thing to know, because it's possible to cheat systematically in a way that completely prevents you from detecting it: just adjust some chances up and some chances down.) $\endgroup$ – whuber Sep 21 '16 at 15:34
  • $\begingroup$ Thanks @whuber for the clarifying questions, I'll do my best: The point of the random coin is mostly setup for the problem of $Q$ and $R$, eg. I have a sequence of observed results and I'd like to know if the values of $Q$ are "correct" in some sense. There is no external motivation for cheating, in fact in the real-world problem that I tried to distill here is the $q$'s belong to a set of prior expectations and I'm trying to provide evidence that the real result of the $R$ deviate from what the expected values are. If all $q$'s were the same then this is binomial, I'm stuck when they differ. $\endgroup$ – Hooked Sep 21 '16 at 16:06
  • $\begingroup$ I made significant changes to clarify your question. Please confirm these changes are accurate, thank you! $\endgroup$ – barrycarter Sep 21 '16 at 16:14
  • $\begingroup$ Thanks @barrycarter, I think you've done a good job at improving the clarity of the question. $\endgroup$ – Hooked Sep 21 '16 at 17:04
  • $\begingroup$ The real-world question is clearer than this artificial problem! Why not just describe the situation you actually face--it doesn't sound any more complicated. $\endgroup$ – whuber Sep 21 '16 at 17:17
2
$\begingroup$

This problem is really similar to that of calibrating the output of a machine learning binary classifier. The output $Q$ is (supposedly) the vector of probabilities of getting heads, and if he is telling the truth it should be exact. So, perform logistic regression of the outcomes $R$ (1 for heads and 0 for tails) as a function of $\log(Q)-\log(1-Q)$. If your friend is telling the truth, then the intercept should be consistent with $0$ and the slope consistent with $1$. Of course, these parameters are correlated, but you can use the covariance matrix to test if the parameters are consistent with $(0,1)$, yielding a score from a $\chi^2$-distribution with two degrees of freedom.

Of course, you can also do a non-parametric test; you can consult many texts on logistic regression to find candidates (e.g. Hosmer-Lemeshow test; Pearson and le Cessie-van Houwelingen-Copas-Hosmer test, etc.). Some of these are implemented in an R package, and I'll add a mention of it when I remember its name. However, actually performing logistic regression should work just fine.

EDIT The rms package provides parametric goodness-of-fit tests for the result of a logistic regression via the residuals.lrm function. There are parametric tests done after performing logistic regression, and have some power to reject cases where the linear hypothesis isn't as good as alternative non-linear hypotheses (though they aren't perfect). Here is an example showing the le Cessie-van Houwelingen-Copas-Hosmer test:

library(rms)
dllq <- rnorm(1000,1,2)
q <- 1/(1+exp(-dllq))
y <- sapply(q,function(qi){rbinom(1,1,qi)})
fit <- lrm(y~dllq,x=TRUE,y=TRUE)
print(fit)
residuals(fit,type="gof")

The output of this may look like the following:

Logistic Regression Model

lrm(formula = y ~ q, x = TRUE, y = TRUE)

                      Model Likelihood     Discrimination    Rank Discrim.    
                         Ratio Test            Indexes          Indexes       
Obs           100    LR chi2      44.51    R2       0.493    C       0.867    
 0             36    d.f.             1    g        2.370    Dxy     0.734    
 1             64    Pr(> chi2) <0.0001    gr      10.697    gamma   0.735    
max |deriv| 1e-06                          gp       0.339    tau-a   0.342    
                                           Brier    0.145                     

          Coef    S.E.   Wald Z Pr(>|Z|)
Intercept -0.0579 0.2852 -0.20  0.8390  
q          1.0072 0.2131  4.73  <0.0001 

Sum of squared errors     Expected value|H0                    SD                     Z 
           14.4577468            14.0041517             0.3717001             1.2203253 
                    P 
            0.2223416 

So, the p-value is 22% which is more than comfortably within the "good fit" region. In combination with the compatibility of the intercept and slope with 0 and 1, this means "actual heads probability consistent with q: no evidence of lying."

I haven't yet found a "non-parametric" version of this test in R, which skips the logistic regression (since your data should already be calibrated) and which I think would be the best test in this case.

SECOND EDIT

I've had no luck finding a native R package that implements non-parametric versions of the standard binomial regression goodness-of-fit tests. However, it's not too hard to calculate them manually. There are three quantities to choose from: $$ G^2 = \sum_i R_i \log(\pi_i) + (1-R_i)\log(1-\pi_i)\\ X^2 = \sum_i \frac{(R_i-q_i)^2}{q_i(1-q_i)}\\ S = \sum_i (R_i-q_i)^2 $$ The first is the deviance, the second the Pearson statistics, and the third has a very long name I've mentioned a couple of times. Each one has been shown to asymptotically follow a normal distribution, if the probabilities $q_i$ are accurate predictors of $R_i$ (see Osius and Rojek). All you have to do is calculate their means and variances. For example, $$ E[S] = \sum_i \left( E[R_i^2] - 2 E[R_i] q_i + q_i^2 \right) = \sum_i q_i(1-q_i) \\ E[(S-E[S])^2] = \sum_i E \left[\left( (R_i-q_i)^2 - q_i(1-q_i) \right)^2 \right] + \sum_{i \neq j} 0 = \sum_i q_i(1-q_i)(1-2q_i)^2 $$ Once you figure these out, you can easily calculate them in R and get z-scores. The three z-scores you get aren't independent, but these different tests are sensitive to different types of errors in the probabilities $q_i$.

$\endgroup$
  • 1
    $\begingroup$ This is a very interesting suggestion. One concern--which you might be able to address--is that regression coefficients of $(0,1)$ do not assure the relationship actually is linear. It sounds rather like what you also need, at a minimum, is a goodness of fit test, too. That leads one to suspect that there must be some more direct and powerful way to address this situation. $\endgroup$ – whuber Sep 21 '16 at 18:07
  • 1
    $\begingroup$ @whuber Thanks for the suggestion; I've added info on some goodness-of-fit tests for logistic regression via the rms package. I'm still unsatisfied, however, because performing the logistic regression shouldn't be necessary and there is a non-parametric version of the goodness-of-fit test in rms; I've written it myself in C++ (to non-parametrically check the calibration of a logistic regression in an orthogonal dataset). However, I can't find one in R so far. $\endgroup$ – jwimberley Sep 21 '16 at 18:37
  • $\begingroup$ In lieu of a package, its simple to do the non-parametric test manually: Defining $S=\sum_i (R_i-q_i)^2$, you can show that $\mu = E[S] = \sum_i q_i(1-q_i)$ and that $\sigma^2 = E[(S-\mu)^2] = \sum_i q_i(1-q_i)(1-2q_i)^2$, which lets you calculate $z=(S-\mu)/\sigma$. The theory (originally due to Osius and Rojek, I believe) is that S is asymptotically normal, so you can turn $z$ into a p-value. $\endgroup$ – jwimberley Sep 21 '16 at 20:33

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.