如何在變量數量發生變化時在循環中編寫公式？

Question

我使用下面的代碼來爲每組觀察值擬合不同的指數曲線，它工作得很好。

p = c(10,20,15,25,20,30,25,35,30,40,25,35,20,30,15,25,10,20) 
v = c(92,110,104,117,123,139,146,162,165,176,160,176,143,163,118,137,92,110) 
group = factor(rep((1:9), each=2)) 

mm = model.matrix(~ 0 + group) 

fit = nls(v ~ drop(mm %*% c(b1, b2, b3, b4, b5, b6, b7, b8, b9))*(1-exp(-k*p)), 
     start = list(k=0.5, b1=1000, b2=2000, b3=3000, b4=4000, b5=5000, b6=6000, b7=7000, b8=8000, b9=9000)) 

summary(fit) 
Formula: v ~ drop(mm %*% c(b1, b2, b3, b4, b5, b6, b7, b8, b9)) * (1 - 
exp(-k * p)) 

Parameters: 
Estimate Std. Error t value Pr(>|t|)  
k 0.10928 0.01374 7.954 4.55e-05 *** 
b1 129.13042 8.01108 16.119 2.20e-07 *** 
b2 126.81086 6.43352 19.711 4.57e-08 *** 
b3 141.74666 5.62817 25.185 6.61e-09 *** 
b4 161.10250 5.06762 31.791 1.04e-09 *** 
b5 174.94417 4.63884 37.713 2.68e-10 *** 
b6 175.73100 5.20655 33.752 6.48e-10 *** 
b7 165.58007 6.01256 27.539 3.26e-09 *** 
b8 146.48962 6.92337 21.159 2.62e-08 *** 
b9 129.13042 8.01108 16.119 2.20e-07 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 5.767 on 8 degrees of freedom 

Number of iterations to convergence: 8 
Achieved convergence tolerance: 5.907e-06 

我需要爲不同的實驗組（在for循環中）多次運行此代碼。棘手的部分是，對於每個實驗，預測因子的數量（即b1，b2等）改變。

有沒有簡單的出路？

Answer 1

您可以使用as.formula從字符串構建公式。

start <- as.list(c(k=0.5, setNames(1:9*1000, paste0("b", 1:9)))) # start values 
groups <- list(g1=1:4, g2=2:4, g3=c(1,6:9))      # some random groupings 

lapply(groups, function(group){         # apply model to groups 
    mm <- mm[, group] 
    groupStr <- sprintf("b%d", group) 
    nls(
     as.formula(paste("v~drop(mm%*%c(", paste(groupStr, collapse=","), "))*(1-exp(-k*p))")), 
     start=start[c("k", groupStr)] 
    ) 
}) 

Answer 2

您有一個部分線性最小二乘模型。使用適當的算法，不僅不需要爲組指定參數，還可以避免使用起始值。

fit1 <- nls(v ~ mm * (1 - exp(-k * p)), algorithm = "plinear", 
      start = list(k=0.5)) 
summary(fit1) 

#Formula: v ~ mm * (1 - exp(-k * p)) 
# 
#Parameters: 
#    Estimate Std. Error t value Pr(>|t|)  
#k    0.10928 0.01374 7.954 4.55e-05 *** 
#.lin.group1 129.13036 8.01107 16.119 2.20e-07 *** 
#.lin.group2 126.81082 6.43351 19.711 4.57e-08 *** 
#.lin.group3 141.74664 5.62816 25.185 6.61e-09 *** 
#.lin.group4 161.10248 5.06761 31.791 1.04e-09 *** 
#.lin.group5 174.94416 4.63883 37.713 2.68e-10 *** 
#.lin.group6 175.73098 5.20655 33.752 6.48e-10 *** 
#.lin.group7 165.58004 6.01256 27.539 3.26e-09 *** 
#.lin.group8 146.48959 6.92337 21.159 2.62e-08 *** 
#.lin.group9 129.13036 8.01107 16.119 2.20e-07 *** 
#--- 
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
# 
#Residual standard error: 5.767 on 8 degrees of freedom 
# 
#Number of iterations to convergence: 9 
#Achieved convergence tolerance: 6.958e-06