用R中的人口增长曲线模型提取个体增长常数

Question

I想直接从我们的增长模型中推导出个人增长率，类似于这个操作和这个操作。

我正在使用一个数据集，该数据集包含人口中大约2000个人的age和体重(wt)测量值。每个个体都由一个唯一的id编号表示。

数据的样本可以找到这里。数据如下所示：

          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

我开发了一个非线性增长曲线来模拟这个数据集的增长(在人口级别)。看起来是这样的：

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

它预测给定age的权重(age)，并具有可修改的参数A、t0和m。我使用nlme回归拟合将该模型与总体级别的数据集进行拟合，其中我将单个id指定为随机效应，并使用pdDiag将每个参数指定为不相关的。(注:在观察个人水平时，需要降低随机效应。)

此代码如下所示：

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

我有我们的人口水平增长模型(nlme.k)，但是我想用它来导出/提取每个增长常数的单个值。

如何使用我的人口增长模型id (**nlme.k**)? 为每个提取单个增长常数?注意，我不需要它是一个使用nlme的解决方案，这就是我用于人口增长模型的模型。

如有任何建议，将不胜感激！

Answer 1

我认为这是不可能的，因为自然如何设计随机效应。根据这个职位，使用部分池估计效果大小(您的增长常数)。这涉及到使用来自其他组的数据点。因此，您无法估计每个组的效果大小(您的个人id)。严格地说(参见这里)，随机效应实际上并不是模型的一部分，而是误差的一部分。

但是，您可以一起估计所有组的R2。如果您希望它在单个级别上(例如，id 1的参数估计)，那么只需在该特定个体的所有数据点上运行相同的模型即可。这给出n个参数集的n个模型。

Answer 2

最后，我们使用了几个循环来完成这个任务。

注意，如果有人想要后台脚本，我们的答案是建立在这个操作中发布的模型的基础上的。我们还将链接到已发布的脚本，当它被发布。

就目前而言，这应该给我们一个关于我们是如何做到这一点的总体想法。

#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

它提供了这样的输出：

> 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