3
$\begingroup$

I wanted to analyze my data below to see the effect of treatments (CC,CCS, CS and SCS), soil health indicators (B.glucosidase, Protein and POX.C) and Locations (M and U) on yield. There are four blocks (1,2,3,4) in my data. For the analysis, I used lm with R as below:

library(finalfit)
library(dplyr)
dependent = "yield"
explanatory = c("B.glucosidase", "Protein", "POX.C", "Location", "treatment")

df %>%
  finalfit.lm(dependent, explanatory)

However, the reviewers are telling me that lm is not suitable for this data and are asking me to use Generalized Linear Mixed Models using PROC GLIMMIX in SAS. I am not clear about how GLIMMIX in SAS differs from lm or glm in R and would like someone here to share some knowledge explaining why I have to resort to GLIMMIX for this analysis. Also, I am not familiar to SAS and wanted to use GLIMMIX alternative in R (perhaps glmer?). Could someone please help me understand this clearly. Thank you for your help in advance!

data:

Location    treatment   block   B.glucosidase   Protein POX.C   yield
M   CC  1   1.600946    6212.631579 810.3024    5156
M   CC  2   1.474084    5641.403509 835.5242    5157
M   CC  3   1.433078    4392.280702 856.206 5551
M   CC  4   1.532492    7120.701754 759.8589    5156
M   CCS 1   1.198667    5305.964912 726.2298    4804
M   CCS 2   1.193193    4936.842105 792.6472    4720
M   CCS 3   1.214941    5383.157895 724.7165    4757
M   CCS 4   1.360981    6077.894737 699.3266    5021
M   CS  1   1.853056    5769.122807 500.9153    4826
M   CS  2   1.690117    5016.842105 634.8698    4807
M   CS  3   1.544357    5060.350877 637.9536    4475
M   CS  4   1.825132    5967.017544 648.8814    4596
M   SCS 1   1.695409    5576.842105 641.0357    4669
M   SCS 2   1.764123    5174.035088 623.3822    4588
M   SCS 3   1.903743    5655.438596 555.2834    4542
M   SCS 4   1.538684    5468.77193  520.8119    4592
U   CC  1   0.845077    7933.333333 683.3528    5583
U   CC  2   1.011463    7000    595.9173    5442
U   CC  3   0.857032    6352.982456 635.4315    5693
U   CC  4   0.989803    8153.684211 672.4234    5739
U   CCS 1   0.859022    6077.894737 847.2944    5045
U   CCS 2   0.919467    4939.649123 745.5665    4902
U   CCS 3   1.01717 5002.807018 778.3548    5006
U   CCS 4   0.861689    6489.122807 735.8141    5086
U   CS  1   0.972332    4694.035088 395.2647    4639
U   CS  2   0.952922    5901.052632 570.4148    4781
U   CS  3   0.804431    4303.859649 458.0383    4934
U   CS  4   0.742634    6768.421053 535.3851    4857
U   SCS 1   1.195837    6159.298246 678.0293    4537
U   SCS 2   1.267285    6090.526316 670.7419    4890
U   SCS 3   1.08571 4939.649123 335.2923    4842
U   SCS 4   1.20097 5262.45614  562.5674    4608

or read the data in R as:

df <- read.table(header= TRUE, text=
                   "Location    treatment   block   B.glucosidase   Protein POX.C   yield
M   CC  1   1.600946    6212.631579 810.3024    5156
M   CC  2   1.474084    5641.403509 835.5242    5157
M   CC  3   1.433078    4392.280702 856.206 5551
M   CC  4   1.532492    7120.701754 759.8589    5156
M   CCS 1   1.198667    5305.964912 726.2298    4804
M   CCS 2   1.193193    4936.842105 792.6472    4720
M   CCS 3   1.214941    5383.157895 724.7165    4757
M   CCS 4   1.360981    6077.894737 699.3266    5021
M   CS  1   1.853056    5769.122807 500.9153    4826
M   CS  2   1.690117    5016.842105 634.8698    4807
M   CS  3   1.544357    5060.350877 637.9536    4475
M   CS  4   1.825132    5967.017544 648.8814    4596
M   SCS 1   1.695409    5576.842105 641.0357    4669
M   SCS 2   1.764123    5174.035088 623.3822    4588
M   SCS 3   1.903743    5655.438596 555.2834    4542
M   SCS 4   1.538684    5468.77193  520.8119    4592
U   CC  1   0.845077    7933.333333 683.3528    5583
U   CC  2   1.011463    7000    595.9173    5442
U   CC  3   0.857032    6352.982456 635.4315    5693
U   CC  4   0.989803    8153.684211 672.4234    5739
U   CCS 1   0.859022    6077.894737 847.2944    5045
U   CCS 2   0.919467    4939.649123 745.5665    4902
U   CCS 3   1.01717 5002.807018 778.3548    5006
U   CCS 4   0.861689    6489.122807 735.8141    5086
U   CS  1   0.972332    4694.035088 395.2647    4639
U   CS  2   0.952922    5901.052632 570.4148    4781
U   CS  3   0.804431    4303.859649 458.0383    4934
U   CS  4   0.742634    6768.421053 535.3851    4857
U   SCS 1   1.195837    6159.298246 678.0293    4537
U   SCS 2   1.267285    6090.526316 670.7419    4890
U   SCS 3   1.08571 4939.649123 335.2923    4842
U   SCS 4   1.20097 5262.45614  562.5674    4608
")
$\endgroup$
  • 1
    $\begingroup$ Briefly, the reviewers are asking you to use a hierarchical model (aka random effects model). Assuming block is something like a plot of land (I'm not familiar with your data or your field, and you shouldn't assume everyone is either), then it looks like multiple treatments to the same 4 blocks. Linear or generalized linear models don't account for the fact that you're treating the same blocks, and that responses within each block might be correlated. Mixed models, which PROC GLIMMX fits, do. R has equivalent packages, e.g. lme4. $\endgroup$ – Weiwen Ng Oct 31 '19 at 19:03
5
$\begingroup$

The reviewers are asking for you to run a mixed effects model. This is because your data are clustered, due to repeated measurements in each block. Each block is likely to have it's own characteristics and therefore measurements in each block are likely to be more similar to each other than to measurements in other blocks. In other words, measurements in each block are likely to be correlated. This can be accounted for by fitting fixed effects or random effects for block, or with generalised estimating equations. The reviewers are asking for random effects, that is, a mixed-effects model. It is a little odd that they specify PROC GLIMMIX as this is usually used for generalised mixed models, but it can also fit models where the conditional distribution of the dependent variable is normal.

Since you only have 4 blocks it is questionable whether this is sufficient to be able to estimate a variance for the random effect of block from only 4 observations. I would suggest that treating block as a fixed effect and using lm() in R would be fine.

But if you must follow their advice, then you can use lmer() from the lme4 package. If you do need a generalised mixed model then you can use glmer() from the same package.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks for your answer. In that case, I was wondering if I could use something like this cmod_lme4C_L <- glmer(yield~ Location + treatment + B.glucosidase + Protein + POX.C + (1|block),data=df)? $\endgroup$ – MAPK Nov 4 '19 at 19:40
  • 1
    $\begingroup$ You would need to specify the family parameter if using glmer and that will depend on your data, but other than that, probably yes. $\endgroup$ – Robert Long Nov 4 '19 at 19:48
  • 1
    $\begingroup$ @MAPK, actually, since you are labeling the blocks 1, 2, 3, 4 in each location, you'd have to specify in the model that the blocks are nested within the location. That is, block 1 at u isn't the same block as block 1 at m. A viable approach may simply be to label the blocks 1 through 8. $\endgroup$ – Sal Mangiafico Nov 10 '19 at 13:01
  • $\begingroup$ @SalMangiafico That makes sense. Thanks! $\endgroup$ – MAPK Nov 10 '19 at 23:43
  • $\begingroup$ @SalMangiafico is quite correct. Another way to deal with it is to specify the random intercepts as (1|Location:block) $\endgroup$ – Robert Long Nov 12 '19 at 18:24
0
$\begingroup$

You can do this with a small adjustment to the code you already have:

library(finalfit)
library(dplyr)
dependent = "yield"
explanatory = c("B.glucosidase", "Protein", "POX.C", "Location", "treatment")

df %>%
  finalfit(dependent, explanatory, random_effect = "(1 | Location / block)")

This runs a nested random intercept model with lme4::lmer() as described above. The t-statistic in the output is converted to a p-value, but take care with this as it assumes a lot more data than you have.

You may want to bootstrap the confidence intervals and/or p-value of the model. Come back if you do as that is a different question.

| cite | improve this answer | |
$\endgroup$

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.