Extracting individual growth constants using population growth curve model in R

Question

I would like to derive individual growth rates from our growth model directly, similar to this OP and this OP.

I am working with a dataset that contains the age and weight (wt) measurements for ~2000 individuals in a population. Each individual is represented by a unique id number.

A sample of the data can be found here. Here is what the data looks like:

          id   age    wt
        1615     6      15  
        3468     32     61  
        1615     27     50  
        1615     60     145    
        6071     109    209  
        6071     125    207  
       10645     56     170   
       10645     118    200

I have developed a non-linear growth curve to model growth for this dataset (at the population level). It looks like this:

wt~ A*atan(k*age - t0) + m

which predicts weight (wt) for a given age and has modifiable parameters A, t0, and m. I have fit this model to the dataset at the population level using a nlme regression fit where I specified individual id as a random effect and used pdDiag to specify each parameter as uncorrelated. (Note: the random effect would need to be dropped when looking at the individual level.)

The code for this looks like:

nlme.k = nlme(wt~ A*atan(k*age - t0) + m, 
data = df, 
fixed = A+k+t0+m~1, 
random = list(id = pdDiag(A+t0+k+m~1)), #cannot include when looking at the individual level
start = c(A = 99.31,k = 0.02667, t0 = 1.249, m = 103.8), #these values are what we are using at the population level # might need to be changed for individual models
na.action = na.omit,
control = nlmeControl(maxIter = 200, pnlsMaxIter = 10, msMaxIter = 100))

I have our population level growth model (nlme.k), but I would like to use it to derive/extract individual values for each growth constant.

How can I extract individual growth constants for each id using my population level growth model (nlme.k)? Note that I don't need it to be a solution that uses nlme, that is just the model I used for the population growth model.

Any suggestions would be appreciated!

Answer 1

We ended up using a few loops to do this.

Note that our answer builds off a model posted in this OP if anyone wants the background script. We will also link to the published script when it is posted.

For now - this is should give a general idea of how we did this.

#Individual fits dataframe generation
yid_list  <- unique(young_inds$squirrel_id)
indf_prs <- list('df', 'squirrel_id', 'A_value', 'k_value', 'mx_value', 'my_value', 'max_grate', 'hit_asymptote', 'age_asymptote', 'ind_asymptote', 'ind_mass_asy', 'converge') #List of parameters 
ind_fits <- data.frame(matrix(ncol = length(indf_prs), nrow = length(yid_list))) #Blank dataframe for all individual fits 
colnames(ind_fits) <- indf_prs    

#Calculates individual fits for all individuals and appends into ind_fits 

