从nls模型中提取循环模型参数

Question

我们有这个for loop

tenIDs <- unique(ten_squirrel$squirrel_id)
tenIDs
[1] 5241 3102 2271 3119 3216

#for loop to run through all the squirrels

for (i in tenIDs){
#Creating our dataframe for squirrel_id "i"
Individual_DFs <- ten_squirrel %>% filter (squirrel_id %in% i) 

#Fit model for squirrel_id "i"
nls.floop <- nls(wt ~ A*atan(k*age - t0) + m, 
    data = Individual_DFs, 
    start = list(A = 102.8, k = 0.02, t0 = 0.751, m = 82.06))

#Show resulting output
print(nls.floop) 
}

它提供了这样的输出：

Nonlinear regression model
  model: wt ~ A * atan(k * age - t0) + m
   data: Individual_DFs
        A         k        t0         m 
100.69638   0.03493   1.78392 123.87479 
 residual sum-of-squares: 401.1

Number of iterations to convergence: 7 
Achieved convergence tolerance: 5.341e-06
Nonlinear regression model
  model: wt ~ A * atan(k * age - t0) + m
   data: Individual_DFs
        A         k        t0         m 
140.23662   0.01953   0.54546  63.33266 
 residual sum-of-squares: 215.8

Number of iterations to convergence: 12 
Achieved convergence tolerance: 4.367e-06
Nonlinear regression model
  model: wt ~ A * atan(k * age - t0) + m
   data: Individual_DFs
        A         k        t0         m 
 70.76447   0.04409   2.04846 101.03060 
 residual sum-of-squares: 146

Number of iterations to convergence: 9 
Achieved convergence tolerance: 5.725e-06
Nonlinear regression model
  model: wt ~ A * atan(k * age - t0) + m
   data: Individual_DFs
       A        k       t0        m 
94.35265  0.03234  1.30053 96.80194 
 residual sum-of-squares: 199.5

Number of iterations to convergence: 8 
Achieved convergence tolerance: 7.996e-06
Nonlinear regression model
  model: wt ~ A * atan(k * age - t0) + m
   data: Individual_DFs
       A        k       t0        m 
75.60633  0.04844  2.06589 98.25557 
 residual sum-of-squares: 481.2

Number of iterations to convergence: 9 
Achieved convergence tolerance: 4.24e-06

我们只在5个模型中显示这一点，但是我想从每个模型中提取A、k、t0和m，这样我们就有了这样的输出：

 tenIDs  A           k          t0       m
 3216    75.60633    0.04844    2.06589  98.25557 
 3119    94.35264    0.03234    1.30053  96.80194
 2271    70.76447    0.04409    2.04846  101.03060 
 3102    140.23656   0.01953    0.54546  63.33272
 5241    100.69638   0.03493    1.78392  123.87479 

如何为每个模型以及相应的tenIDs**?**提取这些值？

Answer 1

模型的输出包含函数...$m$getPars()，得到拟合的参数。

如果没有示例数据集，我无法验证，但是您应该通过以下修改代码获得所需的输出：

tenIDs <- unique(ten_squirrel$squirrel_id)
tenIDs

# define start to obtain the number and name of fit parameters
start <- list(A = 102.8, k = 0.02, t0 = 0.751, m = 82.06)
# create empty data.frame to store IDs and parameters
params <- data.frame(matrix(nrow = length(tenIDs), ncol = 1+length(start)))
names(params) <- c("tenIDs",names(start))

#for loop to run through all the squirrels
# (changed to 'seq_along()' to get consecutive 'i's
#  for easier indexing of the 'param'-data.frame)
for (i in seq_along(tenIDs)){
  #Creating our dataframe for squirrel_id "i"
  Individual_DFs <- ten_squirrel %>% filter (squirrel_id %in% tenIDS[i]) 
  
  #Fit model for squirrel_id "i"
  nls.floop <- nls(wt ~ A*atan(k*age - t0) + m, 
                   data = Individual_DFs, 
                   start = start)
  
  # store IDs
  params[i,1] <- tenIDs[i]
  # store fit parameters
  params[i,2:ncol(params)] <- nls.floop$m$getPars()
}

Answer 2

最后，我们做了类似于上面建议的事情。如果有人想要后台脚本，我们的答案是建立在这个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