If you would be willing to send data (or an example/subset of data) in attachment, I could try to understand what you did and maybe give some advise.

Using comparative tests is appropriate.

Tow or One - Way ANOVA Ordinary or RM

Thank you all so much for your kind replies!

Jos Feys I attach a picture where I show you the kind of data I'm currently dealing with. I've invented some data and wrote it on a table in exactly the same format I'm using for my actual data.

As you can see, for each silique, I obtain the total number of seeds if do "aborted" + "normal". The thing I was considering is, in R, to create a binomial matrix with two columns, (success = normal, failures = aborted) and analyze this with whatever model that fit properly.

So, you agree this would be correctly analyzed if I used the binomial distribution? I've been reading to check also which test should I use. As Philip D DeWitt was saying, I was reading into those tests and some more to come to a conclusion (since I've never done anything like that!). The only thing I'm familiar with is ANOVA and Linear Regression.

The thing is we want also to report this as %. I understand binomial is also the proper distribution for this, right?

Again, thank you all so much for your replies.

Nico.

In attachment, an R script for the analysis I would suggest for the data you gave (see above). I suggest you try it piece by piece.

Jos Feys Thank you so much for the time you dedicated to my issues. I've run the script you provided me in attachment. I believe I was able to follow it (mostly) and learned about some new packages/plots/commands I knew nothing about! So, first of all, thank you for that.

Regarding my data set, I believe this analysis may not be the best fit, probably not due to execution (obviusly, I absolutely trust your knowledge) but because, in our field, transforming data (logarithm) or using non - parametric tests are not apreciated for reviewers. Also, I strongly believe that I must use a test that considers random effects, since my data is not independent (although I believe this is not entirely necessary for no - parametric, but I'm not so sure). Lastly, if I understand correctly, in your analysis you consider siliques (fruits) as repeatead measures (that's what I understand both from the script and the plot) and that would not be my case, since siliques (the units where I count aborted seeds) are all different in all my plants (at random, you might say).

I would love if anyone else could comment on what I'm saying, since I'm kind of new in plant molecular biology. I believe (and I might be wrong) that the answer I'm looking for has more to do with what Philip D DeWitt was saying earlier, but I will keep on looking.

Jos, again, I can't thank you enough.

Nico.

Indeed, transformations on data are not recommended (to say the least). Amongst many others: Glass, Peckham and Sanders (1972), in a very extensive review, have argued against it.

Ref.:

Glass, G.V., Peckham, P.D., & Sanders, J.R. (1972). Consequences of failure to meet assumptions underlying fixed effects analyses of variance and covariance. Review of Educational Research , 42(3), 237-288.

I used the Silique variable as repeated measures (a Within factor) because these measures are dependent and crossed with the plants. If this would be a Between factor, you would need much more plants. In any case, it would not be possible to analyse the data in the prueba example with a Between x Between factorial design anova. I add the new script in attachment, in which I include the non-parametric analyses 'just for fun' or just the consolidate the parametric anova.

An overdispersed binomial model is the most appropriate choice here, e.g. a mixed effect binomial or a beta-binomial model. Repeated measures ANOVA or beta regression on collapsed proportions are applicable too, and have the advantage of not requiring inclusion of a correlation structure (if you collapse data from each plant the between plant variation is accounted for at the observation level), but they are misspecified models for this kind of data.

In R you can fit a mixed effect logistic regression using glmer from lme4 package, and a beta-binomial model using betabin from aod package. If you chose a mixed effect logistic, you may need to include a slique factor nested in plant; i suggest you try a model with slique factor and one without and choose the one that best fits the data.

Alberto Ferrari Thank you so much for your input, I'll start reading into these models and distributions. I think the answer I'm looking for is there, and includes a binomial approach and a nested model, as you were saying.

Thanks again!

A simple way to test the differences in counts is to test for the association with Genotype, with the chi-square test, collapsing over Siliques. One can collapse over Siliques in the 'prueba' example because there was no Silique effect nor an interaction with Genotype with the parametric ANOVA. See also Figure Rplot.png abova. The chi-square is very significant:

Fisher's Exact Test for Count Data ------------------------------------------------------------ Alternative hypothesis: two.sided p = 3.932838e-15

R script in attachment.

It is better to put the dependent variable (counts) in the rows and the independent variable (genotype) in the columns. See the adjusted script. A multiple comparison test is added. (mult1 is different form wt and from mult2)

I implemented the model suggested by @Alberto Ferrari. I am not familiar with such models. Yet, it seems to yield the same results and conclusions. I wrote and R script. See attachment. The input file (nico.csv) is also attached. Make sure to set the correct working directory in R in order to read it.

Jos Feys Thanks again! You're very kind. These days I will be checking in all you sent me. Once I finish, I will tell you how it went!