for (i in 1:length(yid_list)) {

yind_df <-young_inds%>%filter(squirrel_id %in% yid_list[i]) #Extracts a dataframe for each squirrel 
ind_fits[i , 'squirrel_id'] <- as.numeric(yid_list[i]) #Appends squirrel i's id into individual fits dataframe 
sex_lab  <- unique(yind_df$sex) #Identifies and extracts squirrel "i"s sex 
mast_lab <- unique(yind_df$b_mast) #Identifies and extracts squirrel "i"s mast value 

Hi_dp <- max(yind_df$wt) #Extracts the largest mass for each squirrel 
ind_long <- unique(yind_df$longevity) #Extracts the individual death date 

#Sets corresponding values for squirrel "i" 
if (mast_lab==0 && sex_lab=="F") { #Female no mast 
    ind_fits[i , 'df'] <- "fnm" #Squirrel dataframe (appends into ind_fits dataframe)
    df_asm   <- af_asm  #average asymptote value corresponding to sex 
    df_B_guess   <- guess_df[1, "B_value"]  #Inital guesses for nls fits corresponding to sex and mast sex and mast 
    df_k_guess   <- guess_df[1, "k_value"] 
    df_mx_guess  <- guess_df[1, "mx_value"] 
    df_my_guess  <- guess_df[1, "my_value"]
    ind_asyr     <- indf_asy  #growth rate at individual asymptote 

} else if (mast_lab==0 && sex_lab=="M") { #Male no mast 
    ind_fits[i , 'df'] <- "mnm"  
    df_asm <- am_asm
    df_B_guess   <- guess_df[2, "B_value"]  
    df_k_guess   <- guess_df[2, "k_value"] 
    df_mx_guess  <- guess_df[2, "mx_value"] 
    df_my_guess  <- guess_df[2, "my_value"]
    ind_asyr     <- indm_asy 

} else if (mast_lab==1 && sex_lab=="F") { #Female mast
    ind_fits[i , 'df'] <- "fma" 
    df_asm <- af_asm
    df_B_guess   <- guess_df[3, "B_value"]  
    df_k_guess   <- guess_df[3, "k_value"] 
    df_mx_guess  <- guess_df[3, "mx_value"] 
    df_my_guess  <- guess_df[3, "my_value"]
    ind_asyr     <- indm_asy 

} else if (mast_lab==1 && sex_lab=="M") { #Males mast 
    ind_fits[i , 'df'] <- "mma"
    df_asm <- am_asm 
    df_B_guess   <- guess_df[4, "B_value"]  
    df_k_guess   <- guess_df[4, "k_value"] 
    df_mx_guess  <- guess_df[4, "mx_value"] 
    df_my_guess  <- guess_df[4, "my_value"]
    ind_asyr     <- indf_asy 

} else { #If sex or mast is not identified or identified improperlly in the data
    print("NA") 

} #End of if else loop


#Arctangent 

#Fits nls model to the created dataframe
nls.floop <- tryCatch({data.frame(tidy(nls(wt~ B*atan(k*(age - mx)) + my, #tryCatch lets nls have alternate results instead of "code stopping" errors
    data=yind_df, 
    start = list(B = df_B_guess, k = df_k_guess, mx = df_mx_guess, my = df_my_guess),
    control= list(maxiter = 200000, minFactor = 1/100000000))))
},
error = function(e){  
    nls.floop <- data.frame(c(0,0), c(0,0)) #Specifies nls.floop as a dummy dataframe if no convergence
},
warning = function(w) {
    nls.floop <- data.frame(tidy(nls.floop)) #Fit is the same if warning is displayed   
}) #End of nls.floop

#Creates a dummy numerical index from nls.floop for if else loop below 
numeric_floop <-  as.numeric(nls.floop[1, 2])
#print(numeric_floop) #Taking a look at the values. If numaric floop...  
# == 0, function did not converge on iteration "i"
# != 0, function did converge on rapid "i" and code will run through calculations  

if (numeric_floop != 0) {

results_DF <- nls.floop

ind_fits[i , 'converge'] <- 1 #converge = 1 for converging fit 

#Extracting, calculating, and appending values into dataframe  
B_value   <- as.numeric(results_DF[1, "estimate"])  #B value 
k_value   <- as.numeric(results_DF[2, "estimate"])  #k value
mx_value  <- as.numeric(results_DF[3, "estimate"])  #mx value 
my_value  <- as.numeric(results_DF[4, "estimate"])  #my value 
A_value <- ((B_value*pi)/2)+ my_value  #A value calculation
ind_fits[i , 'A_value']  <- A_value 
ind_fits[i , 'k_value']  <- k_value
ind_fits[i , 'mx_value'] <- mx_value
ind_fits[i , 'my_value'] <- my_value #appends my_value into df
ind_fits[i , 'max_grate']     <-  adr(mx_value, B_value, k_value, mx_value, my_value) #Calculates max growth rate

}    
    } #End of individual fits loop

Which gives this output:

> head(ind_fits%>%select(df, squirrel_id, A_value, k_value, mx_value, my_value))
   df squirrel_id  A_value    k_value mx_value  my_value
1 mnm         332 257.2572 0.05209824 52.26842 126.13183
2 mnm        1252 261.0728 0.02810033 42.37454 103.02102
3 mnm        3466 260.4936 0.03946594 62.27705 131.56665
4 fnm         855 437.9569 0.01347379 86.18629 158.27641
5 fnm        2409 228.7047 0.04919819 63.99252 123.63404
6 fnm        1417 196.0578 0.05035963 57.67139  99.65781

Note that you need to create a blank dataframe first before running the loops.

Answer 2

I think this is not possible due to the nature on how random effects are designed. According to this post the effect size (your growth constant) is estimated using partial pooling. This involves using data points from other groups. Thus you can not estimate the effect size of each group (your individual id). Strictly speaking (see here) random effects are not really a part of the model at all, but more a part of the error.

However, you can estimate the R2 for all groups together. If you want it on an individual level (e.g. parameter estiamtes for id 1), then just run the same model only on all data points of this particular individual. This give you n models with n parameter sets for n individuals.