Metóda najmenších štvorcov v nelineárnom modeli

Podobne ako pre lineárny model, aj teraz môžeme využiť funkcie dostupné v R, pomocu ktorých vieme vypočítať odhady parametrov. V tomto prípade bude ale minimalizácia súčtu štvorcov nelineárny optimalizačný problém, na výpočet ktorého môžeme použiť napríklad Gauss-Newtonovu metódu.

S použitím tejto metódy sú spojené dve otázky:

Fitovanie modelu

Uvažujme Michaelisov-Mentenov model \[ Y(x)=\frac{V_m x}{x+K}+\varepsilon_x \] Tento model popisuje napríklad rýchlosť chemickej reakcie v závisloti od koncentrácie subtrátu alebo účinok lieku v závislosti od jeho dávky.

Načítajme jednoduché dáta L.minor, ktoré popisujú chemickú reakciu. Parameter \(V_m\) vyjadruje horizontálnu asymptotu modelu a paramater \(K\) koncentráciu substrátu, pri ktorej je rýchlosť reakcie v polovici medzi 0 a \(V_m\).

library(nlstools)
## 
## 'nlstools' has been loaded.
## IMPORTANT NOTICE: Most nonlinear regression models and data set examples
## related to predictive microbiolgy have been moved to the package 'nlsMicrobio'
data(L.minor)
plot(rate ~ conc, data = L.minor, ylab = "Uptake rate (weight/h)", xlab = Substrate ~ concentration ~ (mmol ~ m^-3), pch=16)

Fitovanie týchto dát môžeme urobiť jednoducho pomocou príkazu nls:

L.minor.m1 <- nls(rate ~ Vm * conc/(K + conc), data = L.minor, start = list(K = 20, Vm = 120), trace = TRUE)
## 624.3282 :   20 120
## 244.546 :   15.92382 124.57148
## 234.5198 :   17.25299 126.43877
## 234.3595 :   17.04442 125.96181
## 234.3533 :   17.08574 126.04671
## 234.3531 :   17.07774 126.03017
## 234.3531 :   17.0793 126.0334
## 234.3531 :   17.07899 126.03275

Rozdiel oproti lm je v tom, že je v modeli nutné explicitne menovať parametre, pre ktoré potom v argumente start dodáme počiatočné hodnoty.

Podobne ako v lm, podrobný výstup dostaneme pomocou

summary(L.minor.m1)
## 
## Formula: rate ~ Vm * conc/(K + conc)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## K    17.079      2.953   5.784  0.00117 ** 
## Vm  126.033      7.173  17.570 2.18e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.25 on 6 degrees of freedom
## 
## Number of iterations to convergence: 7 
## Achieved convergence tolerance: 8.18e-06

Ďalšie hodnoty, ktoré by nás mohli zaujímať, sú reziduálny súčet štvorcov

deviance(L.minor.m1)
## [1] 234.3531

hodnota logaritmu vierohodnostnej funkcie

logLik(L.minor.m1)
## 'log Lik.' -24.86106 (df=3)

a fitované hodnoty

fitted(L.minor.m1)
## [1]  18.06066  28.56474  38.52679  71.09461  78.03806  87.78424  91.62683
## [8] 116.28685
## attr(,"label")
## [1] "Fitted values"

Vykreslenie fitovanej krivky

plot(rate ~ conc, data = L.minor, ylim = c(10, 130), ylab = "Uptake rate (weight/h)", xlab = Substrate ~ concentration ~ (mmol ~ m^-3))
concVal <- with(L.minor, seq(min(conc), max(conc), length.out = 100))
lines(concVal, predict(L.minor.m1, newdata = data.frame(conc = concVal)))
abline(h = coef(L.minor.m1)[2], lty = 2)

Kvalitu odhadu môžeme zobraziť aj vrstevnicovým grafom, do ktorého zakreslíme úrovne hodnôt \(RSS(K, V_m)\). Na to najskôr vytvoríme objekt obsahujúci informácie potrebné k zostrojeniu grafu pomocou nlsContourRSS.

