Ilustração do método Newton-Raphson em regressão não linear
Trajetória obtida com os valores consecutivos nas iterações do método Newton-Raphson para estimação de parâmetros em um modelo de regressão não linear.
Em regressão não linear, as estimativas dos parâmetros são obtidas por emprego de algum método numérico. Um dos métodos empregados é o de Newton-Raphson. Tal método é usado para encontrar raízes de uma função. Portanto, não é um método de maximação/mínimização de funções, mas usa o fato equivalente de que a derivada no ponto crítico de uma função é zero.
Documentação sobre o método Newton-Raphson existe em abundância. O que motivou fazer esse post foi não ter encontrado nenhum que exemplificasse aplicação para estimação de parâmetros em um modelo de regressão não linear. O estimador de mínimos quadrados para um modelo de regressão é aquele que torna mínima a função
O método Newton-Raphson vai procurar as raízes da derivada dessa função com relação à cada parâmetro. Então, se 
O método se baseia na aproximação de primeira ordem da função 
O valor de 
sendo que
é a matriz de derivadas que tem dimensão 
Esse procedimento numérico foi implementado no script abaixo para estimação dos parâmetros do modelo Michaelis-Menten e o modelo Exponencial reparametrizado que desenvolvi na minha Dissertação de Mestrado. No início eu monto as funções separadas e no final eu faço uma função que recebe os ingredientes mínimos e faz a estimação dos parâmetros. No final do script tem um procedimento gráfico para seleção dos chutes iniciais e acompanhar a trajetória durante as iterações. Até a próxima ridícula.
#-----------------------------------------------------------------------------
# método Newton-Raphson
# objetivo é mínimizar a soma de quadrados dos resíduos mas iremos
# encontrar as raízes da derivada dessa função pelo método de Newton-Raphson
#-----------------------------------------------------------------------------
# parâmetros, modelo e dados simulados deles
theta <- c(A=10, B=2) # parâmetros
x <- 1:12
f <- function(theta, x){
## função que retorna os valore preditos para theta e x fornecidos
## theta: vetor de parâmetros A e B do modelo Michaelis-Menten
## x: vetor de valores da variável independente
haty <- theta[1]*x/(theta[2]+x)
return(haty)
}
set.seed(222); e <- rnorm(x,0,0.1)
y <- f(theta=theta, x=x)+e
plot(y~x)
curve(f(theta, x), add=TRUE, col=2)
#-----------------------------------------------------------------------------
# função objetivo, derivada primeira e segunda com relação à theta
S0 <- function(theta, f, y, x){
## função que retorna a soma de quadrados dos resíduos para
## theta: vetor de parâmetros do modelo f
## f: função do modelo de regressão não linear que depende de theta e x
## y e x: vetores dos valores das variáveis dependente e independente
sqr <- sum((y-f(theta, x))^2)
return(sqr)
}
S0(theta=theta, f=f, y=y, x=x)
# derivada primeira, dimensão 1 x p
D(expression((y-A*x/(B+x))^2), "A") # derivada de S0 com relação a A
D(expression((y-A*x/(B+x))^2), "B") # derivada de S0 com relação a B
# derivada segunda, dimensão p x p
D(D(expression((y-A*x/(B+x))^2), "A"), "A") # derivada de S0 com relação a A e A
D(D(expression((y-A*x/(B+x))^2), "B"), "B") # derivada de S0 com relação a B e B
D(D(expression((y-A*x/(B+x))^2), "A"), "B") # derivada de S0 com relação a A e B
D(D(expression((y-A*x/(B+x))^2), "B"), "A") # derivada de S0 com relação a B e A
#-----------------------------------------------------------------------------
# para ter essas funções vamos usar a deriv3()
ff <- deriv3(~(y-A*x/(B+x))^2,
c("A","B"),
function(x, y, A, B){ NULL })
ff(x=x, y=y, theta[1], theta[2])
#-----------------------------------------------------------------------------
# mas as matrizes são 1 x p e p x p, então falta fazer a soma em n
fff <- function(x, y, A, B){
## função que retorna gradiente e hessiano avaliado em theta
## y e x: vetores dos valores das variáveis dependente e independente
## A e B: valores dos parâmetros do modelo Michaelis-Menten
m <- ff(x, y, A, B)
g <- attr(m, "gradient"); G <- apply(g, 2, sum)
h <- attr(m, "hessian"); h <- matrix(h, ncol=2*2)
H <- apply(h, 2, sum); H <- matrix(H, 2, 2)
return(list(G=G, H=H))
}
fff(x=x, y=y, theta[1], theta[2])
#-----------------------------------------------------------------------------
# agora é só fazer o Newton-Raphson
th0 <- c(9,2); tol <- 1e-10; niter <- 100;
dis <- 1; i <- 0
while(dis>tol & i<=niter){
GH <- fff(x=x, y=y, th0[1], th0[2])
G <- GH[[1]]
iH <- solve(GH[[2]])
th1 <- th0-iH%*%G
dis <- sum((th0-th1)^2)
i <- i+1
print(c(i, dis, th1))
th0 <- th1
}
n0 <- nls(y~A*x/(B+x), start=c(A=9, B=2))
coef(n0)
iH/vcov(n0) # hessiano é proporcional à matriz de covariância dos parâmetros
#-----------------------------------------------------------------------------
# vamos incluir todos os elementos em uma única função!
newton <- function(formula, data, start, tol=1e-10, niter=100, trace=TRUE){
## formula: expressão da resposta como função das covariáveis e parâmetros
## data: data.frame que contém dos dados
## start: vetor ou lista nomeados de valores iniciais para os parâmetros
## tol: norma mínima entre valores consecutivos dos parâmetros
## niter: número máximo de iterações do algortimo
## trace: se TRUE imprime o histórico de iterações
form <- paste("~(",paste(formula[2:3], collapse="-"),")^2", sep="")
vars <- all.vars(formula)
nvars <- length(vars)
pars <- 1:(nvars+1)
l <- as.list(pars)
names(l)[1:nvars] <- vars
f <- as.function(l)
body(f) <- NULL
ff <- deriv3(as.formula(form), names(start), f)
th0 <- unlist(start)
i <- 0; dis <- 100; iter <<- matrix(th0, ncol=length(start))
while(dis>tol & i<=niter){
m <- do.call(ff, c(as.list(data), as.list(th0)))
g <- attr(m, "gradient"); G <- apply(g, 2, sum)
h <- attr(m, "hessian"); h <- matrix(h, ncol=length(start)^2)
H <- apply(h, 2, sum); H <- matrix(H, nrow=length(start))
iH <- solve(H)
th1 <- th0-iH%*%G
dis <- sum((th0-th1)^2)
i <- i+1
iter <<- rbind(iter, t(th1))
if(trace==TRUE) print(c(i, dis, th1))
th0 <- th1
}
th0 <- as.vector(th0)
names(th0) <- names(start)
## retorna o último valor iterado e a uma matriz com o histórico
return(list(theta=th0, iter=iter))
}
#-----------------------------------------------------------------------------
# vamos estudar vetor de chutes presentes nessa lista
start.l <- list(c(A=10,B=2.5), c(A=10,B=0.1), c(A=14,B=2.5),
c(A=15.13, B=1.7), c(A=18.4,B=0.13))
data <- data.frame(x=x, y=y)
formula <- y~A*x/(B+x)
start <- start.l[[4]]
nr <- newton(formula, data, start)
nr$theta
sqe <- function(A, B, y, x){ hy <- (A*x)/(B+x); sum((y-hy)^2) }
SQE <- Vectorize(sqe, c("A", "B"))
A.grid <- seq(-3,15,l=100)
B.grid <- seq(-3,15,l=100)
sqe.surf <- outer(A.grid, B.grid, SQE, data$y, data$x)
#png("f023.png", w=500, h=500)
contour(A.grid, B.grid, sqe.surf, levels=(1:45)^2, asp=1,
xlab=expression(theta[0]), ylab=expression(theta[1]), col="gray70")
for(i in seq(nrow(iter)-1)){
arrows(iter[i,1], iter[i,2],
iter[i+1,1], iter[i+1,2],
col=2, length=0.1)
Sys.sleep(0.1)
}
#dev.off()
#-----------------------------------------------------------------------------
# com o mouse selecione regiões da superfície para serem valores iniciais
loc <- locator()
loc <- as.data.frame(loc)
names(loc) <- c("A","B")
start.l <- split(loc, f=1:nrow(loc))
#-----------------------------------------------------------------------------
# ver como fica com o modelo Exponencial
n0 <- nls(y~A*(1-exp(-log(2)*x/B)), start=c(A=9, B=2))
coef(n0)
data <- data.frame(x=x, y=y)
formula <- y~A*(1-exp(-log(2)*x/B))
start <- list(A=5, B=0.5)
nr <- newton(formula, data, start)
nr$theta
sqe <- function(A, B, y, x){ hy <- A*(1-exp(-log(2)*x/B)); sum((y-hy)^2) }
SQE <- Vectorize(sqe, c("A", "B"))
A.grid <- seq(-10,20,l=100)
B.grid <- seq(-10,20,l=100)
sqe.surf <- outer(A.grid, B.grid, SQE, data$y, data$x)
contour(A.grid, B.grid, sqe.surf, levels=(1:45)^2, asp=1,
xlab=expression(theta[0]), ylab=expression(theta[1]), col="gray70")
for(i in seq(nrow(iter)-1)){
arrows(iter[i,1], iter[i,2],
iter[i+1,1], iter[i+1,2],
col=2, length=0.1)
Sys.sleep(0.1)
}
#-----------------------------------------------------------------------------
