1
$\begingroup$

I have a completely randomized block design with 4 treatments and 3 replicates per treatment (12 parcels). Each replicate is placed within a block. In each parcel I have measured the number of infested leaves and the total amount of leaves at different times in the year. The total number of leaves is not equal among parcels. I want to test whether there are differences on the proportion of infested leaves among treatments considering that I have repeated measures in time of the same plants and a blocking factor.

Alternatively, I am interested in knowing whether the variation in the proportion (proportion of the last sampling minus proportion of the first sampling) Is different among treatments.

Which approach do you suggest me? Which test is more appropriate considering the two options?

Edit 19/02/2024: Following the received advices, I include the info on the dataset:

'data.frame':   2645 obs. of  18 variables:
 $ ID_pianta                : chr  "_Pianta_1S3C" "_Pianta_2S3C" "_Pianta_3S3C" "_Pianta_4S3C" ...
 $ Data_rilievo             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Coltura                  : chr  "Spinacio" "Spinacio" "Spinacio" "Spinacio" ...
 $ Nome_parcella            : chr  "S3C" "S3C" "S3C" "S3C" ...
 $ Blocco                   : num  3 3 3 3 3 3 3 3 3 3 ...
 $ Trattamento              : chr  "Controllo" "Controllo" "Controllo" "Controllo" ...
 $ Infestate                : num  0 1 0 0 0 1 0 0 1 0 ...
 $ Totale_foglie            : num  5 5 5 4 6 3 5 9 6 6 ...
 $ Numero_foglie_minate     : num  0 1 0 0 0 1 0 0 3 0 ...
 $ Percentuale_foglie_minate: num  0 0.2 0 0 0 ...
 $ N_mine_serpentina        : num  NA NA NA NA NA NA NA NA NA NA ...
 $ N_mine_macchia           : num  NA 1 NA NA NA 1 NA NA 3 NA ...
 $ Totale mine              : num  NA 1 NA NA NA 1 NA NA 3 NA ...
 $ Foglie_sane              : num  5 4 5 4 6 2 5 9 3 6 ...
 $ log_Numero_foglie_minate : num  0 0.693 0 0 0 ...
 $ log_Foglie_sane          : num  1.79 1.61 1.79 1.61 1.95 ...

and also some informative plots Percentage mined leaves among treatments:

library(flexplot)
flexplot(Percentuale_foglie_minate~Trattamento,spread="quartiles", jitter=c(.1, 0), data = data_infestate)

Percentage mined leaves

Distribution of mined leaves:

flexplot(Numero_foglie_minate~1, data = data_infestate)

Distribution of mined leaves

Then, I applied the following model to test whether the infested leaves "Numero_foglie_minate" and the healthy leaves "Foglie_sane" differ among treatments

MODEL 1

   summary(model_interazione)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: cbind(Numero_foglie_minate, Foglie_sane) ~ Trattamento + Blocco +      Data_rilievo + (1 | ID_pianta)
   Data: data_infestate
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))

     AIC      BIC   logLik deviance df.resid 
  4066.1   4142.5  -2020.0   4040.1     2632 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6876 -0.4813 -0.3308  0.1900  4.8004 

Random effects:
 Groups    Name        Variance Std.Dev.
 ID_pianta (Intercept) 1.307    1.143   
Number of obs: 2645, groups:  ID_pianta, 478

Fixed effects:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -2.98258    0.18485 -16.135  < 2e-16 ***
TrattamentoLavanda   -0.43253    0.19115  -2.263   0.0237 *  
TrattamentoRosmarino -0.40539    0.20356  -1.992   0.0464 *  
TrattamentoTimo      -0.22146    0.18686  -1.185   0.2360    
Blocco2              -0.33823    0.17387  -1.945   0.0517 .  
Blocco3              -0.04462    0.17187  -0.260   0.7951    
Data_rilievo9         0.03357    0.12291   0.273   0.7848    
Data_rilievo16       -0.01192    0.12264  -0.097   0.9226    
Data_rilievo22        0.08879    0.12123   0.732   0.4639    
Data_rilievo29       -0.03263    0.12380  -0.264   0.7921    
Data_rilievo36       -0.32984    0.12850  -2.567   0.0103 *  
Data_rilievo52       -0.56837    0.12778  -4.448 8.67e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) TrttmL TrttmR TrttmT Blocc2 Blocc3 Dt_rl9 Dt_r16 Dt_r22 Dt_r29 Dt_r36
TrttmntLvnd -0.468                                                                      
TrttmntRsmr -0.411  0.463                                                               
TrattamntTm -0.512  0.495  0.463                                                        
Blocco2     -0.441 -0.034 -0.127  0.018                                                 
Blocco3     -0.436 -0.057 -0.153 -0.002  0.546                                          
Data_riliv9 -0.324 -0.007  0.046 -0.007 -0.046 -0.044                                   
Data_rilv16 -0.331 -0.004  0.044 -0.003 -0.040 -0.039  0.539                            
Data_rilv22 -0.329 -0.009  0.037 -0.016 -0.043 -0.040  0.546  0.551                     
Data_rilv29 -0.326 -0.003  0.039 -0.015 -0.041 -0.041  0.536  0.542  0.553              
Data_rilv36 -0.317  0.004  0.047 -0.003 -0.041 -0.046  0.515  0.519  0.531  0.531       
Data_rilv52 -0.318  0.001  0.046 -0.010 -0.038 -0.036  0.520  0.525  0.536  0.540  0.518

Seems there are significant effects for treatments and also for the last sampling dates. I have used the package DHARMa for some diagnostics and this is what I get:

> library(DHARMa)
> # Simulate standardized residuals
> simulated_residuals <- simulateResiduals(model_interazione)
> # Plot diagnostic plots
> plot(simulated_residuals)

Residuals analysis

Edit 20/02/2024

Following suggestions, I add the results of the glmer with interaction between treatment and time. I had to increase the number of interactions to eliminate convergence issues:

MODEL 2

    > model_interazione <- glmer(cbind(Numero_foglie_minate, Foglie_sane) ~ Trattamento*Data_rilievo + Blocco + (1 | ID_pianta),
+                            data = data_infestate,family=binomial, 
+                            glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))
> summary(model_interazione)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: cbind(Numero_foglie_minate, Foglie_sane) ~ Trattamento * Data_rilievo +      Blocco + (1 | ID_pianta)
   Data: data_infestate
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))

     AIC      BIC   logLik deviance df.resid 
  4070.8   4253.1  -2004.4   4008.8     2614 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6316 -0.4806 -0.3263  0.1601  4.4388 

Random effects:
 Groups    Name        Variance Std.Dev.
 ID_pianta (Intercept) 1.313    1.146   
Number of obs: 2645, groups:  ID_pianta, 478

Fixed effects:
                                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)                         -2.946505   0.225251 -13.081  < 2e-16 ***
