





Dear all,


I’m sitting
here with a huge data set from a monitoring program. After that I proved
several ways I have ended in hurdle models.


The
variable of inters at study is the amount of malformed embryos (malformed
embryo/fecund) in an amphipod living in the Baltic Sea. I have 3 categorical
explanatory variables (YEAR, SPECIES and STATIONS) with different numbers of
levels. The stations
are distributed along the coast, covering different bottoms, at different depth
and different length from the coast. The target of the study is to show the
status of reference stations and find if there are any differences between
stations, year and species (two amphipods).


The data is
right skewed, over-disperse with many (true) zeros and the number of amphipods
per station are not the same (unbalance). The embryos are dependent of each
female and each female has from a few to more than 50 embryos, that’s the
reason why I use weights.


Here is the
model: 


Library(pcsl)


hurdle(formula = Tot.dameggs ~ fSPECIES * fSTATN *
fYEAR, data = dat,   weights = FECUND,
dist = c("negbin"), zero.dist = c("binomial"))


Pearson residuals:

    Min      1Q 
Median      3Q     Max 

-6.2241
-3.1106 -2.2523  0.8126 61.0310 


Count model coefficients (truncated negbin with log
link):


                               Estimate Std.
Error z value Pr(>|z|)    

(Intercept)                     0.84196    0.07309 
11.520  < 2e-16 ***

fSPECIESP                      -0.18040    0.08426 
-2.141 0.032270 *  

fSTATN6020                      0.02653   
0.09140   0.290 0.771645    

fSTATN6022                     -0.14544    0.09863 
-1.475 0.140325    

fSTATN6025                      0.11754    0.09841  
1.194 0.232345    

fYEAR2013                      -0.76104    0.10263 
-7.415 1.21e-13 ***

fSPECIESP:fSTATN6020            0.16119    0.11293  
1.427 0.153485    

fSPECIESP:fSTATN6022            0.49364    0.12025  
4.105 4.04e-05 ***

fSPECIESP:fSTATN6025            0.42522    0.11910  
3.570 0.000357 ***

fSPECIESP:fYEAR2013            -0.14249    0.13694 
-1.041 0.298091    

fSTATN6020:fYEAR2013            0.52760    0.15009  
3.515 0.000440 ***

fSTATN6022:fYEAR2013            0.78145    0.13136  
5.949 2.70e-09 ***

fSTATN6025:fYEAR2013           -0.47691    0.16669 
-2.861 0.004222 ** 

fSPECIESP:fSTATN6020:fYEAR2013 -0.50072    0.19299 
-2.595 0.009471 ** 

fSPECIESP:fSTATN6022:fYEAR2013  0.36693   
0.17071   2.149 0.031604 *  

fSPECIESP:fSTATN6025:fYEAR2013  0.46022   
0.20135   2.286 0.022276 *  


Log(theta)                     -0.16826    0.04955 
-3.396 0.000684 ***

 

Zero hurdle model coefficients (binomial with logit
link):


                               Estimate Std.
Error z value Pr(>|z|)    

(Intercept)                    -1.11365    0.06567 -16.958  < 2e-16 ***

fSPECIESP             
         1.70977    0.08862 
19.293  < 2e-16 ***

fSTATN6020                      0.45981    0.08699  
5.286 1.25e-07 ***

fSTATN6022                      0.68621    0.09540  
7.193 6.35e-13 ***

fSTATN6025                     -0.23911    0.09041 
-2.645  0.00817 ** 

fYEAR2013                      -0.08152    0.08878 
-0.918  0.35852    

fSPECIESP:fSTATN6020           -0.87573    0.12014 
-7.289 3.12e-13 ***

fSPECIESP:fSTATN6022           -1.13412    0.12848 
-8.827  < 2e-16 ***

fSPECIESP:fSTATN6025           -0.14570    0.12456 
-1.170  0.24212    

fSPECIESP:fYEAR2013            -0.91383    0.12601 
-7.252 4.10e-13 ***

fSTATN6020:fYEAR2013           -0.23149    0.13579 
-1.705  0.08825 .  

fSTATN6022:fYEAR2013            0.80387    0.12509  
6.426 1.31e-10 ***

fSTATN6025:fYEAR2013                  -0.21088   
0.13462  -1.566  0.11724   


fSPECIESP:fSTATN6020:fYEAR2013  1.17060   
0.18432   6.351 2.14e-10 ***

fSPECIESP:fSTATN6022:fYEAR2013  1.75151   
0.17903   9.783  < 2e-16 ***

fSPECIESP:fSTATN6025:fYEAR2013  1.42613   
0.18108   7.876 3.39e-15 ***


---Signif. codes: 
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 


Theta: count = 0.8451


Number of iterations in BFGS optimization: 25 


Log-likelihood: -2.912e+04 on 33 Df


I read some
chapters in “Mixed effects models and experimental in ecology with R” but I still
have some questions.


1)      I’m not sure that weights are
properly used for the binomial part.

2)      It is possible to compare the
stations and interaction between them (post hoc test), as I have understand the     model compare always the different variables with the intercept that are a baseline
(e.g. STATN1,YEAR1 and SPECIES1). 

3)      Is it possible to calculate the power
of a hurdle model by use the interclass correlation coefficient?

 

 Thanks in advance

/Matias

  		 	   		   		 	   		  
	[[alternative HTML version deleted]]

