1
$\begingroup$

I have an experiment that is designed as 6 blocks of 4 plots each, with two treatments (W_add and P_add) plus combination of treatments and control. The data are flux measurements taken during 9 campaigns.

I want to see if the treatments had a significant effect on the fluxes over time.

Here is my code:

library(nlme)
library(car)

setwd("[...]")

data<-read.csv("season_flux.csv")

flux = data$CO2_NEE #CO2_ER, CO2_NEE, GEP, or CH4

model = lme(flux ~ W_add * P_add * DOY, random = ~1|block/plot, 
                     corr = corGaus (form = ~DOY|block/plot, nugget = TRUE), 
                     data = data, na.action = na.omit, method = "REML")

contrasts(data$W_add) <- contr.sum
contrasts(data$P_add) <- contr.sum

Anova <- Anova(model, type = "III", test.statistic = "F")
summary(Anova)
print(Anova)

And here is the head of my data

  DOE       DOY   plot type CO2_ER  CO2_NEE        GEP        CH4   W_add  P_add
1 0.4831486 179   25    D 4.615219       NA -0.1560684 -2.4904458    NW    NP
2 0.4885537 179   25    L       NA 4.459151         NA         NA    NW    NP
3 0.4948533 179   26    D 4.178014       NA -1.3169272 -0.9756073     W    NP
4 0.5003533 179   26    L       NA 2.861087         NA         NA     W    NP
5 0.5069753 179   27    D 3.524708       NA -0.9260853 -1.2356268    NW     P
6 0.5123025 179   27    L       NA 2.598623         NA         NA    NW     P
  treatment block   PAR soil_moist soil_T_2 soil_T_5
1         C     1    NA   21.80000      4.6      3.0
2         C     1 206.2   21.80000      4.6      3.0
3         W     1    NA   32.13333      4.2      3.2
4         W     1 278.8   32.13333      4.2      3.2
5         P     1    NA   26.93333      4.2      2.6
6         P     1 303.4   26.93333      4.2      2.6

And the summary of my data:

      DOE               DOY             plot       type        CO2_ER      
 Min.   : 0.4832   Min.   :179.0   Min.   :25.00   D:214   Min.   :-1.044  
 1st Qu.:20.4576   1st Qu.:199.0   1st Qu.:31.00   L:182   1st Qu.: 2.320  
 Median :34.5007   Median :212.0   Median :37.00           Median : 4.050  
 Mean   :36.5580   Mean   :214.2   Mean   :36.57           Mean   : 4.376  
 3rd Qu.:54.6842   3rd Qu.:233.0   3rd Qu.:43.00           3rd Qu.: 5.268  
 Max.   :68.7054   Max.   :246.0   Max.   :48.00           Max.   :23.932  
                                                           NA's   :186     
    CO2_NEE            GEP                CH4           W_add    P_add   
 Min.   :-2.796   Min.   :-13.2704   Min.   :-8.03760   NW:199   NP:195  
 1st Qu.: 1.187   1st Qu.: -2.9642   1st Qu.:-2.56166   W :197   P :201  
 Median : 2.273   Median : -1.3474   Median :-1.85111                    
 Mean   : 2.541   Mean   : -2.1185   Mean   :-2.07423                    
 3rd Qu.: 3.587   3rd Qu.: -0.1575   3rd Qu.:-1.15698                    
 Max.   :10.662   Max.   :  0.0000   Max.   : 0.02807                    
 NA's   :214      NA's   :217        NA's   :189                         
 treatment     block            PAR           soil_moist       soil_T_2     
 C :101    Min.   :1.000   Min.   : 149.8   Min.   : 7.40   Min.   : 4.200  
 P : 98    1st Qu.:2.000   1st Qu.: 429.4   1st Qu.:13.98   1st Qu.: 7.000  
 W : 94    Median :4.000   Median : 659.5   Median :16.30   Median : 8.600  
 WP:103    Mean   :3.515   Mean   : 713.0   Mean   :16.38   Mean   : 8.624  
           3rd Qu.:5.000   3rd Qu.:1012.4   3rd Qu.:18.37   3rd Qu.: 9.900  
           Max.   :6.000   Max.   :1594.2   Max.   :32.13   Max.   :17.700  
                           NA's   :24                                       
    soil_T_5     
 Min.   : 2.600  
 1st Qu.: 5.900  
 Median : 6.800  
 Mean   : 6.809  
 3rd Qu.: 7.700  
 Max.   :11.000  
 NA's   :2     

