4
$\begingroup$

I am relatively new to R, and I cannot figure out an acceptable statistical test to determine whether the number of dead aphids in a population depends on the starting number of aphids and the time. I have four columns representing the number of dead aphids in a population of aphids in separate Petri dishes.

How do I decide on a test, and what commands can I use to produce these results and a plot displaying the differences between treatments and time points?

Here is the dput of the dataset for reproduction:

structure(list(result_after_time_hours = structure(c(1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), levels = c("2", 
"18"), class = "factor"), four_aphids_dead = c(3, 2, 1, 2, 2, 
1, 1, 3, 1, 3, 2, 3, 2, 3), eight_aphids_dead = c(2, 2, 3, 1, 
2, 3, 1, 4, 3, 4, 8, 5, 2, 4), sixteen_aphids_dead = c(3, 3, 
5, 3, 2, 3, 3, 8, 10, 6, 8, 11, 8, 10), thirtytwo_aphids_dead = c(8, 
4, 8, 5, 6, 10, 6, 19, 18, 22, 20, 13, 15, 11)), row.names = c(NA, 
-14L), class = c("tbl_df", "tbl", "data.frame"))

I've tried following online flowcharts for tests, but I cannot figure out how to follow them well, and the basic test do not seem to work as expected.

I'd like to answer the question of whether the number of aphids in a population affects the feeding behavior of lacewing larvae based on the proportion of what's eaten, with a hypothesis of "lacewing larvae will eat a higher proportion of aphids in a larger population compared to a smaller population."

$\endgroup$
5
  • 1
    $\begingroup$ Insect person here ! What do the 4 columns represent ? You want to test if the difference in mortality of aphids of different populations is influenced by what ? $\endgroup$
    – CaroZ
    Feb 11 at 19:19
  • $\begingroup$ each column represents a different petri dish that contains 4, 8,16 or 32 aphids and a single lacewing larva. The values recorded are the number of aphids from this original population that have been killed by the lacewing after a 2 hour period and then an 18 hour period. I'd like to test whether there is a statistical difference between the amount of aphids the lacewing killed depending on how many aphids were in the petri dish, and then seeing if there is a difference between these results after 2 or 18 hours. $\endgroup$
    – A-okay
    Feb 11 at 19:24
  • $\begingroup$ How many populations do you have in total ? $\endgroup$
    – CaroZ
    Feb 11 at 20:44
  • 3
    $\begingroup$ It's difficult to suggest any statistical test when you seem undecided on the research question you want answer. But besides, your data remind me of McElreath's tadpole example, so that might be a good place to start. $\endgroup$
    – Durden
    Feb 11 at 20:47
  • $\begingroup$ each of the rows is a separate population, to act as replicates. 2 and 18 hours are the same populations as each other. so 7 sets of 4 aphids, 7 sets of 8 aphids etc. $\endgroup$
    – A-okay
    Feb 11 at 20:49

2 Answers 2

8
$\begingroup$

I would organize the dataset differently.

The easiest would be to have one line per Petri dish and per time point.

Each Petri dish should have a unique identifier, meaning this identifier will come back 2 times in your ID column, since you have 2 time points per Petri dish, meaning that each Petri dish would have 2 lines, one for the number of dead aphids at 2 hours, one line for the number of dead aphids at 18 hours.

Then there should be a time column which contains the information about the time point (2 or 18 hours).

Then finally, one column with the number of dead aphids, and one column with the number of live aphids. And finally, a column containing the starting number of aphids there were in the dish.

The Petri dish should be a random factor in a GLM for a binomial distribution, since it is repeated twice (2 time points). Then it is possible to analyse the proportion of dead aphids with cbind in the glmer function of the lme4 package. It should look a bit like :

glmer(cbind(number_of_dead_aphids, number_of_live_aphids) ~ 
    time*starting_number_of_aphids + (1|Petri_dish)), 
    data=your_data,family=binomial

Edit: to go further, you should then try to know whether the interaction is significant, and which effects are in general. For this you could build all possible models:

glmer(cbind(number_of_dead_aphids, number_of_live_aphids) ~ 
    time+starting_number_of_aphids + (1|Petri_dish)), 
    data=your_data,family=binomial

