1

I am attempting to fit an analytical model to my observed data on human somatic growth. The generalized model is a 3-parameter logarithmic growth curve where:

s= B0+B1*t+B2*log(t),

where s is a growth parameter whether it be a bone length or stature and t is an age.

I have attempted to run through various iterations of code to both form a likelihood function and and maximize teh returned parameters. To be completely honest, I am totally unsure if I am writing this correctly, but the most recent iteration of my attempts looks as such using a dataframe title "cedar":

cedar.f<-function(b){sum((cedar$FLEN~b[0]+b[1]*cedar$MINAGE+b[2]*log(cedar$MINAGE))^2)}
cedar.optimx<-optimx(c(0,0,0), fn = cedar.f, control = list(all.methods=TRUE, save.failures=TRUE, maxit=5000, maximize=TRUE))

cedar$MINAGE=c(2.5,0.5,6,0.01,0.01,0.01,1,0.01,12,0.01,0.01,1,12,1,4.5,1,4.5,0.01,7.8,11,4,7.5,9,0.25,0.01,0.01,0.01,0.1,1,1,0.01,0.01)
cedar$FLEN=c(167,150,300,54,60,78,152,72, 452,84,81,163,450,137,271,151,261,73,330,378,226,319,378,89,78,89,65,81,144,122,  50,55

Above, I have provided my attempt at the function and the optimization. I have received several errors in my attempts at doing this including:

argument "cedar" is missing ,no default
non-numeric argument to binary operator

Really, I am here to ask what recommendations anyone may have in writing this function so as to maximize the best fit of the data to the analytical human growth curve. If I am going about this all wrong / writing the function wrong, any help would be appreciated.

Thank you all.

Chris Wolfe
  • 11
  • 2

1 Answers1

1
cedar <- data.frame(MINAGE=c(2.5,0.5,6,0.01,0.01,0.01,1,0.01,12,0.01,0.01,1,12,1,4.5,1,4.5,0.01,7.8,11,4,7.5,9,0.25,0.01,0.01,0.01,0.1,1,1,0.01,0.01),
         FLEN=c(167,150,300,54,60,78,152,72,452,84,81,163,450,137,271,151,261,73,330,378,226,319,378,89,78,89,65,81,144,122,50,55))

# Sum of squared errors
# Minus sign is for function minimization    
cedar.f <- function(b) { 
    -sum( (cedar$FLEN - (b[1] + b[2]*cedar$MINAGE + b[3]*log(cedar$MINAGE)))^2 )
}

library(optimx)
cedar.optimx <- optimx( c(1,1,1), fn = cedar.f, 
    control = list(all.methods=TRUE, save.failures=TRUE, maxit=5000, maximize=TRUE))

#                 p1       p2       p3          value fevals gevals niter convcode kkt1 kkt2 xtimes
# BFGS        120.4565 24.41910 11.25419  -7.674935e+03     25      8    NA        0 TRUE TRUE   0.00
# CG          120.4565 24.41910 11.25419  -7.674935e+03   1072    298    NA        0 TRUE TRUE   0.15
# Nelder-Mead 120.4714 24.41647 11.25186  -7.674947e+03    258     NA    NA        0 TRUE TRUE   0.02
# L-BFGS-B    120.4565 24.41910 11.25419  -7.674935e+03     17     17    NA        0 TRUE TRUE   0.01
# nlm         120.4564 24.41910 11.25417  -7.674935e+03     NA     NA    12        0 TRUE TRUE   0.01
# nlminb      120.4565 24.41910 11.25419  -7.674935e+03     21     48    13        0 TRUE TRUE   0.02
# spg         120.4565 24.41910 11.25419  -7.674935e+03     99     NA    92        0 TRUE TRUE   0.06
# ucminf      120.4564 24.41910 11.25417  -7.674935e+03     10     10    NA        0 TRUE TRUE   0.00
# Rcgmin            NA       NA       NA -8.988466e+307     NA     NA    NA     9999   NA   NA   0.00
# Rvmmin            NA       NA       NA -8.988466e+307     NA     NA    NA     9999   NA   NA   0.00
# newuoa      120.4565 24.41910 11.25419  -7.674935e+03    118     NA    NA        0 TRUE TRUE   0.01
# bobyqa      120.4565 24.41910 11.25419  -7.674935e+03    142     NA    NA        0 TRUE TRUE   0.02
# nmkb        120.4566 24.41907 11.25421  -7.674935e+03    213     NA    NA        0 TRUE TRUE   0.03
# hjkb          1.0000  1.00000  1.00000  -1.363103e+06      1     NA     0     9999   NA   NA   0.00

Alternatively, model coefficients can be estimated using a simple linear model:

fitlm <- lm(FLEN~MINAGE+log(MINAGE), data=cedar)
coef(fitlm)

# Intercept)      MINAGE log(MINAGE) 
#  120.45654    24.41910    11.25419

The estimated function can be plotted as follows:

optpar <- as.matrix(cedar.optimx[1,1:3])
estim_fun <- function(x, b=optpar) { 
    b[1] + b[2]*x + b[3]*log(x)
}
curve(estim_fun, from=min(cedar$MINAGE), to=max(cedar$MINAGE))

enter image description here

Marco Sandri
  • 23,289
  • 7
  • 54
  • 58
  • Thank you so much! I have been attempting for quite some time to find the best way to maximize the parameters. Was I thinking beyond the scope of my question by using an MLE versus a "typical" linear/nonlinear model? Also, quick follow up, what if age is not a discrete parameter, but instead is a distribution with a mean and upper and lower bounds. For example, the data is actually formatted with a MINAGE, MEANAGE, and MAXAGE. Is there a way to design a function that encapsulates this error range? – Chris Wolfe Sep 30 '17 at 17:47
  • What is the best way to plot this curve given the current parameters? – Chris Wolfe Sep 30 '17 at 22:51
  • @ChrisWolfe About plotting the estimated function: I added some code to my answer. About your problem with MINAGE, MEANAGE, and MAXAGE: I suggest to open a new question on https://stats.stackexchange.com/ where you explain in detail your problem. – Marco Sandri Oct 01 '17 at 08:31
  • Thank you so much. Will do. – Chris Wolfe Oct 01 '17 at 15:13