1
$\begingroup$

I have a dataset and need to test if the values are statistical different to 0.33 (=1/3). Flies have the choice between 3 types of berries. So, if they do not care about the type they would lay 1/3 of their eggs in each typ of berry and that would be the null-hypothesis. I want to know, if they layed signifikcant more than 1/3 of their eggs in one specific type of berry.

With this model I would test if the values are different from 0.5:

fail<-(Summe-SummeF)
c1<-cbind(SummeF,fail)
m1<-glm(c1 ~ 1, family=quasibinomial(link= "logit"), data = data1)

But I need a test against 0.33. As far as I understand the answers to other questions regarding this topic I can change the test by including an offset command, but I quit do not understand what exactly to include to test against 0.33.

Many thanks for any help you can provide.

$\endgroup$
  • 2
    $\begingroup$ Could you clarify why you are using a quasibinomial glm for what looks like a straight test of proportions? Why not binom.test or chisq.test? If you do particularly want to use a glm for whatever reason, why a quasi model? $\endgroup$ – Glen_b Oct 21 '18 at 7:11
  • $\begingroup$ Sorry, that is the model I am using so far to test binomial values. I learned that I have to use a quasi model, if the diepersal is not within 0.5 and 1. I have to admit, aswering your question further is far beyond my state of knowledge. $\endgroup$ – R. Kienzle Oct 22 '18 at 13:31
2
$\begingroup$

Edit: This answer addresses the simple case in which there are three counts and the distribution of these counts are tested against a null hypothesis of equal distribution. Comments by the OP suggest the set-up of the experiment is more complicated than this.

Aside from the multinomial test mentioned by @a_statistician , you could also use a chi-square goodness-of-fit test. This test is probably more common, but also may be inappropriate when there are low expected counts.

observed    = c(10, 20, 30)
theoretical = c(1/3, 1/3, 1/3)

chisq.test(x = observed,
           p = theoretical)

   ### Chi-squared test for given probabilities
   ### 
   ### X-squared = 10, df = 2, p-value = 0.006738

Probably more beneficial would be looking at the multinomial confidence intervals for the proportions.

if(!require(DescTools)){install.packages("DescTools")}

library(DescTools)

observed    = c(10, 20, 30)

MultinomCI(observed,
           conf.level=0.95,
           method="sisonglaz")


   ###            est    lwr.ci    upr.ci
   ### [1,] 0.1666667 0.0500000 0.3079439
   ### [2,] 0.3333333 0.2166667 0.4746106
   ### [3,] 0.5000000 0.3833333 0.6412772

Plotting these gives us a basis for comparison.

if(!require(DescTools)){install.packages("DescTools")}
if(!require(ggplot2)){install.packages("ggplot2")}

library(DescTools)
library(ggplot2)

MCI = MultinomCI(observed,
           conf.level=0.95,
           method="sisonglaz")

MCID = as.data.frame(MCI)

MCID$Berry = c("B1", "B2", "B3")

ggplot(MCID, aes(x = Berry, y = est)) +
    geom_bar(stat = "identity",
            color = "black",
             fill  = "gray50",
             width =  0.6) +
    geom_errorbar(aes(ymin  = lwr.ci, ymax  = upr.ci),
             width = 0.2) +
    theme_bw() +
    ylab("Proportion")

enter image description here

Sources:

(Caveat: I am the author of these pages.)

http://rcompanion.org/handbook/H_03.html

http://rcompanion.org/handbook/H_02.html

$\endgroup$
  • $\begingroup$ With the chi-square goodness of fit test for multinomial distribution, I again get a simular error message as for the goodness of fit test for multinomial distribution, that x and p have to have the same number of elements. The model for looking at the multinomial confidence intervals for the proportions, is runing but the plot looks very weird and I have to admit, that this suggestion might be beyond my understanding. $\endgroup$ – R. Kienzle Oct 22 '18 at 13:58
  • $\begingroup$ Yes, x and p need to have the same number of elements, which is three in your case. Each element in x should be a single number, the count for that berry. $\endgroup$ – Sal Mangiafico Oct 22 '18 at 14:13
  • $\begingroup$ Ok and how do I include, that I have several replicates? Or do you suggest just testing the overall sum of eggs per berry within the 10 replicates? $\endgroup$ – R. Kienzle Oct 22 '18 at 14:51
  • $\begingroup$ In that case I would recommend a generalized linear model for count data, like Poisson or negative binomial, where each replicate is an observation. But it might make sense to sum them up into three counts. $\endgroup$ – Sal Mangiafico Oct 22 '18 at 15:30
  • $\begingroup$ The problem with a poisson model is, that I have in addition 4 different treatments with 1,3, 6 and 10 berries of each type. Since the flies tend to lay more eggs if they have accsses to more berries I would have to work with eggs per berry by dividing the sums through 3,6 or 10. Therefore I get data with not whole numbers and the poission model will not run. That is why I wanted to go with binomial. $\endgroup$ – R. Kienzle Oct 22 '18 at 15:46
1
$\begingroup$

If the probability of the event was 0.33, then the log-odds of the event (the linear predictor for a binomial and quasibinomial model) is:

$$\mbox{logit} 0.33 = -0.7081851$$

To fit the model in R:

data1$os <- qlogis(0.33)
m1<-glm(c1 ~ offset(os) + 1, family=quasibinomial(link= "logit"), data = data1) 

The statistical significance of the intercept term will tell you whether the risk is different from that value.

The quasibinomial works in spite of the constraint of the model probabilities, because of the quasilikelihood. But the multinomial likelihood is preferred here.

To do a multinomial model is quite easy. It's just a special case of a log-linear model:

n <- c(10, 20, 30)
p <- c(0.33)
en <- p * c(60, 60, 60)
x <- factor(0:2)
fit <- glm(n ~ x + offset(log(en)), family=poisson)

test the statistical significance of the x group indicator term.

$\endgroup$
  • 1
    $\begingroup$ This is a good solution, especially considering that the experimental design is actually more complex than indicated in the original post. It would also allow comparisons among the individual treatments. I might add the following code: library(car); Anova(fit); library(emmeans); EM = emmeans(fit, ~ x, type="response"); EM; pairs(EM) $\endgroup$ – Sal Mangiafico Oct 22 '18 at 20:14
0
$\begingroup$

Need to install package EMT.

library(EMT)

observed <- c(n1,n2,n3) ## n1, n2, n3 are observed numbers of eggs on each type of berry

prob <- c(0.33333, 0.33333, 0.33333) 

out <- multinomial.test(observed, prob)   

Example:

multinomial.test(c(10,20,30),c(0.33333,.33333,0.33333))

 Exact Multinomial Test, distance measure: p

    Events    pObs    p.value
      1891   1e-04     0.0059

It is called as goodness of fit test for multinomial distribution.

$\endgroup$
  • $\begingroup$ Thank you so much for your answer. It seams that this is exactly the test I was looking for. But now I still have a little problem, because I get an error message, when I thry to run the model. This is my input, like you suggested: observed<- c(data1$SummeF, data1$SummeC, data1$SummeD) prob<- c(0.33333, 0.33333, 0.33333) out<- multinomial.test(observed, prob) And this the error message: Error in multinomial.test(observed, prob) : Observations and probabilities must have same dimensions. Thank you for your patience! $\endgroup$ – R. Kienzle Oct 22 '18 at 13:07
  • $\begingroup$ Observed has to counts, as in the example. $\endgroup$ – Sal Mangiafico Oct 22 '18 at 14:09
  • $\begingroup$ Maybe you forgot two dollar symbols observed<- c(data1SummeF,data1SummeC, data1$SummeD). Between data1 and SummerF, add a dollar symbol. I cannot type dollar symbol. $\endgroup$ – user158565 Oct 22 '18 at 14:20
  • $\begingroup$ No, the dollar symbols are in. The just did not survived the copy-paste, cause I can not type them either. $\endgroup$ – R. Kienzle Oct 22 '18 at 14:59
  • $\begingroup$ just type in c(data1SummeF,data1SummeC, data1$SummeD) to see if you can get three numbers. (With dollars.) $\endgroup$ – user158565 Oct 22 '18 at 15:49

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.