нелинейная регрессия в R с несколькими наборами данных

Я изучаю R и в настоящее время использую его для нелинейной регрессии (которую я также изучаю). У меня есть два набора данных (продолжительность операции на разных машинах), и я могу найти хорошую нелинейную регрессию для каждого из этих наборов. Теперь я хотел бы найти лучшую регрессию, которая минимизирует сумму обеих остаточных сумм квадратов.

Вот что у меня есть:

A <- c(1:5)
B <- c(100, 51, 32, 24, 19)
C <- c(150, 80, 58, 39, 29)

df <- data.frame (A,B,C)

f <- B ~ k1/A + k2
g <- C ~ k1/A + k2

n <- nls(f, data = df, start = list(k1=10, k2=10))
p <- nls(g, data = df, start = list(k1=10, k2=10))
n
#Nonlinear regression model
#  model: B ~ k1/A + k2
#   data: df
#     k1      k2 
#101.595  -1.195 
# residual sum-of-squares: 2.619

#Number of iterations to convergence: 1 
#Achieved convergence tolerance: 2.568e-07

p
#Nonlinear regression model
#  model: C ~ k1/A + k2
#   data: df
#     k1      k2 
#148.044   3.593 
# residual sum-of-squares: 54.19

#Number of iterations to convergence: 1 
#Achieved convergence tolerance: 1.803e-07

Константы k1 и k2 (конечно) разные для обоих наборов (B и C), мне интересно, как мне удалось найти конкретный k1 и конкретный k2, которые дают «лучшее» решение для обоих наборов данных.

Надеюсь, мое объяснение будет понятно. В противном случае то, что я пытаюсь найти, иногда (по крайней мере, здесь) называется глобальной нелинейной регрессией.

РЕДАКТИРОВАТЬ: я также хотел бы знать, как я могу сказать R, чтобы избежать отрицательных значений для определенного параметра. В этом случае я хотел бы, чтобы k2 был положительным.


person claf    schedule 04.09.2013    source источник


Ответы (1)


Если вам нужны идентичные параметры, вы должны просто объединить свои данные:

df2 <- data.frame(Y=c(df$B,df$C), X=rep(df$A, 2))
p <- nls(Y ~ k1/X + k2, 
         data = df2, 
         start = list(k1=10, k2=10), 
         lower = c(0, 0), 
         algorithm = "port")
summary(p)

#  Formula: Y ~ k1/X + k2
#  
#  Parameters:
#    Estimate Std. Error t value Pr(>|t|)    
#  k1  124.819     18.078   6.904 0.000124 ***
#    k2    1.199      9.781   0.123 0.905439    
#  ---
#    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#  
#  Residual standard error: 16.59 on 8 degrees of freedom
#  
#  Algorithm "port", convergence message: both X-convergence and relative convergence (5)

Изменить:

Если вы хотите, чтобы один параметр был одинаковым, а другой менялся, вы можете использовать модель смешанных эффектов. Однако я не знаю, как указать ограничения для этого (я считаю, что это непростая задача, но, возможно, ее можно решить путем перепараметризации).

library(nlme)

library(reshape2)
df3 <- melt(df, id.vars="A")

r <- nlme(value ~ k1/A + k2, 
          data = df3, 
          start = c(k1=10, k2=10), 
          fixed = k1 + k2 ~1,
          random = k2 ~ 1|variable)

summary(r)
#  Nonlinear mixed-effects model fit by maximum likelihood
#  Model: value ~ k1/A + k2 
#  Data: df3 
#  AIC      BIC    logLik
#  83.11052 84.32086 -37.55526
#  
#  Random effects:
#    Formula: k2 ~ 1 | variable
#                k2 Residual
#  StdDev: 12.49915 7.991013
#  
#  Fixed effects: k1 + k2 ~ 1 
#         Value Std.Error DF   t-value p-value
#  k1 124.81916  9.737738  7 12.818086  0.0000
#  k2   1.19925 11.198211  7  0.107093  0.9177
#  Correlation: 
#         k1    
#  k2 -0.397
#  
#  Standardized Within-Group Residuals:
#    Min         Q1        Med         Q3        Max 
#  -1.7520706 -0.5273469  0.2746039  0.5235343  1.4971808 
#  
#  Number of Observations: 10
#  Number of Groups: 2 

coef(r)
#          k1        k2
#  B 124.8192 -10.81835
#  C 124.8192  13.21684
person Roland    schedule 04.09.2013
comment
Итак, я предполагаю, что «lower=c(0,0)» — это минимальные значения для k1 и k2? Как насчет алгоритма порта? хорошо, я погуглю это :) - person claf; 04.09.2013
comment
Всегда ли следует выполнять нелинейную регрессию для нескольких наборов данных? Что делать, если я хочу, чтобы k1 был общим для обоих наборов данных, а k2 может быть конкретным? - person claf; 04.09.2013
comment
@claferri Да, lower указывает ограничения. Если вам нужны нижние или верхние границы, вы должны использовать алгоритм порта с nls. Читать help("nls"). - person Roland; 04.09.2013
comment
@ Роланд с nlme, можно ли указать, что k2 должен быть положительным? Не могли бы вы также кратко объяснить синтаксис «фиксированный» и «случайный»? - person claf; 04.09.2013
comment
Кажется, НАСТОЛЬКО подавляю @Roland в своих комментариях, пришлось поставить пробел между @ и Roland... странно! - person claf; 04.09.2013
comment
@claferri Это подавлено, потому что я все равно получаю уведомления о комментариях к моему ответу. Как я писал в ответе, я не знаю, как явно указать ограничения в nlme. Если вы не знакомы с моделями смешанных эффектов, здесь не место их объяснять. Синтаксис объясняется в документации к функции. - person Roland; 04.09.2013
comment
На самом деле, я пытаюсь изменить свою модель, чтобы добавить параметр k3, который должен быть специфичным для набора данных, например k2, но я не могу найти способ добавить его в «фиксированный» и «случайный», документация nlme недостаточно ясна. для меня, я продолжаю получать ошибки. - person claf; 04.09.2013
comment
Что-то вроде nlme(value ~ k1/A + k2 + k3*A, data = df3, start = c(k1=10, k2=10, k3=10), fixed = k1 + k2 + k3 ~1,random = k2 + k3 ~ 1|variable)? - person Roland; 04.09.2013
comment
Это то, что я пробовал, я думаю, ошибка означает что-то еще. Все равно спасибо за помощь! - person claf; 04.09.2013
comment
Трудно сказать без сообщения об ошибке. Это могут быть проблемы сходимости, переобучения или ряд других проблем. - person Roland; 04.09.2013
comment
Я пытаюсь nlme(value ~ k1/(k3*A) + k2, data = df3, start = c(k1=10, k2=10, k3=10), fixed = k1 + k2 + k3 ~1,random = k2 + k3 ~ 1|variable) и думаю, что это то, что я не могу делать с nlme. Ошибка заключается в том, что коэффициент уменьшения на каждом этапе снижается ниже минимального шага до PNLS (моя ошибка на французском языке, поэтому я не уверен в переводе PNLS). - person claf; 04.09.2013
comment
@claferri Эта проблема не имеет ничего общего с nlme. k1/(k3*A) + k2 можно упростить до k4/A + k2. Математически невозможно разделить k1 и k3. - person Roland; 04.09.2013
comment
Это то, о чем я догадывался, но я подумал, что, может быть, тот факт, что k1 был глобальным, а k3 - нет, может изменить проблему... Мой математический фон кажется слабым :). Надеюсь, теперь я понимаю, почему это ничего не меняет! :) - person claf; 05.09.2013