2
$\begingroup$

I have ecological data (transects as the unit) from inside marine reserves and matched control transects that I would like to test the difference across 18 separate response variables. These are 3 overall response variables: biomass, density, and richness, and 3 response variables (biomass, density and total Length) for each of 5 of individual fish families. However the data points are also grouped within sites (4 transects per site) and marine reserve (8 separate reserves). I want to test the overall difference for these 18 response variables, as well as the difference for each individual marine reserve. To do this I think I need to run mixed effects models with site nested in reserve as random effects for the overall outcomes, and mixed effects models with site only for each individual reserve outcome.

I have tried a number of things listed below but am running into several problems.

  1. For the overall outcomes I tried running these in nlme. The data is quite skewed for some of the 18 outcomes, with many 0's, so for these I have tried log transforming. However I get error # 1 when trying to log transform (see below)

  2. From there I tried running the models without log transforming them. Many of the residual plots were trumpet shaped as a result but I wanted to see if I could move forwards. The models ran but when I tried to put the main output data into a table format I get error # 2 (see below) for specific outcomes.

  3. If I rerun the models with 1|site instead of 1|name/site I can compile the data to a table, but the p-values of most of the outcome variables with both 1|name/site and 1|site are 1.

My first guess was that many of the outcome variables have a high number of 0's, which is causing problems.

Would it be possible to provide feedback on: A) The best models to run for this data B) How to compile this data from all 18 outcome variables into a single table C) An explanation of why the p - value is 1, and why this is leading to errors when trying to create the output table.

I am guessing this has to do with some groups of the random effects having no data, but I am unsure.

Thank you

```# Input data
matched.cases3 = read.csv('matched.cases2FHR.csv')
matched.cases3=matched.cases3%>% mutate(TorC=as.factor(TorC))
head(matched.cases3)

# Run all for one example from 18, Scarbm. TorC is Treatment or Control (0 or 1)

#try to log transform. See error #1
Scarbm.lme1 = lme(log(Scarbm)~TorC, 
random=~1|name/site,data=matched.cases3, method='REML')
#Log transform error: Error in MEEM(object, conLin, control$niterEM) : 
NA/NaN/Inf in foreign function call (arg 1)

#Run without log transforming
Scarbm.lme1 = lme(Scarbm~TorC, random=~1|name/site,data=matched.cases3, 
method='REML')
Scarbm.lme2 = lme(Scarbm~TorC, random=~TorC|name/site,data=matched.cases3, 
method='REML')

#check AIC fixed intercept is better
AIC(Scarbm.lme1,Scarbm.lme2)


#check residual plot, trumpet shaped
plot(Scarbm.lme1)


#Summary 
summary(Scarbm.lme1)
Linear mixed-effects model fit by maximum likelihood
Data: matched.cases3 

Random effects:
Formula: ~1 | name
        (Intercept)
StdDev:  0.05367244

  Formula: ~1 | site %in% name
         (Intercept) Residual
StdDev:    288.8733 334.5428

Fixed effects: Scarbm ~ TorC 
 Correlation: 
      (Intr)
TorC1 -0.357

Standardized Within-Group Residuals:
   Min         Q1        Med         Q3        Max 
-2.2320898 -0.3453362 -0.1419344  0.2529601  5.9144762 

Number of Observations: 250
Number of Groups: 
          name site %in% name 
             8             32 

AIC  BIC  logLik
3686.832    3704.439    -1838.416   



 Value Std.Error DF t-value p-value
(Intercept) 358.2354    59.57017    217 6.013671    0
TorC1           0.0000          42.48697    217 0.000000    1


#Compile to table (for all 18 outcome variables). See error #2 for when 
run as 1|name/site
lme=Scarbm.lme1
p=data.frame(summary(lme)$tTable[,"p-value"])
t=data.frame(summary(lme)$tTable[,"t-value"])
df=data.frame(summary(lme)$tTable[,"DF"])
ste=data.frame(summary(lme)$tTable[,"Std.Error"])
p=p[2,]
t=t[2,]
df=df[2,]
ste=ste[2,]
mean=emmeans(lme, ~TorC) %>% as.data.frame
mean=mean[c("emmean", "lower.CL","upper.CL")]
meanT=mean[2,]
meanC=mean[1,]
r2=data.frame(r.squaredGLMM(lme))
int=intervals(lme)
int=data.frame(int$fixed)
int=int[2,]
response=c("Scarbm")
response=data.frame(response)
management=c("FHR")
management=data.frame(management)
ScarbmTC=cbind(management,response,int,meanT,meanC,ste,df,t,p,r2)
#Error in intervals.lme(lme) : cannot get confidence intervals on var-cov 
components: Non-positive definite approximate variance-covariance Consider 
'which = "fixed"'

#rerun with random effect as 1|site only, it works. 
Overall=rbind(bmTC,densTC,richTC,AcanthbmTC,AcanthdensTC,AcanthLETC,ScarbmTC,ScardensTC,ScarLETC,LethbmTC,LethdensTC,LethLETC,LutbmTC,LutdensTC,LutLETC,SerrbmTC,SerrdensTC,SerrLETC)
Overall
#Example for just Scarbm row of Overall
management response lower est. upper emmean lower.CL ste df t p
FHR       Scarbm    -83.5972502 -9.077611e-14   83.5972502   
358.182621  235.455253 480.909989   42.4145630  217 -2.140211e-15    
1.000000e+00    
  1. Log transform error: Error in MEEM(object, conLin, control$niterEM) : NA/NaN/Inf in foreign function call (arg 1)

  2. intervals.lme(lme) : cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance Consider 'which = "fixed"'

$\endgroup$
2
  • 1
    $\begingroup$ Do you expect your data to result from two processes (first process: presence/absence, second process: amount)? If yes, you need an approach for zero-inflated data (which basically means two coupled models for the two processes). There is also a high chance that some of your dependent variables require use of a GLMM. And finally, this looks like a data exploration exercise. Your approach has unfortunate implications for resulting p-values. $\endgroup$
    – Roland
    Commented Aug 19, 2019 at 14:40
  • 1
    $\begingroup$ A lot of people use log(y+1) which doesn’t result in all those missing values. It isn’t necessarily theoretically justified, but on the other hand, if it results in more symmetrically distributed residuals, it would be a more desirable model than the one without a transformation. All models are wrong, but some are more useful than others. $\endgroup$
    – Russ Lenth
    Commented Aug 19, 2019 at 18:47

1 Answer 1

1
$\begingroup$

You could use a two-part mixed effects model for semi-continuous data. This combines a logistic regression for the dichotomous indicator that the outcome is zero or not with a log-normal model for the continuous part.

Two-part mixed effects models are available in the GLMMadaptive package I’ve written - you can find a sample analysis here.

$\endgroup$

Your Answer

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