Thanks again!

Jos Feys, from the script it looks like you used a poisson model. Its applicability to binomial data, where counts have an upper bound, is debatable (in my experience it causes a HUGE inflation of type I errors). Also you treated Silique as a fixed effect with three levels and added an interaction term that seems unnecessary, since in my understanding Silique is a random effect nested in Plant.

If I am not mistaken the model should go like this:

model

Thanks Alberto. I must admit I am not familiar (to say the least) with such models. My knowledge stops with what I read in Andy Field's book 'Discovering Statistics Using R'. I indeed treated Silique as crossed with Plants, because each level of Silique co-occures with every Plant, as I understand. See:

https://www.theanalysisfactor.com/the-difference-between-crossed-and-nested-factors/

But, I agree, maybe Nicolás wanted Siliques to be nested. If so, he would need more plants for each level of Siliques.

In attachment the R script (and data) for the logistic regression suggested by Alberto Ferrari. I have not been able to make "Siliques" a random effect.

In attachment an analysis on the proportions aborted seeds with the fruits (siliques) as observational units (instead of the plants). In the former analyses I considered the siliques as 'fixed'. I assigned a sequence (1,2,3) to them. I considered these to be repeated measures over time or location. E.g. fruit 1 in the beginning of the harvest season, 2 in the middle, and 3 at the end of the harvest season. Or, e.g., 1 on top of the plant, 2 in the middel, and 3 at the bottom of the plant. I should have asked you before whether this was reasonable. It probably is much better to take the fruits themselves as ‘subjects’ or observational units. In this way they are considered to be random; as they should. (Alberto made me think ...)

Dear Jos Feys , Personally, I tried the line of code I pasted in my earlier response on your data and it worked.

Hello everyone! I'm starting my week with all the info that you so kindly sent me. I'll be looking into it in the next few hours, so later I'll tell you how it went!

I just read about the discussion about my data. I make observations in the siliques (fruits), as I count % of abortions in each one of them. The thing is, siliques I take from one plant are not independent (since they belong to the same plant) but they are completely independent from other siliques from other plants, even from the same genotype. That's the reason why, in my invented data set, I named each one of the siliques with a different name, since I never use one silique more than one time (I hope I'm making myself clear!). That's also the reason a "repeated - measure" model wouldn't apply in this case.

The model that Alberto Ferrari is suggesting seems correct to me, since I want to establish that siliques from each plant are nested (if I count in one silique, this silique is not used for any other plant, and all the siliques from one plant are not independent). I could treat "plant" as I random effect to check if it is important for total variance, although I believe it won't be. Silique is (to my understanding) definitely a random effect.

Jos Feys Alberto Ferrari Again, I can't thank you enough. I'll be trying all you've suggested me! I'll tell you later!

Nico

I edited this post to mention that the idea is to use more plants (at least 5 for each genotype, I could do more).

Hello again!

I've been reading all that you've done over the weekend. I enjoyed the discussion about using a binomial or a poisson aproach! From the courses I've taken, I believe binomial would be correct since the counts that I perform DO have an upper limit (% or proportion).

So, if I wanted to work with binomial, I would use a Generalized Linear Mixed Model (GLMM) in order to include "Genotype" as a fixed effect and both "Plant" and "Silique" as random effects.

From all the reading I've been doing, I believe I could create a few models that included (or didn't) the random effects that I listed earlier (in a similar way to what Jos Feys did in the 2nd "Nico" Rscript). I think it would be something like this:

#No random effect

model1

Yes, I think you got it right. The model has to "know" that all of your siliques are different. If siliques from different plants are coded in the same way, such as in Jos's example, then you need to use the "nested" notation so that the model understands that they are indeed different siliques, but if you assign from the start a unique code to each single silique then the crossed notation and the nested notation become equivalent.

AIC and BIC in GLMM context are a little bit tricky but they are applicable. You could have a look at predictions and residual plots too to be more confident in your final choice. Good job!

I have a script for the models you would like to test. I just added the siliques. Please set the working directory correctly at the top of the script.

I personally prefer the chisquare test.

Dear Alberto Ferrari , I believe you're right about AIC and BIC in GLMMs, since Jos Feys showed in the "nico_sil.R" script that a model considering "plant" and "silique" wouldn't be as a good fit as a model considering only "plant", and that seems odd to me, I don't know what you think. Also, it worries me that in the script, when Jos Feys performed the ANOVA, the table that shows the variance that random effects add to the system, it shows 0 for both random factors.

Dear Jos Feys , thanks again for taking the time. I like the approach you chose in both, I'll have to try with my real data and see (since differences may be not so clear as in this invented data set). Definitely, chisquare test is an insteresting option. To my understanding, it allows me to see if there are significant differences between something I have and something I expect.

