1

I use the following code to fit different exponential curves for each group of observations, and it works just fine. 

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

I need to run this code several times for different sets of experiments (within a for loop). The tricky part is that, for each experiment, the number of predictors (i.e. b1, b2, etc) changes.

Is there an easy way out?

Eduardo Costa
  • 11
  • 1

2 Answers2

0

You can use as.formula to build formulas from strings.

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)]
    )
})
Rorschach
  • 31,301
  • 5
  • 78
  • 129
0

You have a partially linear least-squares model. Use the appropriate algorithm and not only don't you have to specify the parameters for your groups, but you also avoid the need for starting values.

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
Roland
  • 127,288
  • 10
  • 191
  • 288