I get the following output:

Error in lme.formula(flux ~ W_add * P_add * DOY, random = ~1 | block/plot,  :
  nlminb problem, convergence error code = 1
  message = singular convergence (7)

The error does not happen with another flux than CO2_NE or if I use another correlation. This error also does not happen if I delete the last line of my data, which looks like this:

68.70536503,246,48,L,,2.27190302873568,,,W,NP,W,6,450.2,17.0333333333333,9.2,7.1

Why is this happening?

Thanks a lot!

$\endgroup$
2

1 Answer 1

1
$\begingroup$

I have fitted the model using lme4 and found that there is no variation at the block level:

> m0 <- lmer(flux ~ W_add * P_add * DOY  + (1 | block) + (1 | plot), data = data)
> summary(m0)
Linear mixed model fit by REML ['lmerMod']
Formula: flux ~ W_add * P_add * DOY + (1 | block) + (1 | plot)
   Data: data

REML criterion at convergence: 762.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.85542 -0.64230  0.01562  0.52262  2.95350 

Random effects:
 Groups   Name        Variance Std.Dev.
 plot     (Intercept) 1.817    1.348   
 block    (Intercept) 0.000    0.000   
 Residual             2.781    1.668   
Number of obs: 182, groups:  plot, 24; block, 6

Fixed effects:
                   Estimate Std. Error t value
(Intercept)        2.034342   2.491762   0.816
W_addW             3.257081   3.731735   0.873
P_addP             2.614168   3.620830   0.722
DOY                0.002315   0.011295   0.205
W_addW:P_addP     -4.798173   5.210949  -0.921
W_addW:DOY        -0.017029   0.016990  -1.002
P_addP:DOY        -0.012404   0.016438  -0.755
W_addW:P_addP:DOY  0.025649   0.023681   1.083

lme4 fits the model without convergence problems, though it is unable to handle the correlation structure that you used with nlme. However, if we fit the model in nlme without block in the random structure we obtain very similar results to lme4:

> model = lme(flux ~ W_add * P_add * DOY, random = ~1|plot, corr = corGaus (form = ~DOY|plot, nugget = TRUE), data = data, na.action = na.omit, method = "REML")
> summary(model)
Linear mixed-effects model fit by REML
 Data: data 
      AIC      BIC    logLik
  785.879 823.7877 -380.9395

Random effects:
 Formula: ~1 | plot
        (Intercept) Residual
StdDev:    1.345257 1.671597

Correlation Structure: Gaussian spatial correlation
 Formula: ~DOY | plot 
 Parameter estimate(s):
       range       nugget 
3.176300e+00 3.482499e-06 
Fixed effects: flux ~ W_add * P_add * DOY 
                      Value Std.Error  DF    t-value p-value
(Intercept)        1.954598  2.516343 154  0.7767615  0.4385
W_addW             3.319879  3.771597  20  0.8802317  0.3892
P_addP             2.596632  3.656539  20  0.7101338  0.4858
DOY                0.002818  0.011427 154  0.2466313  0.8055
W_addW:P_addP     -4.852672  5.262861  20 -0.9220597  0.3675
W_addW:DOY        -0.017406  0.017206 154 -1.0116291  0.3133
P_addP:DOY        -0.012394  0.016626 154 -0.7454639  0.4571
W_addW:P_addP:DOY  0.025926  0.023957 154  1.0822138  0.2808

The random effects are almost identical and the fixed effects are very similar.

This leads me to believe that lme fails to converge because the optimizer gets stuck due to an estimated variance component being close to zero. I do not know enough about nlme to diagnose this any further.

$\endgroup$
3
  • $\begingroup$ Thank you very much! If I understand correctly, you are saying that there is no point in nesting the data in blocks? I was expecting some spatial variation as the field site was not completely homogeneous. $\endgroup$
    – Lours
    Mar 15, 2019 at 18:13
  • $\begingroup$ That's what my analysis suggests, yes. The variation at the plot level is almost as large as the residual variation, but lme4 estimated zero variation at the block level. This might be due to missing data, so one further approach is to use multiple imputation $\endgroup$ Mar 16, 2019 at 18:52
  • $\begingroup$ Thank you for the explanation! I have a better idea of what I'm doing now ;-) $\endgroup$
    – Lours
    Mar 17, 2019 at 18:13

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.