R语言利用种群生长曲线模型提取种群生长常数

yb3bgrhw 于 2022-12-06 发布在其他

关注(0)|答案(2)|浏览(224)

**我希望直接从我们的增长模型中推导出个人增长率，**类似于此OP和此OP。

我正在处理一个数据集，其中包含一个种群中大约2000个个体的age和体重（wt）测量值。每个个体都由一个唯一的id数字表示。
here是一个数据示例。

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的权重（wt），并具有可修改的参数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），但我想用它来推导/提取每个增长常数的 * 个体 * 值。

**如何使用我的人口水平增长模型（nlme.k）为每个id提取单独的增长常数？**请注意，我不需要它是使用nlme的解决方案，这只是我用于人口增长模型的模型。

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

r

来源：https://stackoverflow.com/questions/71039697/extracting-individual-growth-constants-using-population-growth-curve-model-in-r

2条答案

按热度按时间

iezvtpos1#

我认为这是不可能的，因为随机效应是如何设计的。根据this post，效应大小（您的增长常数）是使用 * 部分合并 * 估计的。这涉及使用来自其他组的数据点。因此，您无法估计每组的效应量（你的个体id）。严格地说（见here），随机效应根本不是模型的一部分，而更多的是误差的一部分。
但是，您可以同时估计所有组的R2。如果您希望在个体水平上进行估计（例如，id 1的参数估计），则只需对该特定个体的所有数据点运行相同的模型。这将为您给予n个模型，其中n个参数集对应n个个体。

赞(0）回复(0）举报 2022-12-06

jfewjypa2#

我们最终使用了几个循环来完成这一任务。
请注意，如果有人想要后台脚本，我们的答案是建立在OP中发布的模型之上的。
现在-这是应该给予我们如何做到这一点的一个大致的想法。

#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

赞(0）回复(0）举报 2022-12-06

我来回答

R语言利用种群生长曲线模型提取种群生长常数

2条答案

相关问题

热门标签

最新问答

R语言 利用种群生长曲线模型提取种群生长常数

2条答案

相关问题

热门标签

最新问答

R语言利用种群生长曲线模型提取种群生长常数