[R-sig-ME] Trying to modify variance structure

[Marina Pastor]

Hello,
We are studying the jellyfish distribution and we would like to know  
if it depends of temperature and salinity. As our database are counts  
and we have a high amount of 0, but also other high numbers around  
400, we used a Negative Binomial Generalised Lineal Model. As the  
filtered volume was not the same for all the counts, we incorporated  
this through an offset.
Checking the residual distribution we realised we were violating the  
homogeneity assumption since the larger the temperature, the larger  
the variation.
We first tried to apply different variance structures through  
varFixed, varPower and varExp functions (nlme package) depending on  
the temperature, but we did not manage.
The reproducible example is below.
- For varFixed(~Temp) we followed the R help example. Why we obtained  
this error?: Error in `$<-.data.frame`(`*tmp*`, VarFixedT, value =  
numeric(0)) :   replacement has 0 rows, data has 10
- For varPower(1, form =~Temp) and varExp(1, form =~ Temp), how to  
choose the best “value” before “form” in the formula (“1” in the  
example)?
We would be very grateful if you could help us, we don’t know how to  
improve our model and is one of the requirements to publish our  
article. Many thanks for your time in advance.
Best wishes,
Marina Pastor


Reproducible example:
library("nlme")
library("MASS")

my.df <- data.frame (Jelly = c(1.13, 10.98, 0.00, 0.97, 0.62, 1.04,  
0.00, 77.83, 4.12, 0.18), Temp = c(24.63, 24.61, 24.63, 25.64, 25.63,  
26.22, 26.17, 25.34, 25.44, 25.09), Sal = c(37.16, 36.79, 38.06,  
38.20, 38.15, 38.26, 38.25, 38.10, 38.07, 37.96), Vol = c(971.0,  
965.5, 835.0, 823.0, 640.0, 1147.0, 1322.0, 912.0, 1018.0, 1095.0))


my.df$VarFixedT <- varFixed(~Temp)
GLMNB_W <- glm.nb(Jelly ~ Temp + Sal +
                     offset(log(Vol)),
                   data= my.df,
                   weights = VarFixedT)

my.df$VarPowerT <- varPower(1, form =~Temp)
GLMNB_W <- glm.nb(Jelly ~ Temp + Sal +
                     offset(log(Vol)),
                   data= my.df,
                   weights = VarPowerT)
summary(GLMNB_W)

my.df$VarExpT <- varExp(1, form =~ Temp)
GLMNB_W <- glm.nb(Jelly ~ Temp + Sal +
                     offset(log(Vol)),
                   data= my.df,
                   weights = VarExpT)
summary(GLMNB_W)

-- 
Marina Pastor
PhD Student
Marine biology and oceanography department
Institut de Ciències del Mar (ICM-CSIC)
Spanish National Research Council
Passeig Marítim de la Barceloneta 37-49, E-08003 Barcelona, Catalonia, Spain
Phone: +34 932309500 (ext. 1113)
E-mail: marinapastor using icm.csic.es