TrattamentoLavanda                  -0.367148   0.290348  -1.265  0.20605    
TrattamentoRosmarino                -0.229666   0.285925  -0.803  0.42184    
TrattamentoTimo                     -0.761728   0.306496  -2.485  0.01294 *  
Data_rilievo9                        0.109010   0.218369   0.499  0.61764    
Data_rilievo16                      -0.001370   0.218764  -0.006  0.99500    
Data_rilievo22                      -0.159518   0.230889  -0.691  0.48964    
Data_rilievo29                      -0.039270   0.223567  -0.176  0.86057    
Data_rilievo36                      -0.437104   0.227100  -1.925  0.05426 .  
Data_rilievo52                      -0.836964   0.234122  -3.575  0.00035 ***
Blocco2                             -0.293181   0.175388  -1.672  0.09460 .  
Blocco3                             -0.008417   0.173220  -0.049  0.96125    
TrattamentoLavanda:Data_rilievo9    -0.348587   0.333961  -1.044  0.29658    
TrattamentoRosmarino:Data_rilievo9  -0.135294   0.344371  -0.393  0.69441    
TrattamentoTimo:Data_rilievo9        0.217114   0.342126   0.635  0.52569    
TrattamentoLavanda:Data_rilievo16   -0.018477   0.326907  -0.057  0.95493    
TrattamentoRosmarino:Data_rilievo16 -0.365334   0.353116  -1.035  0.30086    
TrattamentoTimo:Data_rilievo16       0.307753   0.342737   0.898  0.36922    
TrattamentoLavanda:Data_rilievo22    0.021807   0.341696   0.064  0.94911    
TrattamentoRosmarino:Data_rilievo22 -0.220751   0.361474  -0.611  0.54140    
TrattamentoTimo:Data_rilievo22       1.006785   0.335352   3.002  0.00268 ** 
TrattamentoLavanda:Data_rilievo29   -0.069308   0.336851  -0.206  0.83698    
TrattamentoRosmarino:Data_rilievo29 -0.688504   0.378488  -1.819  0.06890 .  
TrattamentoTimo:Data_rilievo29       0.569437   0.337192   1.689  0.09127 .  
TrattamentoLavanda:Data_rilievo36    0.031248   0.348233   0.090  0.92850    
TrattamentoRosmarino:Data_rilievo36 -0.318984   0.379486  -0.841  0.40059    
TrattamentoTimo:Data_rilievo36       0.657741   0.349228   1.883  0.05964 .  
TrattamentoLavanda:Data_rilievo52   -0.043118   0.362352  -0.119  0.90528    
TrattamentoRosmarino:Data_rilievo52  0.176141   0.373717   0.471  0.63741    
TrattamentoTimo:Data_rilievo52       0.954141   0.345766   2.759  0.00579 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation matrix not shown by default, as p = 30 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

The results between MODEL 1 and MODEL 2 are quite different:

  • MODEL 1: Lavanda and Rosmarino treatments, and some times are significant
  • MODEL 2: Timo treatment and some Timo*Time are significant

Now, the question:

  • Why the models report such different results?
  • Which model shall I trust more? I am puzzled about which variables should be considered as fixed or random effects
  • How to further check the validity of the model in terms of predictive capacities and residual analysis?
$\endgroup$
25
  • 1
    $\begingroup$ It may be more productive to formulate a model for your experimental data. (Once you have a model that fits the data well, you should be able to estimate the difference in proportion of infested leaves between conditions.) One model to consider is a binomial GLM. (It can handle different number of leaves per data point but there is an assumption that the leaves in the same parcel are independent.) $\endgroup$
    – dipetkov
    Feb 18 at 10:44
  • 2
    $\begingroup$ Your data is aggregated: one row represents one experimental unit = plant, with two observations for each plant, number of infected leaves and number of leaves total. That format should work. Please take a look at the thread I linked above, esp. the part cbind(number_of_dead_aphids, number_of_live_aphids). Also, if you are okay with posting the data on CV, that tends to help with getting on-point advice. $\endgroup$
    – dipetkov
    Feb 18 at 20:35
  • 2
    $\begingroup$ You need to replace NA with (presumably) 0 for the infested leaves, because right now all of those observations get ignored, leaving you with 806. Then you should also run the binomial model as @dipetkov explained. The lme4-function for a glmm is glmer . $\endgroup$ Feb 19 at 16:27
  • 1
    $\begingroup$ Good progress! The question also mentions a blocking factor (which can be included in the model as a fixed effect) and a time factor. I'm not sure which column indicates the time. How many time points are there in the data? $\endgroup$
    – dipetkov
    Feb 19 at 21:33
  • 1
    $\begingroup$ I think that the model should have a fixed time effect. (Otherwise you want be able to compare the average treatment effects at the first sampling / the last sampling.) One alternative is to include time as fixed factor(Data_rilievo), which will take 6 degrees of freedom since there are 7 time points. (Another option is to include a spline of time. This would be more complex; the idea is to fit a flexible function for the time effect, not assume it's linear). $\endgroup$
    – dipetkov
    Feb 19 at 22:15

0

Your Answer

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

Browse other questions tagged or ask your own question.