glmer(cbind(number_of_dead_aphids, number_of_live_aphids) ~ 
    time + (1|Petri_dish)), 
    data=your_data,family=binomial

glmer(cbind(number_of_dead_aphids, number_of_live_aphids) ~ 
    starting_number_of_aphids + (1|Petri_dish)), 
    data=your_data,family=binomial

glmer(cbind(number_of_dead_aphids, number_of_live_aphids) ~ 
    1 + (1|Petri_dish)), 
    data=your_data,family=binomial

The last one is the null model. You could then compare the AICs of all these models and chose the one with the lowest AIC/the simplest one with the lowest AIC in a delta AIC of 2.

You could use the effects package to plot(allEffects(your_model)) in order to have an idea of what the graphical representation would look like.

$\endgroup$
2
  • 1
    $\begingroup$ how do I go about analysing what the results from this glmer indicate? and is there a way to visualize this in a graph? $\endgroup$
    – A-okay
    Feb 11 at 21:59
  • $\begingroup$ For the analysis of the results, you should research it a bit. You should first try to know whether the interaction time*starting number of aphids is significant. If yes, it means the dynamics of aphids consumption is different according to the starting number of aphids. To visualize the results, you could first use the effects package and plot(allEffects(model)) $\endgroup$
    – CaroZ
    Feb 12 at 22:01
4
$\begingroup$

Before one thinks about "tests", try drawing a graph.

It's pretty obvious that more aphids are dead in the dishes with more aphids at the start, so I'd look at the proportion of aphids that were dead in each dish. And I'd take the average number of aphids over each experiment.

That gives me these numbers for the proportion dead in each group of experiments

time 4 aphid 8 aphid 16 aphid 32 aphid
2 hr 0.43 0.25 0.20 0.21
18hr 0.61 0.54 0.54 0.53

Then I drew a bar chart. Excel is good at drawing them, or you can do it by hand.

enter image description here

You don't need tests and p-values to immediately notice that the bars for 8,16, and 32 are pretty much the same. It seems clear that the lacewing kills about 1/5 to 1/4 of the aphids in 2 hours and about 1/2 the aphids in 18 hours, regardless of how many there are. It seems that it kills a slightly greater proportion when the supply is low (but the numbers are so small that it is unlikely to be significant)

Doubtless R could do this same graph, but the benefit of R is when you have lots of data, because computers are fast. The danger of R is that you can throw "tests" at your data without thinking. Slowing down and graphing the data can give you time to think about what you are actually interested in, and then whether doing a formal hypothesis test is necessary.

It can also slow you down, and notice any inconsistencies in the data. I noticed that in the second experiment with 4 aphids, there were 2 dead aphids after 2 hours, but only 1 dead aphid after 18 ... something odd!

$\endgroup$
3
  • 4
    $\begingroup$ Definitely good to graph the data before doing a statistical test (see the famous gorilla study biorxiv.org/content/10.1101/2020.07.30.228916v1). However the benefit of R is not just for crunching large datasets. It also produces a script that you can use to reproduce your own data viz/stat analysis later on, or share with others, or apply to a different dataset with the same structure. $\endgroup$
    – qdread
    Feb 12 at 16:45
  • $\begingroup$ @qdread agreed, but I think in this case (and some other questions I've seen here) the use of a test in software is done instead of thinking about the data, not in addition. In this case in particular, the graph quickly answers "the question of whether the number of aphids in a population affects the feeding behavior of lacewing larvae based on the proportion of what's eaten". My answer attempts to emphasise the importance of considering the data, and not reducing everything to a "p" - at least, not in the first instance. $\endgroup$
    – James K
    Feb 14 at 10:33
  • 1
    $\begingroup$ Yes I agree. As a statistician I have to constantly try to convince the scientists I work with not to fetishize statistical tests. But I do see the suggestion of a trend in your barplot at the 2 hour time point, where >40% of aphids are eaten in the 4-aphid treatment but 20% or less in the 8, 16, and 32 treatments. That is potentially a meaningful effect which might be worth testing. $\endgroup$
    – qdread
    Feb 14 at 15:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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