There's another thing we do, where we compare % of abortions we obtain in a specific genotype to a % we would expect if the mutation in the genome would cause a gametophytic or a embryo defect (for example, 25%). Chisquare would definitely be a suitable option there..

Coming back to the system we've been working so far, I will start looking. It may be possible to combine chisquare test with mixed modeling. I'm pretty sure I should consider silique as a random factor somewhere, since we constantly see that it's a variable that adds considerable variance to the system.

If you have any other imput that may be useful please, be my guest. Thanks for all the help, I've exceeded all my expectations! Thanks!

I would report the results of the chisquare test, which is the most appropriate (non-parametric) test for count data. I added some comments in the chisquare R script, explaining expected values and row percentages.

Careful, the chi-square test on proportions "as is" is acceptable only if you can prove that there is no correlation structure in your outcomes, i.e., in our case, that seeds from the same plant and the same silique are not correlated. You see that this is the case when random effects estimates are very close to zero, which is what we see in the example data but probably not in real data. This aspect of the analysis requires great caution; whereas using a misspecified normal model (e.g. repeated measures ANOVA) on binomial data instead of a proper binomial model often gives very similar results, in my experience ignoring clustering structure can dramatically affect inference. Depending on the degree of within-plant correlation, choosing to ignore it may be inconsequential or make your entire analysis wrong.

For what concerns model selection, I'm not very fond of model selection methods such as AIC or BIC. In my opinion, it is more correct to build a model based on knowledge of the data generating process and then proceed to hypothesis testing. In your case, it just seems sensible to include at least a plant effect in the model, unless you have very sound biological reasons to assume that seeds from the same plant/silique have zero correlation (my biological background makes me believe this is unlikely).

I suggest that you have a look at random effect estimates from the model estimated on your data and drop a random effect only if the model does not converge or the corresponding random effect estimate is very very close to zero.

If you don't feel like using multilevel modeling, I'd suggest you just go for ANOVA on collapsed data; just do not ignore within-plant correlation unless you are very confident that it is zero.

Yes, the 3 observations of fruits per plant might not be independent from each other. They are nested within plants. I think one can circumvent this by taking only one fruit per plant as observational unit. The plants themselves were chosen randomly, I guess. I ran an ANOVA and a Chi-square analysis, respectively on the proportions aborted seeds and on the normal vs. aborted seed count, but each time only for the (in the data sheet) first observed fruit per plant.

In attachments the analyses, input data and results for the Logistic regression using GLM (generalized linear models).

This is definitely the most 'modern' approach.

Thanks Alberto for suggesting it.

In fact they all yield much more significant results for the effect of Genotype on the binary responses (normal vs. aborted). The attachments include the SAS way of analyzing logistic regression with proc GENMOD and CATMOD. (You can run SAS online for free on the "SAS on demands for academics" site.)

Hi again!

Dear Alberto Ferrari , thanks for all the advice. I agree that considering siliques and plants as random effects in a binomial models should be the solution, since seeds from the same silique and also siliques from the same plant are in fact correlated. Also, from now on, I'll be very careful with AIC and BIC, thanks again for the advice.

Dear Jos Feys , the idea of taking only one measurement from each plants would eliminate the problem of having to deal with random effects, but I believe it would force me to increase deeply the number of true replicates (plants) and this would be complicated. What do you think?

I've been reading into logistic regression but I think (from my understanding) that this model wouldn't be the best choice. Despite the fact that It deals with binomial data, I'm more interested in comparing levels from my independent variable (genotype) and that sounds more like an ANOVA to me.

I've asked some collegues in my university (that work in a different field, but have knowledge on statistics like you guys) and they have recommended me that I treat plants as replicates, build a model that considers siliques and plants as random effects (GLMM) and enter my continuous variable as a binomial column, where each row shows number of succeses (could be normal seeds) and failures (could be aborted seeds) per silique (fruit). They also said that I could perform an ANOVA "type II Wald" test, that results in a chi-square test.

What do you think about these suggestions? As usual, thank you very much for your kind replies!

Nico

I could be wrong, but I thought "type II Wald" tests are not relevant with a one-way ANOVA (with only one Between factor, i.c. Genotype). Correct me if I'm wrong.

Logistic regression is perfectly fine in this case, as long as it has random effects. Indeed, a binomial GLMM is just a logistic regression with random effects.

But if you just use logistic regression without random effects it's pretty much the same as using a simple chi-square test, and has the same problems with your data.

As I underdstand, if the counts of the seeds are dependent within plants if one takes 3 fruits per plant, one can consider these as repeated measures. I already ran such analyses. They yield a sign. effect of Genotype, including for the Wald type statistic. Mut1 is very sign. different from both WT and Mut2 in this example. See attachments. The SAS catmod procedures only runs when one does not include WT, because this genotype 'condition' does not have enough aborted seed counts.

