[R] NLS -- multiplicative errors and group comparison

From: Derek Ogle <DOgle_at_northland.edu>
Date: Tue, 26 Feb 2008 17:58:33 -0600


Hello,  

I am attempting to fit a non-linear model (Von Bertalanffy growth model) to fish length-at-age data with the purpose of determining if any of the three parameters differ between male and female fish. I believe that I can successfully accomplish this goal assuming an additive error structure (illustrated in section 1 below). My trouble begins when I attempt this analysis using a model with a multiplicative error structure. I believe that I can fit the multiplicative error model to data without attempting to test between males and females (section 2 below), but I am not sure how to declare the model formula if I want to test between males and females (section 3 below). In addition, I assumed what I did in section 2 was correct, separated the data into a dataframe of males and a dataframe of females, and then attempted to fit separate models to each group (Section 4). This was "successful" for males but not for females. The problem with the female group appears to be that my "logging" creates NaN cells. When I put a trace on the nls function it appears to do one iteration and then fails. My assumption is that the nls function must be attempting different parameter estimates that return a NaN but I cannot test this as I'm not sure how to see the new parameters that it is trying. [I tried different starting values but did not find any that corrected this problem.]  

So (1) how do I declare the model formula for using multiplicative errors (i.e., is what I did in Section 2 correct), (2) how do I declare the model formula for comparing two groups and using multiplicative errors, and (3) any suggestions for how to find the "issue" leading to the error for just the female model in section 4.  

Sys.info information is in Section 4 and I am using R 2.6.1. I have not included the data as it is a rather large file.  

Thank you very much for reading this long message and for any help that you can offer.      

### Section 1 ###

> library(xlsReadWrite)

> fwd <- read.xls("FWDrum_SDF.xls",sheet=1)

> fwd$Gm <- fwd$Gf <- rep(0,length(fwd$age))

> fwd$Gf[fwd$sex=="female"] <- 1

> fwd$Gm[fwd$sex=="male"] <- 1

> str(fwd)

'data.frame': 719 obs. of 5 variables:

 $ age: num 1.27 2.25 2.25 2.25 2.25 ...

 $ tl : num 226 208 226 226 234 ...

 $ sex: Factor w/ 2 levels "female","male": 1 1 1 1 1 1 1 1 1 1 ...

 $ Gf : num 1 1 1 1 1 1 1 1 1 1 ...

 $ Gm : num 0 0 0 0 0 0 0 0 0 0 ...  

# starting values

> sLinf <- 405

> sK <- 0.11

> sto <- -5.2

> svb <- list(Linf=sLinf,K=sK,to=sto)
 

# Fit the additive error structure model to both groups combined

> vbla <- nls(tl~Linf*(1-exp(-K*(age-to))),start=svb,data=fwd)

> summary(vbl1)
 

Formula: tl ~ Linf * (1 - exp(-K * (age - to)))  

Parameters:

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

Linf 520.751048 19.560552 26.623 < 2e-16 ***

K 0.061097 0.008506 7.183 1.72e-12 ***

to -6.950346 1.112839 -6.246 7.24e-10 ***

---

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

 

Residual standard error: 83.75 on 716 degrees of freedom

 

Number of iterations to convergence: 9 

Achieved convergence tolerance: 7.525e-06 

 

# Fit the additive error structure model to the separate groups

#  compare this to vbla to see if any parameters differ (fit other
models to see which parameter)


> vbla.gen<- nls(tl~Gf*Linff*(1-exp(-Kf*(age-tof))) +
Gm*Linfm*(1-exp(-Km*(age-tom))), + start=list(Linff=sLinf,Kf=sK,tof=sto,Linfm=sLinf,Km=sK,tom=sto),data=fwd )
> summary(vbla.gen)
Formula: tl ~ Gf * Linff * (1 - exp(-Kf * (age - tof))) + Gm * Linfm * (1 - exp(-Km * (age - tom))) Parameters: Estimate Std. Error t value Pr(>|t|) Linff 692.719734 47.191732 14.679 < 2e-16 *** Kf 0.038532 0.006469 5.957 4.04e-09 *** tof -7.589941 1.235644 -6.142 1.35e-09 *** Linfm 372.224132 13.268491 28.053 < 2e-16 *** Km 0.095814 0.024983 3.835 0.000137 *** tom -8.058334 2.445654 -3.295 0.001033 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 68.15 on 713 degrees of freedom Number of iterations to convergence: 8 Achieved convergence tolerance: 6.589e-06 ### Section 2 ###
> vblm <- nls(log(tl)~log(Linf*(1-exp(-K*(age-to)))),start=svb,data=fwd)

