3

I have the following data

Species <- c(rep('A', 47), rep('B', 23))
Value<- c(3.8711, 3.6961, 3.9984, 3.8641, 4.0863, 4.0531, 3.9164, 3.8420, 3.7023, 3.9764, 4.0504, 4.2305,
          4.1365, 4.1230, 3.9840, 3.9297, 3.9945, 4.0057, 4.2313, 3.7135, 4.3070, 3.6123, 4.0383, 3.9151,
          4.0561, 4.0430, 3.9178, 4.0980, 3.8557, 4.0766, 4.3301, 3.9102, 4.2516, 4.3453, 4.3008, 4.0020,
          3.9336, 3.5693, 4.0475, 3.8697, 4.1418, 4.0914, 4.2086, 4.1344, 4.2734, 3.6387, 2.4088, 3.8016,
          3.7439, 3.8328, 4.0293, 3.9398, 3.9104, 3.9008, 3.7805, 3.8668, 3.9254, 3.7980, 3.7766, 3.7275,
          3.8680, 3.6597, 3.7348, 3.7357, 3.9617, 3.8238, 3.8211, 3.4176, 3.7910, 4.0617)
D<-data.frame(Species,Value)

I have the two species A and B and want to find out which is the best cutoffpoint for value to determine the species.

I found the following question:

R: Determine the threshold that maximally separates two groups based on a continuous variable?

and followed the accepted answer to find the best value with the dose.p function from the MASS package. I have several similar values and it worked for them, but not for the one given above (which is also the reason why i needed to include all 70 observations here).

D$Species_b<-ifelse(D$Species=="A",0,1) 
my.glm<-glm(Species_b~Value, data = D, family = binomial)
dose.p(my.glm,p=0.5)

gives me 3.633957 as threshold:

             Dose        SE
p = 0.5: 3.633957 0.1755291

this results in 45 correct assignments. however, if I look at the data, it is obvious that this is not the best value. By trial and error I found that 3.8 gives me 50 correct assignments, which is obviously better.

Why does the function work for other values, but not for this one? Am I missing an obvious mistake? Or is there maybe a different/ better approach to solving my problem? I have several values I need to do this for, so I really do not want to just randomly test values until I find the best one.

Any help would be greatly appreciated.

OCSee
  • 35
  • 8

1 Answers1

3

I would typically use a receiver operating characteristic curve (ROC) for this type of analysis. This allows a visual and numerical assessment of how the sensitivity and specificity of your cutoff changes as you adjust your threshold. This allows you to select the optimum threshold based on when the overall accuracy is optimum. For example, using pROC:

library(pROC)

species_roc <- roc(D$Species, D$Value)

We can get a measure of how good a discriminator Value is for predicting Species by examining the area under the curve:

auc(species_roc)
#> Area under the curve: 0.778

plot(species_roc)

enter image description here

and we can find out the optimum cut-off threshold like this:

coords(species_roc, x = "best")
#>   threshold specificity sensitivity
#> 1   3.96905   0.6170213   0.9130435

We see that this threshold correctly identifies 50 cases:

table(Actual = D$Species, Predicted = c("A", "B")[1 + (D$Value < 3.96905)])
#>       Predicted
#> Actual  A  B
#>      A 29 18
#>      B  2 21
Allan Cameron
  • 147,086
  • 7
  • 49
  • 87