My mistake! One cannot run an ANOVA on binary data as in the example. Rosenthal & Rosnow (1991, p.521) have argued one can run such an analysis, provided the study is not too small (say, df < 20), and the splits are not too extreme (say 90/10 percent). In the example there are only 2 observation (n = 2) and the proportions often are quite close to 0. (The frequency chi-square analysis is OK, preferable by siliques. The logistic regression, as well as the SAS proc catmod analysis aslo are OK.)

Ref.:

Rosenthal, R. & Rosnow, R.L. (1991). Essentials of Behavioral Research. Methods and Data Analysis. 2nd Ed. , McGraw-Hill, NY.

In attachment an R script for expanding frequencies, followed by the logistic regression.

You can test the differences in binary prortions between Mut1 and Mut2 with nonparametric tests such as the Mann-Whitney and or the Kruskal-Wallis test.

Dear Alberto Ferrari , based on your suggestions, I'll look deeply into Logistic Regression, as a GLMM that may be my best option. Thank you so much for the information. I should try with my original data set, since I have a concern about the great amount of 0s in it. If I tried to analyze my data as percentages (that was the original idea) the binomial distribution may not be the best fit. From what I've been reading (and also what you've told me earlier) beta distributuon may be appropiate.

Dear Jos Feys , as usual, thanks for all the information you're providing me. I'll start by saying that always remember my "true" data set has more plants (~5) and more siliques per plant (~5), so maybe one of the analyses you were suggesting could be performed anyway!

The part where you said "As I understand, if the counts of the seeds are dependent within plants if one takes 3 fruits per plant, one can consider these as repeated measures" left me thinking. It's a very interesting approach, one I had never considered. As always, I'll look into it.

I'll check all the scripts you've sent me. As always, thank you so much,

Nico.

Dear Nicolás Setzes , it is actually the opposite. The binomial distribution accomodates zeroes naturally, whereas the beta distribution is bound to values in the open interval (0,1), i.e. 0 and 1 themselves are excluded. You probably confused the beta distribution and the beta-binomial distribution (that's a very common misunderstanding). the beta-binomial is ok for your data and I mentioned it earlier, the beta is not since it does not tolerate zeroes.

However, beta binomial regression and mixed effect logistic regression, just as any other model, require a reasonable sample size (i.e. number of events) to estimate parameters. You won't be able to fit any model if your sample is entirely or almost entirely made up of zeroes (or ones), because there is no outcome for you to model there. But if you have THAT many zeroes, then every model that exists will fail; not even a simple chi-square test will do the trick, because it needs at least 5 observations per cells to be applied. In that case you will have to resort to a simple Fisher's exact test.

Bu that's a very extreme situation. I'm assuming you have enough events to apply logistic regression.

Just try to fit a GLMM; if you have non convergence then you can try something else.

To test whther the siliques are independent from each other I calculated correlation on the proportions of aborted seed between siliques. None of the Kendall tau values are sign., but with only 6 plants ... See R script 'correlation.r'.

I ran the nonparametric test as with the nparLD R package with SAS proc MIXED. See the SAS program 'nparMIXED.rtf' and its results 'nparMIXED-results.pdf'. This methode is described and illustrated in the Shah-Brunner and Shah&Madden attachments. There is also a SAS macro for it: F1_LD_F1.sas.

Dear Alberto Ferrari , thank you so much for the information. As you mentioned, I was quite confused between a beta distribution and a beta - binomial distribution. Now I understand beta alone is not an option, since I've read the same you're telling me now (that it doesn't include 0 and 1). I'm confident I have enough amount of samples in my "real" dataset (in the case I don't, I could get more). Thanks for your kind suggestions.

Dear Jos Feys , thank you so much for your latest scripts. Now I've learned how to analize possible correlation between variables, I didn't know how to do it in R. From what I can understand, the correlation tests do not agree that much on this dataset, I believe it may have to do with it having only a few replicates. I'll try it on my real dataset. I think the tests will indicate siliques from the same plant are indeed correlated, but that's only my opinion based on my biological background.

For differences between Pearson, Spearman and Kendall, see:

http://d-scholarship.pitt.edu/8056/1/Chokns_etd2010.pdf

or:

https://stats.stackexchange.com/questions/8071/how-to-choose-between-pearson-and-spearman-correlation

In the R script I added scatterplots. But, indeed with only 6 plants ...

Maybe also:

https://stats.stackexchange.com/questions/8071/how-to-choose-between-pearson-and-spearman-correlation

You can test the differences in binary proportion between Mut1 and Mut2 with non parametric tests such as the Mann-Whitney and or the Kruskal-Wallis test.

Haider Raid Talib Thanks for your response, it's one of the things I did!