1

我想为Temp 和 Variety 组合的包中log-logistic的以下数据拟合回归。drc R但是，我的代码抛出以下错误

Error in Temp:Variety : NA/NaN argument

代码：

df2 <-
  structure(list(Temp = c(15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 
15L, 20L, 20L, 20L, 20L, 20L, 20L, 25L, 25L, 25L, 25L, 30L, 30L, 
30L, 30L, 35L, 35L, 35L, 35L, 40L, 40L, 40L, 40L, 15L, 15L, 15L, 
15L, 15L, 15L, 15L, 15L, 15L, 20L, 20L, 20L, 20L, 20L, 20L, 25L, 
25L, 25L, 25L, 30L, 30L, 30L, 30L, 35L, 35L, 35L, 35L, 40L, 40L, 
40L, 40L), Variety = c("FH-142", "FH-142", "FH-142", "FH-142", 
"FH-142", "FH-142", "FH-142", "FH-142", "FH-142", "FH-142", "FH-142", 
"FH-142", "FH-142", "FH-142", "FH-142", "FH-142", "FH-142", "FH-142", 
"FH-142", "FH-142", "FH-142", "FH-142", "FH-142", "FH-142", "FH-142", 
"FH-142", "FH-142", "FH-142", "FH-142", "FH-142", "FH-142", "FH-942", 
"FH-942", "FH-942", "FH-942", "FH-942", "FH-942", "FH-942", "FH-942", 
"FH-942", "FH-942", "FH-942", "FH-942", "FH-942", "FH-942", "FH-942", 
"FH-942", "FH-942", "FH-942", "FH-942", "FH-942", "FH-942", "FH-942", 
"FH-942", "FH-942", "FH-942", "FH-942", "FH-942", "FH-942", "FH-942", 
"FH-942", "FH-942"), Start = c(0L, 24L, 48L, 72L, 96L, 120L, 
144L, 168L, 192L, 0L, 24L, 48L, 72L, 96L, 120L, 0L, 24L, 48L, 
72L, 0L, 24L, 48L, 72L, 0L, 24L, 48L, 72L, 0L, 24L, 48L, 72L, 
0L, 24L, 48L, 72L, 96L, 120L, 144L, 168L, 192L, 0L, 24L, 48L, 
72L, 96L, 120L, 0L, 24L, 48L, 72L, 0L, 24L, 48L, 72L, 0L, 24L, 
48L, 72L, 0L, 24L, 48L, 72L), End = c(24, 48, 72, 96, 120, 144, 
168, 192, Inf, 24, 48, 72, 96, 120, Inf, 24, 48, 72, 96, 24, 
48, 72, Inf, 24, 48, 72, Inf, 24, 48, 72, Inf, 24, 48, 72, 96, 
120, 144, 168, 192, Inf, 24, 48, 72, 96, 120, Inf, 24, 48, 72, 
Inf, 24, 48, 72, Inf, 24, 48, 72, Inf, 24, 48, 72, Inf), Germinated = c(0L, 
0L, 1L, 3L, 3L, 12L, 14L, 12L, 15L, 0L, 11L, 27L, 15L, 3L, 4L, 
2L, 30L, 15L, 13L, 6L, 43L, 7L, 4L, 5L, 48L, 3L, 4L, 0L, 31L, 
21L, 8L, 0L, 0L, 0L, 12L, 13L, 6L, 2L, 1L, 26L, 0L, 10L, 13L, 
11L, 13L, 13L, 11L, 21L, 19L, 9L, 7L, 18L, 23L, 12L, 14L, 23L, 
12L, 11L, 12L, 18L, 11L, 19L)), .Names = c("Temp", "Variety", 
"Start", "End", "Germinated"), row.names = c(NA, -62L), class = "data.frame")


library(drc)
fm2 <-
  drm(
        formula   = Germinated ~ Start + End
      , curveid   = Temp:Variety
  #   , pmodels   = 
  #   , weights   = 
      , data      = df2
  #   , subset    = 
      , fct       = LL.2()
      , type      = "event"
      , bcVal     = NULL
      , bcAdd     = 0
  #   , start     =
      , na.action = na.fail
      , robust    = "mean"
      , logDose   = NULL
      , control   = drmc(
                            constr      = FALSE
                            , errorm      = TRUE
                            , maxIt       = 1500
                            , method      = "BFGS"
                            , noMessage   = FALSE
                            , relTol      = 1e-07
                            , rmNA        = FALSE
                            , useD        = FALSE
                            , trace       = FALSE
                            , otrace      = FALSE
                            , warnVal     = -1
                            , dscaleThres = 1e-15
                            , rscaleThres = 1e-15
                            )
      , lowerl    = NULL
      , upperl    = NULL
      , separate  = FALSE
      , pshifts   = NULL 
      )
4

1 回答 1

1

有几种方法可以做到这一点。一种是使用apply循环其中一个因素。但是，可能最简单和最简洁的方法是创建一个分组因子（让我们称之为grp），基于您想要分组的所有变量的独特组合。然后，您可以简单地使用该列对 drm 进行分组curveid=grp。

df2$grp <- factor(paste(df2$Temp, df2$Variety, sep = ","))
fm <-drm(data = df2, curveid = grp,
          formula = Germinated ~ Start + End,  fct = LL.2(), type = "event", 
          control = drmc(constr = FALSE, errorm = TRUE, maxIt = 1500, method = "BFGS", 
                         noMessage = FALSE, relTol = 1e-07, rmNA = FALSE, useD = FALSE, 
                         trace = FALSE, otrace = FALSE, warnVal = -1, dscaleThres = 1e-15, rscaleThres = 1e-15))


summary(fm)

Model fitted: Log-logistic (ED50 as parameter) with lower limit at 0 and upper limit at 1 (2 parms)

Parameter estimates:

             Estimate Std. Error   t-value p-value
b:15,FH-142  -6.03042    0.78914  -7.64179       0
b:20,FH-142  -4.96487    0.60745  -8.17333       0
b:25,FH-142  -4.43982    0.54905  -8.08634       0
b:30,FH-142  -4.80965    0.60804  -7.91014       0
b:35,FH-142  -5.45881    0.69145  -7.89477       0
b:40,FH-142  -5.43884    0.79770  -6.81819       0
b:15,FH-942  -2.92989    0.42582  -6.88058       0
b:20,FH-942  -3.32332    0.42392  -7.83949       0
b:25,FH-942  -2.89164    0.39282  -7.36123       0
b:30,FH-942  -3.23898    0.45011  -7.19602       0
b:35,FH-942  -2.43731    0.35084  -6.94698       0
b:40,FH-942  -1.95940    0.31456  -6.22908       0
e:15,FH-142 162.33560    6.10484  26.59130       0
e:20,FH-142  64.71878    3.08661  20.96757       0
e:25,FH-142  48.23889    2.68275  17.98115       0
e:30,FH-142  36.38342    2.04231  17.81486       0
e:35,FH-142  35.07645    1.85565  18.90249       0
e:40,FH-142  48.44619    2.21382  21.88352       0
e:15,FH-942 158.03982   13.25750  11.92079       0
e:20,FH-942  83.69820    5.80435  14.41991       0
e:25,FH-942  43.00616    3.49778  12.29527       0
e:30,FH-942  50.04317    3.60182  13.89387       0
e:35,FH-942  39.25428    3.76747  10.41925       0
e:40,FH-942  48.40065    5.74681   8.42217       0
于 2016-07-08T15:02:40.633 回答