L.minor.m1con <- nlsContourRSS(L.minor.m1)
## 0%100%
##  RSS contour surface array returned
plot(L.minor.m1con, nlev = 10)

Okrem vrstevníc tu vidíme aj 95% oblasť spoľahlivosti pre parametre (o tej bližšie na poslednom cvičení).

Hľadanie počiatočných hodnôt \(\theta\).

Okrem ‘uhádnutia’ vhodných počiatočných hodnôt parametra z prechádzajúcich experimentov alebo vizuálne z dát, môžeme pre niektoré modely použiť takzvané self-starter functions, ktoré pre konkrétny model tieto hodnoty vypočítajú z dát. Tieto funkcie sú zabudované v príkaze nls a ich zoznam dostaneme pomocou ?selfStart. Pre Michaelisov-Mentenov model máme funkciu SSmicmen:

L.minor.m2 <- nls(rate ~ SSmicmen(conc, Vm, K), data = L.minor)
summary(L.minor.m2)
## 
## Formula: rate ~ SSmicmen(conc, Vm, K)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## Vm  126.033      7.173  17.570 2.18e-06 ***
## K    17.079      2.953   5.784  0.00117 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.25 on 6 degrees of freedom
## 
## Number of iterations to convergence: 0 
## Achieved convergence tolerance: 2.496e-06

Lineárne parametre v nelineárnom modeli

Všimnime si, že funkciu stredných hodnôt v našom modeli môžeme písať ako \[ \frac{x}{x+K}\cdot V_m, \] parameter \(V_m\) teda v modeli vystupuje lineárne. Pre fixné hodnoty nelineárnych parametrov môžu byť lineárne parametre odhadnuté pomocou štandardnej lineárnej regresie. To znamená, že v skutočnosti potrebujeme menej počiatočných hodnôt.

Pre takúto situáciu máme v nls nastavenie algortihm=“plinear”:

L.minor.m3 <- nls(rate ~ conc/(K + conc), data = L.minor, algorithm = "plinear", start = list(K = 20))
summary(L.minor.m3)
## 
## Formula: rate ~ conc/(K + conc)
## 
## Parameters:
##      Estimate Std. Error t value Pr(>|t|)    
## K      17.079      2.953   5.784  0.00117 ** 
## .lin  126.033      7.173  17.570 2.18e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.25 on 6 degrees of freedom
## 
## Number of iterations to convergence: 8 
## Achieved convergence tolerance: 2.111e-06

Všimnime si, že lineárny parameter \(V_m\) nie je nls explicitne pomenovaný a vo výstupe je označený ako .lin.

Manuálny výpočet odhadov v nelineárnom modeli

Ukážeme si, ako získať odhady parametrov v nelineárnej regresii bez použitia funkcie nls. Dáta pochádzajú z príkladu 3.11 v učebnici.

norma<-function(x) sqrt(t(x)%*%x)

mm <- read.table("http://www.iam.fmph.uniba.sk/ospm/Filova/rm/data.txt", header = TRUE, sep = "")
plot(mm)

# funkcia odozvy
eta <- function(x,theta) theta[1]*x/(x+theta[2])

# vektor odoziev
eta_vec <- function(X,theta){
  N <- length(X)
  res <- rep(0,N)
  for(i in 1:N)
    res[i] <- eta(X[i],theta)
  res
}

# gradient funkcie odozvy
d_eta <- function(x,theta) c(x/(x+theta[2]),-theta[1]*x/(x+theta[2])^2)

# Jacobiho matica
J <- function(X,theta){
  m <- length(theta)
  N <- length(X)
  res <- matrix(rep(0,m*N),nrow=N,ncol=m)
  for(i in 1:N)
    res[i,] <- t(d_eta(X[i],theta))
  res
}

# Informacna matica
M <- function(X,theta) t(J(X,theta))%*%J(X,theta)

# Projektor
P <- function(X,theta) J(X,theta)%*%solve(M(X,theta),t(J(X,theta)))

# Startovacia hodnota
theta_0 <- c(1.13, 2.14)
eta_vec(mm$X, theta_0)
##  [1] 0.0000000 0.2140152 0.3598726 0.4656593 0.5458937 0.6088362 0.6595331
##  [8] 0.7012411 0.7361564 0.7658133 0.7913165 0.8134817 0.8329238 0.8501157
## [15] 0.8654267 0.8791494 0.8915187 0.9027256 0.9129264 0.9222509 0.9308072
J(mm$X, theta_0)
##            [,1]        [,2]
##  [1,] 0.0000000  0.00000000
##  [2,] 0.1893939 -0.08106635
##  [3,] 0.3184713 -0.11460911
##  [4,] 0.4120879 -0.12792839
##  [5,] 0.4830918 -0.13185839
##  [6,] 0.5387931 -0.13121470
##  [7,] 0.5836576 -0.12831383
##  [8,] 0.6205674 -0.12433353
##  [9,] 0.6514658 -0.11989517
## [10,] 0.6777108 -0.11533332
## [11,] 0.7002801 -0.11082865
## [12,] 0.7198953 -0.10647666
## [13,] 0.7371007 -0.10232480
## [14,] 0.7523148 -0.09839303
## [15,] 0.7658643 -0.09468563
## [16,] 0.7780083 -0.09119807
## [17,] 0.7889546 -0.08792098
## [18,] 0.7988722 -0.08484263
## [19,] 0.8078995 -0.08195030
## [20,] 0.8161512 -0.07923117
## [21,] 0.8237232 -0.07667276
M(mm$X, theta_0)
##           [,1]       [,2]
## [1,]  9.007998 -1.3256391
## [2,] -1.325639  0.2251313

Naprogramujeme (pomalú) verziu algoritmu, ktorý používa nls.

# Iteracny vypocet noveho theta
theta_nove <- function(X,Y,theta,mu,lambda){
  m <- length(theta)
  Mtheta <- M(X,theta)
  g <- t(J(X,theta)) %*% (Y-eta_vec(X,theta))
  res <- theta + lambda * solve(Mtheta + mu*diag(rep(1,m))) %*% g
  res
}

# Algoritmus vypoctu MNS odhadu
tolerancia <- 0.0001
iter <- 0
iter_max <- 1000000
krit_zastavenia <- 100000000
theta<-theta_0

while (iter <= iter_max && krit_zastavenia > tolerancia){
  iter <- iter+1
  theta <- theta_nove(mm$X,mm$Y,theta,0,lambda=1)
  #theta <- theta_nove(mm$X,mm$Y,theta,0,lambda=1/iter) # pri tejto dlzke kroku je algoritmus pomaly
  krit_zastavenia <- norma(2*t(J(mm$X,theta))%*%(mm$Y-eta_vec(mm$X,theta)))
  }
theta_hat<-theta
# rezíduá
nu<-mm$Y-eta_vec(mm$X,theta_hat)

# odhad sigma2
N <- length(mm$X)
m<-2
s2<-(norma(nu))^2/(N-m)

#Skutocne hodnoty theta_1=1.2, theta_2=0.8, sigma^2=0.005
theta_hat
##           [,1]
## [1,] 1.1848407
## [2,] 0.7297557
iter
## [1] 9
s2
##             [,1]
## [1,] 0.004181564
# vykreslenie
osX <- seq(from=0, to=10, by=0.01)
Yhat <- eta_vec(osX, theta_hat)
Ytrue <- eta_vec(osX, c(1.2,0.8))

plot(osX, Yhat, type="l", lwd=3)
points(osX, Ytrue, col="red", lwd=3, type="l", lty=2)
points(mm$X, mm$Y)