R help - Feb 2020 - Help to solve modeling problem with gamm

I conducted an experiment where earthworms were subjected to two treatments,
with and without herbicide in the soil. Biomass measurements were taken every 12
days for 398 days and the biomass growth curves as a function of time were
plotted.

There was clearly a non-linear growth pattern such that an additive mixed
effects model was proposed to model the behavior of biomass as a function of
time and treatments.

When plotting the residuals a clear cone-shaped pattern was observed, therefore
a series of additive models were proposed sequentially to deal with violations
of the assumption of homogeneity. Below we can see the models with the following
names: M.1; M.2; M.3; M.4



lmc <- lmeControl (niterEM = 5000, msMaxIter = 1000)

f1 <- formula (Biomass ~ Treat + s (Time, by = Treat))



M.1 <-gamm (f1, random = list (fcajita = ~ 1), method = "REML",
control = lmc, data = Acorticis)



#This first model uses the experimental box factor (i.e. fcajita) as the random
element of the model. This random effects model assumes homogeneity between the
experimental boxes and within them over time



M.2 <-gamm (f1, random = list (fcajita = ~ 1), method = "REML",
control = lmc, data = Acorticis, weights = varIdent (form = ~ 1 | fcajita))



#This second model assumes heterogeneity between boxes, but homogeneity within
each box over time



M.3 <- gamm (f1, random = list (fcajita = ~ 1), method = "REML",
control = lmc, data = Acorticis, weights = varExp (form = ~ Time10))



#The third model assumes homogeneity between boxes but heterogeneity within each
box over time



Finally, we decided to model the heterogeneity using the 'varComb'
function in order to combine the variances where the model allows heterogeneity
between the experimental boxes and heterogeneity within the experimental boxes
over time:



M.4 <- gamm (f1, random = list (fcajita = ~ 1), data = Acorticis, method =
"REML", control = lmc, weights = varComb (varIdent (form = ~ 1 |
fcajita), varPower (form = ~ Time10)))



The first three models executed perfectly and the following values ??of the AIC
indicator were obtained:
> AIC (M.1 $ lme, M.2 $ lme, M.3 $ lme)


        df       AIC



M.1     8        379.6464



M.2    15        309.5736



M.3     9        310.4828



Unfortunately, the execution of the M.4 model failed and the following error
message was obtained:



Error in environment (attr (ret $ lme $ modelStruct $ varStruct,
"formula")) <-. GlobalEnv:

attempt to set an attribute on NULL



A final model I tried was M5:

M.5 <- gamm(f1, random = list(fcajita =~ 1), data = Acorticis, method =
"REML", control = lmc, weights = varComb(varIdent(form = ~1|fcajita),
varExp(form =~ Time10|fcajita)))

and this time I got the following error message:

Error in lme.formula(y ~ X - 1, random = rand, data = strip.offset(mf),  :

  nlminb problem, convergence error code = 1

  message = function evaluation limit reached without convergence (9)

Adem?s: Warning message:

In logLik.reStruct(object, conLin) :

  Singular precision matrix in level -1, block 1



My question is: Could someone help me fix these problems to run the M.4 and M.5
models?


	[[alternative HTML version deleted]]

Wrong list.

Post here:
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

in **plain text** not html.

Bert Gunter

"The trouble with having an open mind is that people keep coming along and
sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )


On Fri, Feb 28, 2020 at 12:25 AM Jos? Antonio Garc?a P?rez <
garci95 at hotmail.com> wrote:
> I conducted an experiment where earthworms were subjected to two
> treatments, with and without herbicide in the soil. Biomass measurements
> were taken every 12 days for 398 days and the biomass growth curves as a
> function of time were plotted.
>
> There was clearly a non-linear growth pattern such that an additive mixed
> effects model was proposed to model the behavior of biomass as a function
> of time and treatments.
>
> When plotting the residuals a clear cone-shaped pattern was observed,
> therefore a series of additive models were proposed sequentially to deal
> with violations of the assumption of homogeneity. Below we can see the
> models with the following names: M.1; M.2; M.3; M.4
>
>
>
> lmc <- lmeControl (niterEM = 5000, msMaxIter = 1000)
>
> f1 <- formula (Biomass ~ Treat + s (Time, by = Treat))
>
>
>
> M.1 <-gamm (f1, random = list (fcajita = ~ 1), method =
"REML", control > lmc, data = Acorticis)
>
>
>
> #This first model uses the experimental box factor (i.e. fcajita) as the
> random element of the model. This random effects model assumes homogeneity
> between the experimental boxes and within them over time
>
>
>
> M.2 <-gamm (f1, random = list (fcajita = ~ 1), method =
"REML", control > lmc, data = Acorticis, weights = varIdent (form =
~ 1 | fcajita))
>
>
>
> #This second model assumes heterogeneity between boxes, but homogeneity
> within each box over time
>
>
>
> M.3 <- gamm (f1, random = list (fcajita = ~ 1), method =
"REML", control > lmc, data = Acorticis, weights = varExp (form = ~
Time10))
>
>
>
> #The third model assumes homogeneity between boxes but heterogeneity
> within each box over time
>
>
>
> Finally, we decided to model the heterogeneity using the 'varComb'
> function in order to combine the variances where the model allows
> heterogeneity between the experimental boxes and heterogeneity within the
> experimental boxes over time:
>
>
>
> M.4 <- gamm (f1, random = list (fcajita = ~ 1), data = Acorticis, method
> "REML", control = lmc, weights = varComb (varIdent (form = ~ 1 |
fcajita),
> varPower (form = ~ Time10)))
>
>
>
> The first three models executed perfectly and the following values ??of
> the AIC indicator were obtained:
>
> > AIC (M.1 $ lme, M.2 $ lme, M.3 $ lme)
>
>
>
>         df       AIC
>
>
>
> M.1     8        379.6464
>
>
>
> M.2    15        309.5736
>
>
>
> M.3     9        310.4828
>
>
>
> Unfortunately, the execution of the M.4 model failed and the following
> error message was obtained:
>
>
>
> Error in environment (attr (ret $ lme $ modelStruct $ varStruct,
> "formula")) <-. GlobalEnv:
>
> attempt to set an attribute on NULL
>
>
>
> A final model I tried was M5:
>
> M.5 <- gamm(f1, random = list(fcajita =~ 1), data = Acorticis, method
> "REML", control = lmc, weights = varComb(varIdent(form =
~1|fcajita),
> varExp(form =~ Time10|fcajita)))
>
> and this time I got the following error message:
>
> Error in lme.formula(y ~ X - 1, random = rand, data = strip.offset(mf),  :
>
>   nlminb problem, convergence error code = 1
>
>   message = function evaluation limit reached without convergence (9)
>
> Adem?s: Warning message:
>
> In logLik.reStruct(object, conLin) :
>
>   Singular precision matrix in level -1, block 1
>
>
>
> My question is: Could someone help me fix these problems to run the M.4
> and M.5 models?
>
>
>         [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
	[[alternative HTML version deleted]]

> There was clearly a non-linear growth pattern such that an additive mixed
effects model was proposed to model the behavior of biomass as a function of
time and treatments.
The presence of non-linearity (in x) does not necessarily mean that
you can't use a "linear model".
In general, linear models work well for biological growth curves with
small to medium sized data.

Technically, this is off topic, however, the notion of what is linear
and what is not linear, is a recurring theme on this mailing list.
> When plotting the residuals a clear cone-shaped pattern was observed
Based on classical statistics, the obvious approach is to transform your data.
This is easy to do, however, some expertise is required to
back-transform parameters/estimates and interpret them.

While you're welcome to seek advice on mailing lists, I'd recommend
you consult a biostatistician.

On 29/02/20 8:11 am, Abby Spurdle wrote:
>> There was clearly a non-linear growth pattern such that an additive
>> mixed effects model was proposed to model the behavior of biomass
>> as a function of time and treatments.
> 
> The presence of non-linearity (in x) does not necessarily mean that 
> you can't use a "linear model". In general, linear models
work well
> for biological growth curves with small to medium sized data.
> 
> Technically, this is off topic, however, the notion of what is
> linear and what is not linear, is a recurring theme on this mailing
> list.
> 
>> When plotting the residuals a clear cone-shaped pattern was
>> observed
> 
> Based on classical statistics, the obvious approach is to transform
> your data. This is easy to do, however, some expertise is required
> to back-transform parameters/estimates and interpret them.
Although it is indeed off topic I would like to emphasise a point that a 
lot of people don't seem to get.  A linear model is

     *linear in the parameters*

not in the predictors.  For example

     y = beta_0 + beta_1 * x + beta_2 * x^2 + E

is a linear model.  It is linear in (beta_0,beta_1,beta_2).  It is not 
of course "linear in x".

cheers,

Rolf Turner

-- 
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

Sorry, I need to apologize.

My statement that "In general, linear models work well for biological
growth curves with small to medium sized data", may not be correct.
Also, I read through you code too quickly, my bad...

If you want to set weights (not sure if that's a good idea or not
here), I'd recommend you set them to a ***numeric vector***.

Your code is trying to set the weights to more complex objects (from
the nlme package), and it's possible that the gamm function is trying
to coerce them to numeric vectors, and that may be causing problems.

For nontrivial weights, you may want to compute weights (ensuring that
they are a numeric vector) in a separate step, and then fit the model
in a second step.


On 2/28/20, Jos? Antonio Garc?a P?rez <garci95 at hotmail.com>
wrote:
> I conducted an experiment where earthworms were subjected to two
treatments,
> with and without herbicide in the soil. Biomass measurements were taken
> every 12 days for 398 days and the biomass growth curves as a function of
> time were plotted.
>
> There was clearly a non-linear growth pattern such that an additive mixed
> effects model was proposed to model the behavior of biomass as a function
of
> time and treatments.
>
> When plotting the residuals a clear cone-shaped pattern was observed,
> therefore a series of additive models were proposed sequentially to deal
> with violations of the assumption of homogeneity. Below we can see the
> models with the following names: M.1; M.2; M.3; M.4
>
>
>
> lmc <- lmeControl (niterEM = 5000, msMaxIter = 1000)
>
> f1 <- formula (Biomass ~ Treat + s (Time, by = Treat))
>
>
>
> M.1 <-gamm (f1, random = list (fcajita = ~ 1), method =
"REML", control > lmc, data = Acorticis)
>
>
>
> #This first model uses the experimental box factor (i.e. fcajita) as the
> random element of the model. This random effects model assumes homogeneity
> between the experimental boxes and within them over time
>
>
>
> M.2 <-gamm (f1, random = list (fcajita = ~ 1), method =
"REML", control > lmc, data = Acorticis, weights = varIdent (form =
~ 1 | fcajita))
>
>
>
> #This second model assumes heterogeneity between boxes, but homogeneity
> within each box over time
>
>
>
> M.3 <- gamm (f1, random = list (fcajita = ~ 1), method =
"REML", control > lmc, data = Acorticis, weights = varExp (form = ~
Time10))
>
>
>
> #The third model assumes homogeneity between boxes but heterogeneity within
> each box over time
>
>
>
> Finally, we decided to model the heterogeneity using the 'varComb'
function
> in order to combine the variances where the model allows heterogeneity
> between the experimental boxes and heterogeneity within the experimental
> boxes over time:
>
>
>
> M.4 <- gamm (f1, random = list (fcajita = ~ 1), data = Acorticis, method
> "REML", control = lmc, weights = varComb (varIdent (form = ~ 1 |
fcajita),
> varPower (form = ~ Time10)))
>
>
>
> The first three models executed perfectly and the following values ??of the
> AIC indicator were obtained:
>
>> AIC (M.1 $ lme, M.2 $ lme, M.3 $ lme)
>
>
>
>         df       AIC
>
>
>
> M.1     8        379.6464
>
>
>
> M.2    15        309.5736
>
>
>
> M.3     9        310.4828
>
>
>
> Unfortunately, the execution of the M.4 model failed and the following
error
> message was obtained:
>
>
>
> Error in environment (attr (ret $ lme $ modelStruct $ varStruct,
"formula"))
> <-. GlobalEnv:
>
> attempt to set an attribute on NULL
>
>
>
> A final model I tried was M5:
>
> M.5 <- gamm(f1, random = list(fcajita =~ 1), data = Acorticis, method
> "REML", control = lmc, weights = varComb(varIdent(form =
~1|fcajita),
> varExp(form =~ Time10|fcajita)))
>
> and this time I got the following error message:
>
> Error in lme.formula(y ~ X - 1, random = rand, data = strip.offset(mf),  :
>
>   nlminb problem, convergence error code = 1
>
>   message = function evaluation limit reached without convergence (9)
>
> Adem?s: Warning message:
>
> In logLik.reStruct(object, conLin) :
>
>   Singular precision matrix in level -1, block 1
>
>
>
> My question is: Could someone help me fix these problems to run the M.4 and
> M.5 models?
>
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>