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

[Marina Pastor]

Thank you Thierry,
Here is the reproducible example (I copied the wrong jelly column),  
but we obtained the same error:
- 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)?
Best wishes,
Marina Pastor

library("nlme")
library("MASS")
my.df <- data.frame (Jelly = c(11, 106, 0, 8, 4, 12, 0, 699, 42, 2),  
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)


Thierry Onkelinx <thierry.onkelinx using inbo.be> escribió:

> Dear Marina,
>
> The Poisson and negative binomial distributions assume counts =
> non-negative integers. Your response variable is not integer. You'll need
> the actual counts instead.
>
> Best regards,
>
> ir. Thierry Onkelinx
> Statisticus / Statistician
>
> Vlaamse Overheid / Government of Flanders
> INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND
> FOREST
> Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
> thierry.onkelinx using inbo.be
> Havenlaan 88 bus 73, 1000 Brussel
> www.inbo.be
>
> ///////////////////////////////////////////////////////////////////////////////////////////
> To call in the statistician after the experiment is done may be no more
> than asking him to perform a post-mortem examination: he may be able to say
> what the experiment died of. ~ Sir Ronald Aylmer Fisher
> The plural of anecdote is not data. ~ Roger Brinner
> The combination of some data and an aching desire for an answer does not
> ensure that a reasonable answer can be extracted from a given body of data.
> ~ John Tukey
> ///////////////////////////////////////////////////////////////////////////////////////////
>
> <https://www.inbo.be>
>
>
> Op di 17 nov. 2020 om 11:20 schreef Marina Pastor <marinapastor using icm.csic.es
>> :
>
>>
>> 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
>>
>> _______________________________________________
>> R-sig-mixed-models using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>


-- 
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