> summary(vblm)
Formula: log(tl) ~ log(Linf * (1 - exp(-K * (age - to)))) Parameters: Estimate Std. Error t value Pr(>|t|) Linf 515.414166 21.239627 24.267 < 2e-16 *** K 0.055431 0.007486 7.404 3.71e-13 *** to -8.289902 1.002402 -8.270 6.51e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.1907 on 716 degrees of freedom Number of iterations to convergence: 8 Achieved convergence tolerance: 2.236e-06 ### Section 3 ### # FIRST TRY -- Fit the multiplicative error structure model to the separate groups
> vblm.gen1<- nls(log(tl)~log(Gf*Linff*(1-exp(-Kf*(age-tof))) +
Gm*Linfm*(1-exp(-Km*(age-tom)))), start=list(Linff=sLinf,Kf=sK,tof=sto,Linfm=sLinf,Km=sK,tom=sto),data=fwd ) Error in numericDeriv(form[[3]], names(ind), env) : Missing value or an infinity produced when evaluating the model In addition: Warning message: In log(Gf * Linff * (1 - exp(-Kf * (age - tof))) + Gm * Linfm * : NaNs produced # SECOND TRY -- Fit the multiplicative error structure model to the separate groups
> vblm.gen2<- nls(log(tl)~log(Gf*Linff*(1-exp(-Kf*(age-tof)))) +
log(Gm*Linfm*(1-exp(-Km*(age-tom)))), start=list(Linff=sLinf,Kf=sK,tof=sto,Linfm=sLinf,Km=sK,tom=sto),data=fwd ) Error in numericDeriv(form[[3]], names(ind), env) : Missing value or an infinity produced when evaluating the model
> THIRD TRY -- Fit the multiplicative error structure model to the
separate groups
> vblm.gen3<- nls(log(tl)~Gf*log(Linff*(1-exp(-Kf*(age-tof)))) +
Gm*log(Linfm*(1-exp(-Km*(age-tom)))), start=list(Linff=sLinf,Kf=sK,tof=sto,Linfm=sLinf,Km=sK,tom=sto),data=fwd ) Error in numericDeriv(form[[3]], names(ind), env) : Missing value or an infinity produced when evaluating the model In addition: Warning message: In log(Linff * (1 - exp(-Kf * (age - tof)))) : NaNs produced ### Section 4 ### # Just males
> fwd.m <- subset(fwd,sex=="male")

> vblm.m <-
nls(log(tl)~log(Linf*(1-exp(-K*(age-to)))),start=svb,data=fwd.m)
> summary(vblm.m)
Formula: log(tl) ~ log(Linf * (1 - exp(-K * (age - to)))) Parameters: Estimate Std. Error t value Pr(>|t|) Linf 363.07557 4.88130 74.381 < 2e-16 *** K 0.11765 0.01184 9.936 < 2e-16 *** to -6.29861 0.75724 -8.318 4.15e-15 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.09035 on 276 degrees of freedom Number of iterations to convergence: 5 Achieved convergence tolerance: 7.262e-06 # Just females
> fwd.f <- subset(fwd,sex=="female")

> vblm.f <-
nls(log(tl)~log(Linf*(1-exp(-K*(age-to)))),start=svb,data=fwd.f) 26.89721 : 405.00 0.11 -5.20 Error in numericDeriv(form[[3]], names(ind), env) : Missing value or an infinity produced when evaluating the model In addition: Warning message: In log(Linf * (1 - exp(-K * (age - to)))) : NaNs produced ### Section 5 ###
> Sys.info()
sysname release version "Windows" "NT 5.1" "(build 2600) Service Pack 2" nodename machine login "CSE229-001" "x86" [[alternative HTML version deleted]] ______________________________________________ R-help_at_r-project.org mailing list 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.
Received on Wed 27 Feb 2008 - 00:02:28 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 27 Feb 2008 - 00:30:16 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive