Regressão não linear
Modelos de regressão linear e não linear
Modelos de regressão linear
Até o presente momento do curso, consideramos modelos lineares nos parâmetros.
Por exemplo:
1) Modelo linear geral:
Yi  0  1 X i1  ...   p1 X i , p1  i
1) Modelo polinomial:
Yi  0  1 X i1  2 X i21  i
1
1) Modelo com variáveis transformadas:
log10 Yi  0  1 X i1  2 exp(X i 2 )   i
Os modelos lineares, podem ser escritos, na forma:
Yi  f (Xi , )  i
Onde Xi é o vetor de observações das variáveis preditoras para o i-ésimo
caso:
 1 
 X 
 i1 
Xi   . 


.


 X i , p 1 
 é o vetor dos parâmetros, e f(Xi,) representa o valor esperado E(Yi), o qual
para o modelo linear é:
'
f (Xi , β)  Xiβ
2
Nos modelos lineares, o problema de estimação dos parâmetros, cai no problema
de resolver um sistema de equações lineares com relação aos coeficientes de
regressão desconhecidos. Existe uma solução única e, portanto, obtemos uma
forma analítica de estimação dos parâmetros. Esta forma é a mesma para
qualquer modelo e qualquer conjunto de dados.
Além disso, como os coeficientes são combinações lineares das observações,
pela teoria estatística, demonstra-se que a distribuição amostral dos coeficientes
de regressão segue uma distribuição t, assim, podemos realizar os testes de
hipóteses, calcular os intervalos de confiança para esses coeficientes.
Modelos de regressão não linear
Existe, entretanto, muitas situações nas quais não é desejável, ou mesmo
possível, descrever um fenômeno através de um modelo de regressão linear.
Ao invés de se fazer uma descrição puramente empírica do fenômeno em
estudo, pode-se, a partir de suposições importantes sobre o problema
(freqüentemente dadas através de uma ou mais equações diferenciais),
trabalhar no sentido de obter uma relação teórica entre as variáveis
observáveis de interesse. O problema, diferentemente do caso linear, é que os
parâmetros entram na equação de forma não linear, assim, nós não podemos
simplesmente aplicar fórmulas para estimar os parâmetros do modelo.
3
Outra vantagem dos modelos não lineares é obter parâmetros que são
facilmente interpretáveis.
Em muitas situações, necessita-se menos parâmetros nos modelos não lineares
do que nos lineares, isto simplifica e facilita a interpretação.
Os modelos não lineares podem ser escritos como:
Yi  f ( Xi ,γ )   i
f(Xi, ) é uma função não linear; os erros, i, tem média zero, variância
constante, e não são correlacionados. Assume-se que os erros apresentam
distribuição normal, são independentes e com variância constante.  é o vetor
de parâmetros do modelo.
Dois exemplos de modelos não lineares.
1) Modelo exponencial
Yi   0 exp( 1 X i )  i (1)
0 e 1são os parâmetros do modelo; Xi são constantes conhecidas (variável
preditora) e i são os termos do erro, independentes, com distribuição normal
de média 0 (zero) e variância 2.
4
Diferenciando f com respeito a 0 e 1 obtemos (usando MAPPLE):
f
 0
 exp( 1X)
f
 1
  0 Xexp( 1X)
Como estas derivadas envolvem pelo menos um dos parâmetros, o modelo é
reconhecido como não linear.
Um modelo exponencial mais geral:
Yi   0   1 exp( 2 X i )  i
(2)
Veja figura.
5
S
c
a
e
t
rp
lo
t
y:=
1
0
0
-5
0
*
e
xp
(-2
*
x)
1
1
0
1
0
0
E(X)
9
0
8
0
7
0
6
0
5
0
0
,0
0
,5
1
,0
1
,5
2
,0
2
,5
3
,0
3
,5
X
Estes modelos exponenciais são muito utilizados em estudos de crescimento,
onde a taxa de crescimento num dado tempo X é proporcional a quantidade de
crescimento restante (final) que ocorre com o aumento do tempo, e 0 representa
o crescimento máximo
6
2) Modelo logístico
0
Yi  1 1 exp( 2 X i )   i
(3)
i são os termos do erro, independentes, com distribuição normal de média 0
(zero) e variância 2. A função esperada é:
0
f ( X, γ )  1 1 exp(
 2 Xi )
y:=
1
0
/(1
+
2
0
*
e
xp
(-2
*
x))
1
2
1
0
8
E(Y)
6
4
2
0
-2
-0
,5 0
,0
0
,5 1
,0
1
,5
X
2
,0 2
,5
3
,0
3
,5
O modelo logístico é
muito usado para
variáveis qualitativas.
Exemplo: acertos na
cache (acerta/não
acerta). Neste caso, os
erros não tem mais
distribuição normal
com variância
constante.
7
Alguns aspectos do uso de modelos não lineares:
• os modelos não lineares tem uma base teórica, os parâmetros dos modelos
fornecem um maior conhecimento sobre o fenômeno em estudo do que os
modelos lineares.
• os modelos não lineares, geralmente fornecem um bom ajuste, com menos
parâmetros do que os modelos lineares.
• A transformação de um modelo não linear em um modelo linear nos
parâmetros, se por um lado facilita o processo de ajuste, implica em fazer
suposições não realísticas sobre o termo dos erros (distribuição normal com
variância constante); além disso, perde-se informação sobre os erros padrões
dos parâmetros originais.
• Além disso, existem modelos que são intrinsicamente não lineares, isto é, não
podem ser linearizados por transformação.
• Embora vamos usar variáveis contínuas como variáveis independentes, não há
razão para que as variáveis independentes, nos modelos não lineares, sejam
contínuas. Ao contrário, podemos fazer uso de variáveis dummy para indicar a
presença ou ausência de um grupo, ou codificar diferenças entre indivíduos
(dados de medidas repetidas).
8
• Estimação de modelos não lineares, é um bom exemplo de que a despeito de
se obter os resultados no computador, não significa que os resultados sejam
corretos ou razoáveis.
A forma geral do modelo não linear
Yi  f ( X i , γ )   i
 X i1 
X 
 i2 
Xi   . 
( q x 1)
 
 . 
 X iq 
 
(4)
 0 
 
 1 
γ  . 
(p x 1)


 . 
 p 1 


Onde f(Xi, ) é a função esperada para o i-ésimo caso.
9
Estimação dos parâmetros
Métodos: »Mínimos quadrados
»Máxima verossimilhança
Importante: nos modelos não lineares não é possível encontrarmos formas
analíticas para os estimadores de mínimos quadrados ou máxima
verossimilhança. Ao invés, métodos numéricos devem ser usados
juntamente com os métodos referidos e, isto, requer cálculos
computacionais intensivos. Sempre usamos softwares computacionais.
Exemplo
Um administrador de um hospital deseja ajustar um modelo de regressão para
estimar o tempo de recuperação depois que o paciente saiu do hospital devido
a uma doença grave. A variável preditora é o número de dias que o paciente
ficou hospitalizado (X), e a variável resposta é um índice de prognóstico para
o tempo de recuperação (Y), onde, valores grandes indicam um bom
prognóstico. A seguir temos os dados e o diagrama de dispersão:
10
Dados para pacientes com doença grave.
Pacientes
Dias hospitalizados
i
Xi
1
2
2
5
3
7
4
10
5
14
6
19
7
26
8
31
9
34
10
38
11
45
12
52
13
53
14
60
15
65
Prognóstico (índice)
Yi
54
50
45
37
35
25
20
16
18
13
8
11
8
4
6
11
S
c
a
tte
rp
lo
t
6
0
5
0
4
0
Prognóstico(índice)
3
0
2
0
1
0
0
-1
0
0
1
0
2
0
3
0
4
0
5
0
6
0
7
0
D
ia
sh
o
s
p
ita
liz
a
d
o
Encontrou-se na literatura que a relação entre a variável preditora e a variável
resposta segue o modelo:
Yi   0 exp( 1 X i )   i
Onde os i são os termos dos erros, independentes, com distribuição normal de
média 0 (zero) e variância 2 (constante). Precisamos estimar os parâmetros 0
e 1.
12
Método de mínimos quadrados na
regressão não linear
Como no modelo de regressão linear geral, o critério de mínimos
n
quadrados é:
Q   (Yi  f ( Xi , γ))2
(5)
i 1
O critério Q deve ser minimizado com respeito aos parâmetros de regressão
não linear 0, 1,..., p-1 para obter as estimativas de mínimos quadrados.
Métodos: 1) procura numérica e 2) equações normais de mínimos quadrados.
A diferença com a regressão linear é que a solução das equações normais
usualmente requer um método numérico iterativo, pois a solução analítica
geralmente não pode ser encontrada.
13
Exemplo: para os dados de pacientes com doença grave, a função
esperada é:
f (X,γ)   0 exp( 1 X )
O critério Q é dado por:
n
Q   (Yi   0 exp( 1 X i ))
2
i 1
Método da máxima verossimilhança:
Vamos considerar que os erros i são independentes, normalmente
distribuídos com variância constante. A função de máxima verossimilhança
é dada por:
n
1

2


L( γ, 2 ) 
exp

Y


exp(

X
)
0
1 i
 2 2  i

(2 2 )n / 2
i 1


1
Maximizar esta função com relação aos parâmetros, é idêntico a minimizar o
somatório na parte do expoente, portanto, chega-se aos mesmos estimadores com
os dois métodos.
14
Solução das equações normais
Para obter as equações normais para um modelo não linear
Yi  f (Xi , γ)  i
Precisamos minimizar o critério Q
n
Q   (Yi  f ( Xi , γ))2
i 1
com respeito aos parâmetros 0, 1,..., p-1-. As derivadas parciais de Q
com respeito aos k é:
n
 f ( Xi , γ ) 
Q
   2(Yi  f ( Xi , γ ))

 k i 1
  k 
15
Igualando-se as derivadas parciais a zero e, substituindo-se k por gk
(estimativas de mínimos quadrados), obtemos o sistema de equações normais
(p equações, k=0,1,...,p-1):
n
 f ( Xi , γ) 
 f ( Xi , γ) 
Yi 
  f ( Xi , g) 
0



i 1
  k  γg i 1
  k  γ g
n
(6)
Onde g é o vetor das estimativas de mínimos quadrados gk:
 g0 
g 
 1
g  . 
( p x 1)
 
 . 
g p -1 
 
As equações normais (6) são não lineares nas estimativas dos parâmetros gk,
portanto, difíceis de serem resolvidas. Dessa forma, vamos precisar de métodos
numéricos para obter uma solução das equações normais iterativamente.
16
Exemplo: para os dados de pacientes com doença grave, a função
esperada para o i-ésimo caso é:
f (Xi ,γ)   0 exp( 1 X i )
As derivadas parciais já foram mostradas anteriormente. Substituindo-se 0 e
1 pelas estimativas de mínimos quadrados g0 e g1, as equações normais (6)
são dadas por:
Y exp(g X )   g exp(g X ) exp(g X )  0
Y g X exp(g X )   g exp(g X ) g X exp(g X )  0
i
i
0
i
1
1
i
0
i
0
1
1
i
i
0
i
1
i
1
i
Procedendo-se a algumas simplificações, obtemos:
Y exp(g X )  g  exp(2 g X )  0
Y X exp(g X )  g  X exp(2 g X )  0
i
i
i
1
1
i
i
0
0
i
1
i
1
i
São equações não lineares nas estimativas dos parâmetros, assim, métodos
numéricos devem ser empregados(métodos iterativos).
17
Método de Gauss-Newton (Procura numérica
direta – Direct numerical search)
Na maioria dos problemas com modelos não lineares, é mais prático encontrar
as estimativas de mínimos quadrados por procedimentos de procura numérica
direta do que, inicialmente, obter as equações normais e, então, usar métodos
numéricos para encontrar a solução dessas equações iterativamente.
O método de Gauss-Newton, também conhecido como método da
linearização, usa uma expansão em série de Taylor para aproximar o modelo
de regressão não linear com termos lineares e, então, aplica mínimos
quadrados ordinário para estimar os parâmetros. Iterações desses passos
geralmente conduzem a uma solução para o problema de regressão não linear.
O método de Gauss-Newton inicia dando-se valores iniciais aos parâmetros 0,
1,..., p-1, denotados por:
g0( 0) , g1( 0) ,...,g (p0)1
Esses valores iniciais podem ser obtidos de estudos anteriores, conhecimentos
teóricos ou por uma grade de valores que minimize (5).
18
Com os valores iniciais dos parâmetros, aproximamos a função esperada f(Xi, )
para os n casos por termos lineares da expansão em série de Taylor, de primeira
ordem, em torno dos valores iniciais gk(0). Obtemos para o i-ésimo caso:
 f ( X i , γ ) 
(0)
f (Xi , γ)  f (Xi , g )   
(


g
k
k )

 k  γ  g ( 0 )
k 0 
p 1
(0)
(7)
Aqui g(0) é o vetor dos valores iniciais dos parâmetros. Observe que as
derivadas, assim como a f, são avaliadas em k=gk(0).
Fazendo-se:
f i 0  f ( Xi , g (0) )
 k( 0)  ( k  g k( 0) )
(7.A)
 f ( Xi , γ ) 
Dik( 0)  



k

 γ g ( 0 )
19
Podemos reescrever a aproximação (7) como:
p 1
f ( Xi , γ)  f i ( 0)   Dik( 0)  k( 0)
(8)
k 0
E uma aproximação para o modelo (4)
Yi  f (Xi , γ)  i
é dada por:
p 1
Yi  fi ( 0)   Dik( 0)  k( 0)   i
(9)
k 0
Passando fi(0) para o lado esquerdo e, denotando a diferença Yi- fi(0) por Yi(0),
p 1
temos:
Yi ( 0)   Dik( 0)  k( 0)   i
i  1,2,...,n
(10)
k 0
Observe que chegamos a uma aproximação para um modelo de regressão
linear.
20
Cada coeficiente de regressão k(0) representa a diferença entre os
verdadeiros parâmetros da regressão e as estimativas iniciais dos
mesmos. Assim, os coeficientes de regressão representam uma correção
que deve ser feita nos coeficientes de regressão iniciais. O propósito de
ajustar o modelo de regressão linear (10) é estimar os coeficientes de
regressão k(0) e usar essas estimativas para corrigir as estimativas iniciais
dos parâmetros de regressão.
O modelo (10) na forma matricial fica:
Y( 0)  D( 0)β( 0)  ε
Y( 0)
nx1
Y1  f1( 0 ) 


.



.



.

( 0) 
Yn  f n 
D( 0 )
nxp
(11)
 D10( 0 ) ... D1(,0p)1 
 .



 .

 .

 ( 0)

(0)
 Dn 0 ... Dn , p 1 
21
  0( 0) 
 . 


(0)
β  . 
( p x 1)
 . 
 ( 0) 
  p 1 
 1 
.
 
ε  . 
( n x 1)
.
 
 n
Observe as similaridades entre o modelo de regressão linear :
Y  Xβ  ε
A matriz D faz o papel da matriz X:
DX
Podemos, portanto, estimar os parâmetros (0) pelo método de mínimos
quadrados ordinários:
b(0)  ( D(0)' D(0) )1 D(0)'Y (0)
Usar um programa de computador que faça regressão múltipla, porém não
esquecer de especificar que não desejamos o intercepto.
22
Nós, então, usamos estas estimativas de mínimos quadrados para obter os
coeficientes de regressão estimados corrigidos gk(1) por meio de (7.A):
gk(1)  gk(0)  bk(0)
Onde gk(1) representa a estimativa corrigida de  k no fim da primeira iteração.
Na forma matricial, temos:
g(1)  g(0)  b(0)
(11.A)
Neste ponto, nós podemos verificar se os coeficientes de regressão corrigidos
representam uma melhoria na direção apropriada. Denotaremos o critério Q,
calculado nos coeficientes de regressão iniciais g(0), por SQE(0), ou seja,
SQE
( 0)
n
n
  (Yi  f ( Xi , g ))   (Yi  fi (0) )2
i 1
(0)
2
i 1
23
No final da primeira iteração, os coeficientes de regressão corrigidos são g(1).
Denotaremos o critério Q, calculado nos coeficientes de regressão g(1), por
SQE(1), ou seja,
SQE
(1)
n
n
  (Yi  f ( Xi , g ))   (Yi  fi (1) )2
i 1
(1)
2
i 1
Se o algoritmo de Gauss-Newton está na direção correta, SQE(1) deverá ser menor
do que SQE(0), pois os coeficientes de regressão no passo (1) deverão ser melhores.
O método de Gauss-Newton repete o procedimento como foi descrito, com g(1)
sendo, agora, usado como valores iniciais. Isto resulta num novo conjunto de
estimativas corrigidas, representadas por g(2), e teremos um novo critério SQE(2).
O processo iterativo continua até que as diferenças entre sucessivas estimativas
dos coeficientes g(s+1)-g(s) e/ou a diferença entre sucessivas soma de quadrados de
erros SQE(s-1)-SQE(s) tornam-se desprezíveis. As estimativas finais dos
coeficientes de regressão são representadas por g e a soma de quadrado dos erros
por SQE.
24
Exemplo: para os dados de pacientes com doença grave, a função é:
Yi   0 exp( 1 X i )   i
Usando o PROC NLIN do SAS, vamos fazer a análise estatística dos dados. O
programa é:
data doenca;
input obs dias
datalines;
1.000
2.000
2.000
5.000
3.000
7.000
4.000
10.000
5.000
14.000
6.000
19.000
7.000
26.000
8.000
31.000
9.000
34.000
10.000 38.000
11.000 45.000
12.000 52.000
13.000 53.000
14.000 60.000
15.000 65.000
;
indice;
54.000
50.000
45.000
37.000
35.000
25.000
20.000
16.000
18.000
13.000
8.000
11.000
8.000
4.000
6.000
proc print data=doenca; run;
proc nlin data=doenca method=gauss maxiter=20;
parms a=56.6646
b=-0.03797;
model indice = a*exp(b*dias);
der.a=exp(b*dias);
der.b=a*dias*exp(b*dias);
output out=doencaou p=predito r=residuo;
run;
Os valores iniciais de a e b, foram obtidos através de
uma regressão linear simples do modelo:
ln Y  ln γ0  γ1 X
25
Output do SAS:
Non-Linear Least Squares Iterative Phase
Method: Gauss-Newton
Iter
A
B
Sum of Squares
0
56.664600
-0.037970
56.086713
1
58.557844
-0.039533
49.463830
2
58.605484
-0.039585
49.459304
3
58.606531
-0.039586
49.459300
4
58.606565
-0.039586
49.459300
NOTE: Convergence criterion met.
Non-Linear Least Squares Summary Statistics
Source
DF Sum of Squares
Mean Square
Regression
Residual
Uncorrected Total
2
13
15
12060.540700
49.459300
12110.000000
6030.270350
3.804562
(Corrected Total)
14
3943.333333
Parameter
A
B
Estimate
Asymptotic
Std. Error
58.60656517
-0.03958645
1.4721603058
0.0017112939
Asymptotic 95 %
Confidence Interval
Lower
Upper
55.426158088 61.786972243
-0.043283475 -0.035889427
26
S
c
a
e
t
rp
lo
t
y
:=
5
8
,6
0
6
5
*e
x
p
(-0
,0
3
9
5
9
*x
)
1
1
0
9
0
Índice
7
0
5
0
3
0
1
0
-1
0
-1
0
0
1
0
2
0
3
0
4
0
5
0
6
0
7
0
D
ia
s
SQErro
49 , 4593
r 2  1  SQTotal

1

 0,9875  98,78%
Corrigdo
3943 ,333
27
Exercício: vamos considerar os dados de pacientes com doença grave.
Aplicar a transformação logarítmica e obter as estimativas iniciais dos
coeficientes de regressão.
A função resposta é:
Yi   0 exp( 1 X i )   i
Aplicando o logaritmo, obtemos:
log Yi  log  0   1 X i
Podemos aproximar o modelo exponencial pelo modelo linear:
Yi '   0  1 X i   i
onde :
Yi '  log Yi
 0  log  0
1   1
28
Com o uso do PROC IML do SAS obtemos:
proc iml;
reset print;
Y={54, 50, 45, 37, 35, 25, 20, 16, 18, 13, 8, 11, 8, 4, 6};
X={1 2, 1 5, 1 7, 1 10, 1 14, 1 19, 1 26, 1 31, 1 34, 1 38, 1 45, 1 52,
1 53, 1 60, 1 65};
YT=log(Y);
XLX=X`*X;
XLXinv=inv(xlx);
b=XLXinv*x`*yt;
b0=4,0371
b1=-0,03797
 g 0( 0)  exp( b0 )  56,6646
 g1( 0)  b1  0,03797
29
A soma de quadrados do erro no passo zero, SQE(0), requer o cálculo da
função de regressão não linear
f ( X,γ )   0 exp( 1 X i )
(12)
para cada caso, utilizando os valores iniciais. Por exemplo, para o primeiro
caso, onde X1=2, obtemos:
f ( X1, g (0) )  f1(0)  g0(0) exp( g1(0) X1 )
 56,6646 * exp(-0,03797(2))  52,5208
Para os 15 casos, temos:
f(0) =
52.520821
46.866338
43.439088
38.76236
33.300409
27.542208
21.11386
17.462918
15.58283
13.387075
10.262533
7.8672587
7.574139
5.8063357
4.8023226
/* valores iniciais */
g00=56.6646; g10=-0.03797;
X2=X[1:15,2];
/* funcao de regressão
*/
f=g00*exp(g10*X2);
30
Para o primeiro caso, Y1=54, portanto, o desvio da resposta esperada é:
Y1(0)  Y1  f1(0)  54  52,5208  1,4792
Y(0)
=
1.4791792
3.133662
1.5609122
-1.76236
1.6995911
-2.542208
-1.11386
-1.462918
2.4171698
-0.387075
-2.262533
3.1327413
0.425861
-1.806336
1.1976774
Y0=Y-f;
/* soma de quadrados do erro no
passo zero */
SQE0=Y0`*Y0;
A soma de quadrados do erro no passo zero, SQE(0), vale:
SQE ( 0)   (Yi  fi ( 0) )   (Yi ( 0) ) 2
 1,47952  ...  1,19772  56,0869
31
Para obter as estimativas dos coeficientes corrigidos, precisamos calcular
D(0). Para obter esta matriz, precisamos das derivadas parciais da função de
regressão (12) calculadas em  = g(0).
Para ilustrar, vamos tomar o caso 1, para o qual X1=2. Assim, o valor das
derivadas parciais em g(0) são:
( 0)
D10
 exp( g1( 0) X 1 )  exp( 0,03797(2))  0,92687
( 0)
D11
 g 0( 0)  X 1 exp( g1( 0) X 1 )  56,6646(2) exp( 0,03797(2))  105,0416
D(0) =
0.9268718
0.8270832
0.7666001
0.6840666
0.5876757
0.4860567
0.3726111
0.3081804
0.2750011
0.2362511
0.1811101
0.138839
0.1336662
0.1024685
0.08475
105.04164
234.33169
304.07361
387.6236
466.20573
523.30196
548.96035
541.35047
529.81623
508.70884
461.81398
409.09745
401.42937
348.38014
312.15097
/*derivadas parciais calculadas em g(0)*/
D0_0=exp(g10*X2);
D1_0=g00*X2#exp(g10*X2);
D0=D0_0||d1_0;
32
Agora, podemos obter as estimativas de mínimos quadrados b(0), fazendo a
regressão de Y(0) sobre as 2 variáveis X na matriz D(0). Continuando com o
nosso programa no IML do SAS obtemos:
b(0) =
1.893244
-0.001563
b0=inv(D0`*D0)*D0`*Y0;
Usando 11.A, obtemos os coeficientes de regressão corrigidos g(1):
g (1)  g (0)  b (0)
 56,6646   1,8932 

   0,001563
0,03797

 

 58,5578 


- 0,03953
/* novas estimativas corrigidas
*/
g0=g00//g10;
g1=g0+b0;
Aqui, chegamos ao final da primeira iteração com:
g0(1)  58,5578 g1(1)  0,03953
A soma de quadrados residual na primeira iteração vale:
33
SQE
(1)
n
  (Yi  fi )
i 1
= 49.46383
(1) 2
f1=g1[1,1]*exp(g1[2,1]*X2);
Y1=Y-f1;
/* soma de quadrados do erro na iteracao 1 */
SQE1=Y1`*Y1;
Observe que houve uma redução nas somas de quadrados dos resíduos.
Continuação do exercício: Faça as próximas três iterações, verifique se foi
encontrado o critério de convergência ((SQE(s)-SQE(s-1)) <0,0001) e escreva o
modelo.
34
proc iml;
reset print;
Y={54, 50, 45, 37, 35, 25, 20, 16, 18, 13, 8, 11, 8, 4, 6};
X={1 2, 1 5, 1 7, 1 10, 1 14, 1 19, 1 26, 1 31, 1 34, 1 38, 1 45, 1 52, 1 53, 1 60, 1 65};
YT=log(Y);
XLX=X`*X;
XLXinv=inv(xlx);
b=XLXinv*x`*yt;
/* valores iniciais */
g00=56.6646; g10=-0.03797;
X2=X[1:15,2];
f=g00*exp(g10*X2);
Y0=Y-f;
/* soma de quadrados do erro no passo zero */
SQE0=Y0`*Y0;
/* derivadas parciais calculadas em g(0)
*/
D0_0=exp(g10*X2);
D1_0=g00*X2#exp(g10*X2);
D0=D0_0||d1_0;
b0=inv(D0`*D0)*D0`*Y0;
/* novas estimativas corrigidas - iteracao 1 */
g0=g00//g10;
g1=g0+b0;
f1=g1[1,1]*exp(g1[2,1]*X2);
/* residuos da iteracao 1 */
Y1=Y-f1;
/* soma de quadrados do erro na iteracao 1 */
SQE1=Y1`*Y1;
/*********************fim da iteracao 1 ****************/
35
/* derivadas parciais calculadas em g(1)
*/
D0_1=exp(g1[2,1]*X2);
D1_1=g1[1,1]*X2#exp(g1[2,1]*X2);
D1=D0_1||d1_1;
/* estimativas corrigidas na iteracao 2 */
b1=inv(D1`*D1)*D1`*Y1;
/* novas estimativas corrigidas - iteracao 2 */
g2=g1+b1;
f2=g2[1,1]*exp(g2[2,1]*X2);
/* residuos da iteracao 2
*/
Y2=Y-f2;
/* soma de quadrados do erro na iteracao 2 */
SQE2=Y2`*Y2;
/***********fim da iteracao 2 *******************/
36
/* derivadas parciais calculadas em g(2)
D0_2=exp(g2[2,1]*X2);
D1_2=g2[1,1]*X2#exp(g2[2,1]*X2);
D2=D0_2||d1_2;
/* estimativas corrigidas na iteracao 3 */
b2=inv(D2`*D2)*D2`*Y2;
g3=g2+b2;
f3=g3[1,1]*exp(g3[2,1]*X2);
/* residuos da iteracao 3 */
Y3=Y-f3;
/* soma de quadrados do erro na iteracao 3
SQE3=Y3`*Y3;
/************fim da iteracao 3 */
/* derivadas parciais calculadas em g(3)
D0_3=exp(g3[2,1]*X2);
D1_3=g3[1,1]*X2#exp(g3[2,1]*X2);
D3=D0_3||d1_3;
/* estimativas corrigidas na iteracao 4 */
b3=inv(D3`*D3)*D3`*Y3;
g4=g3+b3;
f4=g4[1,1]*exp(g4[2,1]*X2);
/* residuos da iteracao 4 */
Y4=Y-f4;
/* soma de quadrados do erro na iteracao 4
SQE4=Y4`*Y4;
/************fim da iteracao 4 */
*/
*/
*/
*/
37
Comentários:
1) A escolha das estimativas iniciais no método de Gauss-Newton é muito
importante, pois uma má escolha pode resultar num número muito grande de
iterações até convergir; pode convergir num mínimo local, ou, mesmo, não
convergir. Bons valores iniciais pode levar a um mínimo global, quando existir
vários mínimos locais.
SQE
b(0)
b(1)
38
2) Para o método de Gauss-Newton ou similares, é uma boa prática utilizar um
outro conjunto de valores iniciais e verificar se chega-se ao mesmo resultado.
3) Algumas propriedades válidas para os modelos lineares, não são para os
modelos não lineares. Por exemplo, a soma dos resíduos não necessariamente
é igual a zero; a soma dos quadrados do erro mais a soma dos quadrados da
regressão, não necessariamente é igual a soma dos quadrados total.
Consequentemente, o coeficiente de determinação pode não ser uma
estatística descritiva importante para os modelos não lineares.
Inferência sobre os parâmetros na regressão
não linear
Na análise de regressão não linear com erros normais, os estimadores de
mínimos quadrados ou de máxima verossimilhança, para qualquer tamanho de
amostra, não tem distribuição normal, não são imparciais e não tem variância
mínima.
As inferências sobre os parâmetros da regressão, no caso não linear,
geralmente são baseadas na teoria das grandes amostras.
39
Esta teoria mostra que os estimadores (de mínimos quadrados ou máxima
verossimilhança) para os modelos de regressão não linear com erros normais, quando
o tamanho da amostra é grande, apresentam distribuição aproximadamente normal, são
aproximadamente não tendenciosos, e aproximadamente variância mínima.
Estimativa de 2
SQE
QME 

n p
 (Y  Yˆ )
i
2
i
n p
(Y  f (X , g))


2
i
i
n p
g é o vetor das estimativas finais dos parâmetros; para os modelos de
regressão não linear, o QME não é um estimador não tendencioso de 2,
porém, o viés é pequeno se o tamanho da amostra for grande.
Teoria das grandes amostras
Teorema: para i independentes N(0,2) e o tamanho da amostra n
razoavelmente grande, a distribuição amostral de g é aproximadamente normal.
O valor esperado do vetor de médias é aproximadamente:
E ( g)  γ
(13)
40
Uma aproximação da estimativa da matriz de variância-covariância dos
coeficientes de regressão é dada por:
s2 (g)  QME(D'D)1
D é a matriz de derivadas parciais calculada nas estimativas finais, g.
Quando a teoria de grandes amostras é aplicável?
Orientações:
» o processo iterativo converge rapidamente;
» calcular algumas medidas: medidas de curvatura de Bates e Watts, medida de
vício de Box;
» estudos de simulação, por exemplo, amostragem Bootstrap verifica se as
distribuições amostrais das estimativas dos parâmetros de regressão não linear
são aproximadamente normal, se as variâncias das distribuições amostrais são
próximas das variâncias para o modelo linearizado, e se o viés em cada
estimativa dos parâmetros é pequeno.
41
Algumas medidas usadas quando os resultados da teoria das grandes amostras
não se aplica:
 Usar outra parametrização do modelo
 Fazer intervalos de confiança Bootstrap
 Aumentar o tamanho da amostra
42
Intervalo de confiança para os parâmetros
De acordo com o teorema 13, temos:
gk   k
~ t (n  p) k  0,1,2,...,p - 1
s( gk )
(14)
Onde t(n-p) é a variável com distribuição t com (n-p) graus de liberdade. De (14)
obtemos:
gk  t (1   / 2; n  p)s( gk )
Onde t(1-/2;n-p) é o (1-/2)100 percentil da distribuição t com (n-p) graus
de liberdade.
Exemplo: vamos considerar os dados de pacientes com doença grave.
Desejamos estimar 1 com um intervalo de 95% de confiança. Temos:
t (0.975;13)  2,160
g1  0,03959
s( g1 )  0,00171
 0,0433  1  0,0359
43
Concluímos, com aproximadamente 95% de confiança, que 1está entre 0,0433 e -0,0359.
Teste de hipóteses
H0 :  k   k 0
Ha :  k   k0
Onde k0 é um valor específico de k. O teste estatístico é:
gk   k 0
t 
s( g k )
*
Regra de decisão:
Se| t* | t(1   / 2; n  p), aceita - se H0 , cc rejeita - se.
Exemplo: vamos considerar os dados de pacientes com doença grave.
Desejamos testar as hipóteses:
H 0 :  0  54
H a :  0  54
44
t* 
O valor p é:
58,6065 54
 3,13
1,472
P(| t | 3,13)  0,007973
Portanto, rejeitamos a hipótese nula.
45