Gauss Moutinho CORDEIRO A ndrea A ndrade P RU DEN T E ...jaguar.fcav.unesp.br/RME/fasciculos/v27/v27_n3/A4_Andreia.pdf · exibem varia»c~ao heterog^enea e muitas vezes do tipo heterosced¶astica

UMA REVISAO D O S MO D EL O S N O RMAIS N AO -L IN EARES

Gauss Moutinho CORDEIRO1

A nd rea A nd rad e P RU DEN T E1

Claric e Garc ia B org e s DEMET RIO2

RESUMO: Os modelos normais nao-lineares continuam recebendo um tratamento

esp ecial mesmo com o surg imento dos modelos lineares g eneralizados. Esses modelos

sao ap licados as mais div ersas areas, tais como, econometria, ag ricultura, ag ronomia,

farmacolog ia, biolog ia, eng enh aria, educacao, q uımica, etc, e a p rincip al caracterıstica

dos mesmos e q ue, em g eral, sao deduz idos a p artir de sup osicoes teoricas, sendo os

p arametros resultantes interp retav eis. A o contrario dos modelos lineares, a q ualidade e,

p rincip almente, a v alidade do ajuste nos modelos nao-lineares sao av aliadas nao so p or

meios de diag nosticos de reg ressao, mas p ela ex tensao do comp ortamento nao-linear.

Modelos nao-lineares com comp ortamento distante do comp ortamento linear p odem

ter resultados assintoticos inv alidados, p rincip almente em situacoes em q ue p eq uenas

amostras sao usadas. N este artig o, tratamos da estimacao e das p rincip ais tecnicas de

diag nostico e medidas de nao-linearidade com a p resenca, tambem, de analises de dados

reais.

P A L A V RA S-C H A V E: A nalise de diag nostico; medida de v ies; medidas de curv atura;

modelos normais nao-lineares; resıduo p rojetado.

1 Introducao

A te o inıc io d a d e cad a d e 7 0 as p rinc ip ais te cnicas d esenv olv id as p ara os

m od e los d e re g re ssao nao-linear se re string iam a sup osic ao d e norm alid ad e p ara

a v ariav e l re sp osta. Em 1 9 7 2 , N e ld e r e W e d d e rb urn am p liaram a d istrib uic ao

1Departamento de Estatıstica e Informatica – DEINFO, Universidade Federal Rural de Pernambuco– UFRPE, C EP: 5 0 1 7 1 -9 0 0 , Recife, Pernambuco, B rasil. E-mail: [email protected] /deaprudente@gm ail.com

2Departamento de C iencias Ex atas, Escola S uperior de A g ricultura L uiz de Q ueiroz – ES A L Q ,Universidade de S ao Paulo – US P, C aix a Postal 9 , C EP: 1 3 4 1 8 -9 0 0 , Piracicaba, S P, B raz il. E-mail:[email protected]

3 6 0 R ev . B ras. B iom ., S ao Paulo, v.2 7 , n.3 , p.3 6 0 -3 9 3 , 2 0 0 9

da variavel resposta para a famılia exponencial de distribuicoes, definindo, entao,os modelos lineares generalizados (ML G). Mesmo assim, os modelos normais nao-lineares (MNNL ) continuam recebendo um tratamento especial, surgindo diversosartigos cientıficos na decada de 70 e nas decadas posteriores. Particularmente,destacam-se os livros de Ratk ow sk y (198 3 , 1990), q ue apresentam varios MNNL ,segundo diversos aspectos.

A forma classica do modelo normal nao-linear e expressa como

Y = g(β;X) + ε = µ(β) + ε, (1)

em q ue os εi’s sao nao correlacionados e normalmente distribuıdos com media zeroe variancia constante σ2, isto e, ε ∼ N(0, Iσ2), µ(β) = g(β;X) e uma funcaodiferenciavel em β, β = (β1, . . . , βp)

T contem os parametros desconhecidos a seremestimados e X = (x(1), . . . ,x(s)) representa a matriz , de dimensoes n×s, dos valoresde s variaveis exploratorias. No entanto, em estudos de crescimento, comuns emq uase todos os campos de pesq uisas, em geral, as medidas sao repetidas nos mesmosindivıduos caracterizando estudos longitudinais. Na maioria desses estudos, os errosexibem variacao heterogenea e muitas vezes do tipo heteroscedastica multiplicativa,pois a variancia cresce proporcionalmente com a resposta no tempo. Maioresdetalhes podem ser encontrados em Richards (195 9), S andland e McGilchrist (1979),Pinheiro e Bates (2000), Davidian e Giltinan (2003 ) dentre outros.

A principal caracterıstica do modelo (1) e q ue a parte fixa, g(β;X), em geral,decorre de um processo determinıstico deduz ido a partir de suposicoes teoricas(q uase sempre eq uacoes diferenciais), sendo os parametros resultantes interpretaveise a parte aleatoria e definida de erros homogeneos. Assim, aproxima-los paraos modelos normais lineares, mesmo q ue sejam alcancados ajustes satisfatorios,prejudicaria bastante a obtencao de estimativas mais realistas dos parametros deinteresse e, mais importante, dificultaria a interpretacao dos parametros.

Um caso particular desses modelos sao os modelos normais parcialmente nao-lineares q ue surgem muito freq uentemente na pratica e, por terem uma estruturasimples, diversas tecnicas estatısticas usuais sao simplificadas. Esses modelos saoexpressos na forma

Y = θ1g1(τ ;X) + · · · + θqgq(τ ;X) + ε = µ(β) + ε,

em q ue θ = (θ1, . . . , θq)T e o vetor de q parametros lineares, τ = (τ1, . . . , τr)

T

e o vetor de r parametros nao-lineares e β = (θT , τT )T e o vetor de p = q + rparametros, X e a matriz q ue contem os valores de s variaveis explicativas e as gj’ssao funcoes diferenciaveis. Mais formalmente, diz -se q ue um parametro θj e linearse a derivada parcial de g(β;X) em relacao a θj , ∂g(β;X)/ ∂θj , nao depende deθj ou, em outras palavras, se ∂2g(β;X)/ ∂θ2

j = 0. Em termos matriciais, isolandotermos lineares dos nao-lineares, esse modelo e escrito na forma

Y = G(τ )θ + ε,

em q ue fica definida a matriz G(τ ) = [g1(τ ;X), . . . , gq(τ ;X)].

Rev. Bras. Biom., Sao Paulo, v.27, n.3, p.360-393, 2009 361

Nem sempre os MNNL sao expressos numa forma parametrica adequada, quefacilite a convergencia rapida dos processos iterativos utilizados na estimacao dosseus parametros, sendo necessario procurar, em muitos casos, uma parametrizacaomais apropriada.

Segundo Cordeiro e Paula (1989b), embora as tecnicas de diagnostico daregressao normal nao-linear sejam simples extensoes das tecnicas da regressao linear,as interpretacoes nao sao diretamente aplicaveis, particularmente, em virtude de osresıduos ordinarios nao terem mais uma distribuicao, aproximadamente, normal.Isso levou ao desenvolvimento de tecnicas especıficas de diagnostico para os MNNL(Cook e Tsai, 1985). Similarmente, as propriedades das somas de quadradoscontidas nas tabelas classicas de analise da variancia (ANOV A) nao sao estendidasdiretamente para o caso nao-linear. Entretanto, alguns pesquisadores continuamconstruindo tais tabelas apos o ajuste do MNNL e utilizam apenas, descritivamente,os valores obtidos para a estatıstica F , razao do quadrado medio da regressaonao-linear pelo quadrado medio do resıduo, calculada pela equacao (5).

Os MNNL sao aplicados nas mais diversas areas, tais como, econometria,agricultura, agronomia, farmacologia, biologia, ecologia, engenharia, educacao,quımica, etc. Na Tabela 1, sao apresentados alguns modelos especiais.

Tabela 1 - Alguns modelos normais nao-lineares

Modelo Componente sistematico

Gallant β1x1 + β2x2 + β3 exp(γx3)

V on-Bertalanff y α{1 − exp[−γ(x − δ)]}

Mitscherlich α[1 − 10−γ(x−δ)]

Logıstico α/{1 + exp(β − γx)}

Gompertz α exp{− exp(β − γx)}

Richards α/[{1 + exp(β − γx)}1/ δ]

Morgan-Mercer-F lodin (MMF ) (βγ + αxδ)/(γ + xδ)

Weibull α − β exp(−γxδ)

Stone α + δ log[x1/(γ + x2)]

Michaelis-Menten αx/(γ + x)

No modelo de Gallant, x1 e x2 representam dois tratamentos e x3 o tempo,que afeta, exponencialmente, a media da variavel resposta. No modelo de V on-Bertalanff y, frequentemente usado na area ecologica, para explicar o comprimentode um peixe, α e o comprimento maximo esperado para a especie, γ e a taxa mediade crescimento e δ e um valor nominal em que o comprimento do peixe e zero. Nomodelo de Mitscherlich, usado em agricultura, α e a producao maxima esperada,γ e o coeficiente de eficacia do adubo e δ e o teor de nutriente existente no solosem adubacao. Nos modelos, logıstico, de Gompertz, de Richards, Morgan-Mercer-

362 Rev. Bras. Biom., Sao Paulo, v.27, n.3, p.360-393, 2009

Flodin (MMF) e Weibull α e a assıntota superior, γ e a taxa media de crescimentoe β esta relacionado com o intercepto. O parametro δ, que aparece nos modelosde Richards, Morgan-Mercer-Flodin (MMF) e Weibull, e utilizado para aumentar afl exibilidade dos mesmos no ajuste dos dados. No modelo de Stone, usado parao estudo de mistura de duas drogas, z1 e z2 representam, respectivamente, asconcentracoes de uma droga ativa e de um reagente, δ e a inclinacao comum darelacao log-dose-resposta. No modelo de Michaelis-Menten, usado no estudo develocidade de reacoes quımicas, x e a concentracao do substrato, α e a velocidademaxima da reacao e γ e a constante de Michaelis.

Na Secao 2, apresentam-se o processo iterativo de Newton-Raphson paracalcular a estimativa de mınimos quadrados de β e alguns resultados assintoticos deinferencia. As principais tecnicas de diagnostico usadas para modelos de regressaonormal nao-linear sao apresentadas na Secao 3. Medidas de nao-linearidade,correcao do vies e melhoramento dos testes sao abordados na Secao 4 . Finalmente,algumas aplicacoes sao descritas na Secao 5.

2 Estimacao

Sejam Y1, . . . , Yn variaveis aleatorias independentes com a estruturaespecificada em (1) e suponha que as observacoes sejam representadas pelo vetory = (y1, . . . , yn)T . A estimativa do vetor de parametros β e calculada pelometodo dos mınimos quadrados que, nesse caso, coincide com o metodo de maximaverossimilhanca, pois o modelo tem respostas normais independentes com a mesmavariancia. A estimativa de β e, entao, calculada, minimizando-se a funcaoquadratica

S(β) =

n∑

i= 1

{yi − µi(β)}2.

Derivando-se S(β) em relacao a βr, obtem-se

∂S(β)

∂βr= 2

n∑

i= 1

{yi − µi(β)}∂µi

∂βr.

A estimativa β do vetor de parametros β e obtida igualando-se∂S(β)

∂βra zero

para r = 1, . . . , p. Em geral, as equacoes∂S(β)

∂βr= 0 nao sao lineares e tem que ser

resolvidas, numericamente, por processos iterativos do tipo Newton-Raphson.Considerando-se que se deseja obter a solucao do sistema de equacoes

∂S(β)/∂β = 0 e, usando-se a versao multivariada do metodo de Newton-Raphson,tem-se

β(m+ 1) = β(m) + (J(m))−1 ∂S(β)

∂β

(m)

,


sendo β(m) e β(m+1) os vetores de parametros estimados nos passos m e (m + 1),respectivamente, {∂S(β)/∂β}(m) o vetor de derivadas de primeira ordem avaliadono passo m, e (J(m))−1 a inversa da negativa da matriz de derivadas parciais desegunda ordem de S(β), com elementos −∂2S(β)/∂βr∂βs, avaliada no passo m.Portanto,

β(m+1) = β(m) + {X(m)T X(m)}−1X(m)T {y − µ(β(m))}, (2)

m = 0, 1, . . . , em que X = {∂µi(β)/∂βr} e a matriz J acobiana da transformacaode µ(β) em β, tem posto completo p e pode, tambem, ser denominada de matrizmodelo local. O processo (2) e repetido ate obter a convergencia, definindo-se,

entao, β = β(m+1).Dentre os muitos existentes, um criterio para verificar a convergencia poderia

ser

p∑

r=1

(β

(m+1)r − β

(m)r

β(m)r

)2

< ξ,

tomando-se para ξ um valor suficientemente pequeno. A convergencia de (2), emgeral, depende dos valores iniciais para os parametros do vetor β. Isso podeevitar que problemas relacionados com a estrutura parametrica do modelo, taiscomo, nao-linearidade acentuada (Secao 4) e/ ou mau condicionamento da matriz

X, prejudiquem a convergencia do processo iterativo. Em Sousa (1986 ), ha umadescricao detalhada do metodo de Newton-Raphson e de outros metodos iterativosusuais em regressao normal nao-linear. Ratkowsky (1983) sugere algumas tecnicaspara calcular valores iniciais para os parametros de β para os modelos descritosna Secao 1. A estimacao dos parametros β por metodos Bayesianos pode serencontrada em Eaves (1983) e Gelman et al (2004), dentre outros.

2 .1 Resu ltados assintoticos

O logaritmo da funcao de verossimilhanca do modelo (1), como funcao de β,e expresso na forma

`(β) = (2π σ2)−n/2 exp{−S(β)/2σ2}.

A estimativa de maxima verossimilhanca β e obtida pelo processo iterativo (2),sendo consistente e tendo, assintoticamente, uma distribuicao normal p variada demedia β e estrutura de variancia-covariancia K−1 = σ2(XT X)−1 (J ennrich, 196 9).Analogamente a regressao linear, a estimativa mais usual para σ2 e expressa pors2 = S(β)/(n − p), em que S(β) e a soma de quadrados dos resıduos do modeloajustado. Logo, um intervalo de confianca 100(1−α)% para βj , sera formado peloslimites

βj ± tα /2 × (κjj)1/2,


em que tα/2 e o quantil (1−α/2) de uma distribuicao t de Student com n− p grausde liberdade e κjj e a estimativa do elemento (j, j ) de K−1.

Uma regiao de, aproximadamente, 100(1 − α)% de confianca para β foiproposta por Beale (1960), sendo formada pelos contornos de S(β) tais que

S(β) = S(β){

1 +p

n − pFp,n−p(α)

}. (3)

Em particular, se `(β) for aproximadamente quadratica, a regiao de confianca (3)e bem aproximada por

(β − β)T (X

T X)(β − β) ≤ s2pFp,n−p(α),

em que Fp,n−p(α) e o quantil (1 − α) de uma distribuicao F eX e a matriz X

avaliada em β. Essa ultima expressao e uma adaptacao da regiao de confianca daregressao normal linear.

Para testar a hipotese H : β ∈ B, em que B e um subconjunto do espacoparametrico, adota-se, usualmente, a estatıstica da razao de verossimilhancas,expressa por

w = n log{S(β)/S(β)}, (4)

em que S(β) e a soma de quadrados dos resıduos para o modelo ajustado em H.Sob essa hipotese, a estatıstica (4) tem, assintoticamente, distribuicao χ2 com p−mgraus de liberdade, em que m = dim(B). Johansen (1983) mostra que a estatıstica(4) e, assintoticamente, equivalente a estatıstica

n

n∑

i=1

{µi(β) − µi(β)}2/S(β),

que e mais facil de ser calculada. A aproximacao qui-quadrado para a distribuicaonula de (4) pode ser melhorada pela correcao de Bartlett (Cordeiro, 1983, 1987;Cordeiro e Paula, 1989a).

Uma estatıstica alternativa para testar H e expressa como

F =n − p

p − m

S(β) − S(β)

S(β), (5)

que sob a hipotese H tem, assintoticamente, distribuicao F com p − m e n − pgraus de liberdade. Logo, deve-se rejeitar H, para um nıvel de significancia α,se F ≥ Fp−m,n−p(α), em que Fp−m,n−p(α) e o ponto crıtico da distribuicao Fcorrespondente. Esse resultado vale, tambem, quando a variavel resposta nao enormal, havendo contudo algumas condicoes adicionais de regularidade.

Segundo Souza (1998) a adequacidade do ajuste de um modelo nao-linear a umconjunto de observacoes pode ser medida pelo quadrado do coeficiente de correlacaoentre os valores observados e preditos, calculada com a utilizacao da formula:


R2 = r2yy =

[∑n

i=1 yiyi −( n

i=1yi)(

n

i=1yi)

n

]2

[∑n

i=1 y2i −

( n

i=1yi)2

n

][∑n

i=1 y2i −

( n

i=1yi)2

n

] = 1 −SQ R

SQ T o ta lc(6)

em que SQ R e a soma de quadrados do resıduo e SQ T o ta lc e a soma de quadradostotal corrigida pela media. Porem, a ultima igualdade de (6) so se verifica se amatriz Jacobiana tiver uma coluna de uns. Vale ressaltar que existe muito uso eabuso do R2 em modelos de regressao nao-linear, pois podem-se ter modelos comR2 = 0, 99 e haver uma enorme discrepancia entre os valores observados e preditos.

3 Tecnicas de diagnostico

Exceto com relacao aos resıduos, as tecnicas mais usuais de diagnostico emregressao normal nao-linear sao simples adaptacoes da regressao linear. Algumasdessas tecnicas serao apresentadas nesta secao.

3.1 Matriz de projecao

No modelo normal nao-linear, para a deteccao de pontos mais afastados dosdemais, possivelmente, pontos influentes, utiliza-se a matriz de projecao ortogonalno subespaco tangente a µ, expressa como

H = {hij} = X(XT X)−1XT ,

em que X deve ser avaliada em β. Ao contrario da regressao linear, essa e umamatriz de projecao local, pois depende de β e deve ser estimada em β. Mesmoassim, o criterio hii ≥ 2p/n continua sendo adotado como guia para detectar pontossuspeitos de serem influentes.

3.2 Resıduo projetado

Segundo Cordeiro e Lima Neto (2004), os resıduos ordinarios no modelo

nao-linear sao definidos por ri = yi − µi(β), i = 1, . . . , n. A distribuicao dessesresıduos e matematicamente intratavel, principalmente para pequenas amostras.Alem disso, esses resıduos, em geral, tem esperanca diferente de zero e distribuicaodependendo fortemente dos valores ajustados, o que pode conduzi-los a naorefletirem exatamente a distribuicao dos erros. Logo, nesse caso, os criterios dediagnostico da regressao normal linear podem falhar. Por exemplo, um resıduomuito diferente de zero, que segundo os criterios da regressao linear seria umponto aberrante, pode agora nao ser, caso seu valor esperado seja, tambem,substancialmente diferente de zero.


Define-se, a seguir, um novo resıduo, que apesar de algebricamente maiscomplexo, tem propriedades mais proximas das propriedades correspondentes doresıduo ordinario da regressao normal linear.

Seja C(X) o subespaco gerado pelas colunas de X. De agora em diante,

usa-se a notacao C(X) para representar o subespaco ortogonal a C(X), ou seja, o

subespaco gerado pelas colunas ortogonais aquelas de X.Definem-se as matrizes p × p, Wi = {∂2µi/∂βr∂βs} de derivadas de segunda

ordem para i = 1, . . . , n. Define-se ainda W como uma matriz n × (p × p) com nfaces cuja i-esima face e, simplesmente, igual a Wi.

Expandindo µ(β) em serie de Taylor em torno de β ate segunda ordem, Cooke Tsai (1985) apresentaram a seguinte aproximacao para r

r ∼= (I − H)ε − X

n∑

i=1

riWiδ −1

2(I − H)δT Wδ, (7)

em que δ = β − β. Na expansao (7), δT Wδ representa δT Wiδ para produzir ai-esima componente de r.

Uma aproximacao quadratica para r e obtida substituindo a aproximacaolinear para r ∼= (I − H)ε e δ ∼= XT ε no segundo e terceiro termos da expansao(7), implicando que

E(r) ∼= (I − H)f

eCov(r, µ(β)) ∼= NNT σ2 − Var(r),

em que f e um vetor n × 1 de elementos fi = − 12σ2tr(Wi), i = 1, . . . , n, N e

uma matriz n×n cujas colunas formam uma base ortonormal em C(X) (subespaco

gerado pelas colunas ortogonais a X) e Var(r) = NNT σ2+ parte positiva. Logo, a

covariancia entre r e µ(β) tende a ser negativa, o que pode dificultar a interpretacaodos graficos padroes baseados em r.

Sejam os vetores colunas n × 1 nao-nulos calculados como:

x11 =

(∂2µ1

∂β21

, . . . ,∂2µn

∂β21

)T

, x12 =

(∂2µ1

∂β1∂β2, . . . ,

∂2µn

∂β1∂β2

)T

, . . . ,

xpp =

(∂2µ1

∂β2p

, . . . ,∂2µn

∂β2p

)T

.

Existem p(p + 1)/2 desses vetores, embora considerem-se apenas os q nao-nulos(q ≤ p(p + 1)/2).

Define-se a matriz W∗ como uma matriz n × (p × p) que pode ser calculada

de W projetando cada coluna de derivadas de segunda ordem xjk sobre C(X), isto

e, sobre o subespaco complementar gerado por X. Logo, se T = (x11, . . . , xpp) e

uma matriz n × q e S = (I − H)T, com H = X(XT X)−1XT , tem-se que H1 =S(ST S)−1ST e o operador de projecao ortogonal em C(W∗). Assim, I − H1 e ooperador de projecao ortogonal em C(W∗).


Seja agora uma matriz V do tipo n × (p + q) definida como V = (X,S). Ooperador de projecao ortogonal em C(V) e, simplesmente, H2 = V(VT V)−1VT e,entao, I − H2 e o operador de projecao ortogonal em C(V).

Utilizando (7), Cook e Tsai (1985) definiram o resıduo projetado

(I − H2)r = (I − H)ε − (I − H1)ε. (8)

O primeiro termo de (8) e a aproximacao linear para o resıduo ordinario r,enquanto o segundo termo reflete a perda de informacao necessaria para remover oscomponentes nao-lineares de (7). Se q = posto(H1) for pequeno em relacao a n−p,entao essa perda, tambem, sera pequena. Independente disso, se a medida de nao-linearidade intrınseca γI N , descrita detalhadamente na Secao 4, for significativa,isto e, γI N > 2F−1/2, o ganho sera substancial.

De (8) mostra-se, facilmente, que

E{(I − H2)r} = 0, Var{(I − H2)r} = σ2(I − H2)

eE{rT (I − H2)r} = σ2tr(I − H2).

Logo, uma estimativa alternativa para σ2 e expressa por

σ2 =rT (I − H2)r

tr(I − H2).

Os resıduos projetados superam os resıduos ordinarios em diversos aspectose muitas das tecnicas de diagnostico utilizadas na regressao linear sao, tambem,aplicaveis aos mesmos. Por exemplo, os graficos de (I − H2)r versus variaveisexplanatorias nao incluıdas no modelo podem revelar como esses termos aparecemno componente sistematico.

E importante lembrar que os operadores utilizados dependem de β e,portanto, na pratica e preciso substituir essas quantidades pelas respectivasestimativas. Claramente, r esta em C(X), quando X e avaliado em β; logo,

(I − H2)r = (I − H1 − H)r = (I − H1)r, sendo H1r os valores ajustados da

regressao linear de r sobre (I − H)xkj , k , j = 1, . . . , p.Na regressao linear, mesmo para erros nao-correlacionados e de variancia

constante, os resıduos sao correlacionados e tem variancias diferentes. Definem-se,entao, os resıduos estudentizados que mesmo correlacionados, apresentam mediazero e variancia constante e igual a um. O i-esimo resıduo ordinario estudentizadoe expresso por

ti =ri

s(1 − hii)1/2, i = 1, . . . , n.

Os resıduos projetados estudentizados tem esperanca nula e variancia σ2(I −H2) e podem ser, entao, definidos, como

si ={(I − H2)r}i

σ{(I − H2)}1/2ii

, i = 1, . . . , n. (9)


Cook e Tsai (1985) mostram, para um exemplo particular, o grafico de (ti−si)versus os valores de uma unica variavel explanatoria e concluem que diferentesdiag nosticos sao ob tidos se os criterios utilizados para si forem, tamb em, adotadospara os resıduos ordinarios estudentizados ti, i = 1 , . . . , n .

P ara avaliar se os erros εi’s tem distrib uicao, aproximadamente, normal,assim como para detectar se h a pontos ab errantes e/ ou infl uentes, o g rafi co de

prob ab ilidades dos resıduos projetados ordenados s(i) versus Φ−1( i−3/8

n+ 1/4

)pode ser

util, sendo Φ(.) a funcao de distrib uicao acumulada da normal padrao. A analisedos resıduos em (9 ) procede-se, similarmente, ao modelo normal linear.

3.3 Medidas de influencia

A s medidas de infl uencia para o modelo normal nao-linear sao b aseadas nareg ressao linear. A unica diferenca, que pode ser relevante, e a sub stituicao

da estimativa β(i) pela estimativa correspondente β1

(i), que e ob tida iniciando o

processo iterativo (2 ) em β sem a i-esima ob servacao e considerando a estimativade um passo. C omo o metodo de N ew ton-R aph son utiliza em cada passo uma

aproximacao quadratica para `(β), a estimativa β1

(i) pode nao estar muito proxima

de β(i), se `(β) nao for, localmente, quadratica. E ntretanto, varios estudos desimulacao tem mostrado que essa aproximacao e sufi ciente para detectar os pontosinfl uentes.

C ook e W eisb erg (1 9 8 2 ) mostraram que essa estimativa de um passo pode serexpressa como

β1

(i) = β −(XT

X)−1

(1 − hii)xiri,

em que xi e a i-esima linh a de X. P ara calcular o valor numerico de β1

(i), a matriz

X deve ser avaliada em β. L og o, β1

(i) depende de quantidades correspondentes aoi-esimo ponto e de quantidades conh ecidas que envolvem todas as ob servacoes.

A distancia de C ook e expressa por

Di = (β − β(i))T (XT

X)(β − β(i))/ ps2, (1 0 )

em que s2 foi defi nido na S ecao 2 .1 . U sando as formulas anteriores de ti e β1

(i) naexpressao (1 0 ), ob tem-se a forma aproximada

D1i =

t2ip

hii

(1 − hii)para i = 1 , . . . , n .

O s criterios de calib racao para a reg ressao normal linear podem serestendidos para o modelo normal nao-linear desde que os contornos de S(β) =∑n

i= 1{yi − µi(β)}2 sejam, aproximadamente, elıpticos. Isso porque em muitos

Rev. B ra s. B io m ., S a o P a u lo , v .2 7 , n .3 , p .3 6 0 -3 9 3 , 2 0 0 9 3 6 9

problemas de regressao normal nao-linear, as regioes de confianca usuais paraβ podem ser seriamente viesadas (B eale, 196 0) e o vies pode depender daparametrizacao escolhida (B ates e Watts, 1980). Logo, escolher uma parametrizacaoadequada pode ser importante para detectar pontos influentes.

O grafico de D1i versus a ordem das observacoes e um procedimento usual,

devendo-se dar atencao aqueles pontos com o D1i correspondente mais afastado dos

demais. Se o interesse e detectar pontos influentes nas estimativas individuais βj ,

j = 1, . . . , p, sugere-se o grafico de ∆ iβj = (βj − β(i)j)/DP (βj) versus a ordem dasobservacoes.

3.4 Grafico da variavel adicionada

O grafico da variavel adicionada pode revelar como as observacoes,conjuntamente, estao influenciando na obtencao da estimativa do parametro queesta sendo incluıdo no modelo. G iltinan et al. (1988) mostraram que esse graficopode ser estendido para a classe dos modelos normais nao-lineares, mas de umaforma um pouco diferente. No modelo normal nao-linear, faz mais sentido incluirum novo parametro na parte sistematica, que em muitos casos pode significar umainteracao, do que incluir uma nova variavel exploratoria.

Suponha, entao, a media nao-linear µ(β) para o modelo reduzido e o preditornao-linear µ(β, γ) com um parametro escalar γ a ser incluıdo no modelo. Seja xγ

um vetor n× 1 com as derivadas parciais de µ(β, γ) em relacao a γ. G iltinan et al.

(1988) sugerem o grafico de r = y − µ(β) versus (I − H)xγ , em que H e a matrizde projecao correspondente ao modelo reduzido e xγ e o vetor xγ computado sob ahipotese H : γ = 0. A estimativa γ e igual a estimativa do parametro da regressaolinear simples, passando pela origem, de y−µ(β) sobre (I− H)xγ . Logo, o graficoproposto pode revelar como as observacoes estao contribuindo para essa relacao ecomo estao se afastando dela.

4 Medidas de nao-linearidade

O principal objetivo das medidas de nao-linearidade e verificar se o grau denao-linearidade de um problema de estimacao nao-linear e suficientemente pequenopara que as tecnicas usuais de estimacao, desenvolvidas para a regressao linear,sejam utilizadas como uma boa aproximacao para o modelo nao-linear.

A primeira tentativa relevante no sentido de desenvolver uma medida de nao-linearidade foi de B eale (196 0). Entretanto, G uttman e M eeter (196 5 ) mostraramque a medida proposta por B eale tende a subestimar o verdadeiro grau de nao-linearidade do modelo. Uma outra contribuicao importante foi a de B ox (197 1)que calculou o vies de ordem n−1 do estimador de maxima verossimilhanca (EM V )do vetor β de um modelo normal nao-linear. Entretanto, foi somente no inıcioda decada de 80, que surgiu o trabalho mais relevante nessa area. B ates e Watts(1980), utilizando alguns conceitos de geometria diferencial, desenvolveram duas


medidas de curvatura para os modelos normais nao-lineares. Essas medidas indicam,respectivamente, o grau de nao-linearidade intrınseca de um modelo e o grau denao-linearidade aparente ou devida a parametrizacao utilizada.

Ratkowsky (1983) comparou algumas formas parametricas para diversosmodelos normais nao-lineares por meio de simulacoes e utilizou as medidas de Box ede Bates e Watts. Para se ter uma ideia mais clara dos conceitos de nao-linearidadeintrınseca e de nao-linearidade aparente, serao comparados, a seguir, um modelolinear e um modelo nao-linear para o caso de n = 2 e p = 1.

Considere, inicialmente, o modelo linear simples Yi = βxi + εi, i = 1, 2, emque x e uma covariavel qualquer e β um parametro desconhecido. Nesse caso, oespaco de estimacao tem dimensao igual a um, sendo formado pelos pontos

Xβ =

(x1

x2

)β, β ∈ R,

ou seja, e uma reta em R2. Alem disso, para qualquer conjunto de solucoes β(1), β(2),. . ., tais que β(i+1) − β(i) = ∆, em que ∆ e uma constante arbitraria, as solucoespossıveis para Xβ serao tais que

Xβ(i+1) − Xβ(i) =

(x1

x2

)∆, i = 1, 2, . . . ,

ou seja, se as solucoes para β forem igualmente espacadas, entao os valores ajustadoscorrespondentes serao, tambem, igualmente espacados.

Considere agora o modelo normal nao-linear Yi = xβi + εi, i = 1, 2 e os dados

apresentados em Ratkowsky (1983, pag. 7)

y = (2, 5 10)T e X = (2 3)T .

Nesse caso, o espaco de estimacao nao e mais uma reta, e sim uma curva aoredor da estimativa de maxima verossimilhanca β = 2, 05. A curva correspondenteaos pontos (2β 3β)T com β variando em espacamentos iguais a 0, 5 e apresentadana F igura 1. Note que os pontos do espaco de estimacao nao sao igualmenteespacados como ocorre no modelo linear.

Assim, quanto mais essa curva se afasta da reta tangente em β maior sera o queBates e Watts (1980) chamam de “ nao-linearidade intrınseca” do modelo, e quantomais desiguais forem os espacamentos entre os pontos do espaco de estimacao,maior sera o que ambos autores chamam de “ nao-linearidade aparente” causadapela parametrizacao do modelo.

Portanto, a nao-linearidade de um modelo pode ser devida a duas causas. Aprimeira e a curvatura real do modelo ou intrınseca, como definem Bates e Watts,que e invariante com qualquer tipo de reparametrizacao. A segunda e a curvaturadevida a forma como os parametros aparecem no modelo. Essa ultima pode sereliminada ou pelo menos reduzida atraves da reparametrizacao. Para ilustrar essefato, considere o modelo normal nao-linear descrito anteriormente com a seguintereparametrizacao:

Yi = xlo g φi + εi, i = 1, 2,


Figura 1 - Representacao da curva (2β 3β)T com β variando em espacamentosiguais a 0,5 (Ratkowsky, 1983).

em que φ = exp(β). A Figura 2 mostra os pontos da curva (2log φ 3log φ)T comespacamentos iguais a 1,0 para φ. Nota-se que os espacamentos entre os pontoscorrespondentes sao praticamente iguais, indicando que o grau de nao-linearidadeaparente foi, substancialmente, reduzido com essa reparametrizacao. Entretanto, acurvatura do espaco de estimacao continua com a mesma forma anterior, como erade se esperar.

4.1 Medidas de curvatura de Bates e Watts

Considere o modelo de regressao normal nao-linear definido na Secao 1.Uma reta no espaco parametrico passando por β, pode ser expressa, usando umparametro escalar b, por

β(b) = β + bh,

em que h = (h1, . . . , hp)T e um vetor de valores nao-nulos. Essa reta gera uma

curva, sobre o espaco de estimacao, definida por

µh(b) = µ(β + bh).

A tangente a essa curva no ponto b = 0 e expressa na forma

µh

= Xh, (11)

em que X e a matriz J acobiana da transformacao de µ(β) em β = β. O conjuntode todas as combinacoes lineares da forma (11) e, tambem, denominado de plano

tangente em µ(β).


Figura 2 - Representacao da curva (2log φ 3log φ)T com φ variando em espacamentosiguais a 1,0 (Ratkowsky, 1983).

A aceleracao da curva µh e estimada por

µh = hT Wh,

em que W e uma matriz de dimensao n × (p × p) com i-esima face expressa comoWi = (∂2µi/∂βr∂βs), i = 1, . . . , n e r, s = 1, . . . , p. Portanto, cada elemento

indexado por i e da forma hT Wih, para i = 1, . . . , n.O vetor de aceleracao µh pode ser decomposto em tres componentes. O

primeiro componente µI Nh

determina a variacao na direcao do vetor de velocidadeinstantanea µh normal ao plano tangente, enquanto o segundo e o terceirocomponentes, cuja norma sera representada por µP E , determinam, respectivamente,a variacao na direcao de µh paralela ao plano tangente e a variacao na velocidade doponto movel. Esses componentes foram transformados, por Bates e Watts (1980),nas seguintes curvaturas:

A - Curvatura intrınseca – definida por

K I Nh

= ‖µI N ‖/‖µh‖2.

B - Curvatura devida a parametrizacao – definida por

K P Eh

= ‖µP E ‖/‖µh‖2.

Essas curvaturas podem ser padronizadas de modo que fiquem invariantesa mudancas de escala. Isso e obtido multiplicando K I N

he K P E

hpor s

√p com

s = {S(β)/n − p}1/2. T em-se, portanto, as curvaturas padronizadas

γI Nh = s

√p K I N

he γP E

h = s√

p K P Eh

. (12)


As medidas relativas (12) podem ser usadas nao somente para comparardiferentes parametrizacoes de um determinado modelo, mas, tambem, diferentesconjuntos de dados para o mesmo modelo ou para modelos diferentes.

As medidas de nao-linearidade de Bates e Watts (1980) sao definidas comosendo as curvaturas maximas

γIN = maxh{KINh } e γPE = maxh{KPE

h }.Bates e Watts sugerem o criterio

γIN ≥ 2F−1/2 e γPE ≥ 2F−1/2

como guia para indicar se o modelo ajustado tem, respectivamente, curvaturaintrınseca e curvatura aparente acentuada, em que F e o quantil (1 − α) de umadistribuicao F com p e n − p graus de liberdade.

Para o calculo dessas medidas e preciso, inicialmente, decompor a matriz X

num produto de duas matrizes Q e R, isto e, X = QR, sendo Q uma matriz n× nortogonal e R uma matriz n × p definida por

R =

[R

0

],

sendo R uma matriz p × p triangular superior e inversıvel. As matrizes Q e R

podem ser calculadas a partir da decomposicao de Businger e Golub (1965).

A seguir, deve-se obter a matriz U = LT WL, sendo L = R−1. Os elementosde U sao vetores n×1 representados por Uk j , k , j = 1, . . . , p. D efine-se, entao, o queBates e Watts chamam de matriz de aceleracao A = QT U de dimensao n × p × p.O (k , j)-esimo elemento dessa matriz e um vetor n × 1 expresso na forma QT Uk j .A matriz A pode ser decomposta como

A =

QT U11 · · · QT U1p

...QT Up1 · · · QT Upp

,

em que QT Uk j = (ak j1, . . . , ak jn)T . A i-esima face de A e expressa na forma

Ai =

a11i · · · a1pi

...ap1i · · · appi

, i = 1, . . . , n.

Sejam APE a matriz constituıda das p primeiras faces de A e AIN a matrizconstituıda das ultimas n − p faces de A. Entao, as medidas de nao-linearidadeserao formuladas como

γIN = maxh‖hT AINh‖ e γPE = maxh‖hT APEh‖,sendo ‖h‖ = 1. Para efetuar os calculos nao ha, em geral, formulas explıcitas, sendonecessario recorrer a algum processo iterativo. Sousa (1986) descreve a obtencao deγIN e γPE por meio de um processo iterativo proposto por Bates e Watts (1980).


4.2 Vies de ordem n−1 de β

Cox e Snell (1968) deduziram o vies de ordem n−1 do estimador de maximaverossimilhanca do vetor de parametros β em uma classe geral de modelos queinclui o modelo normal nao-linear como um caso especial. Box (1971), utilizando

esse trabalho, obteve uma aproximacao b ∼= E(β−β) expressa, em forma matricial,como

b = (XT X)−1XT d, (13)

em que d e um vetor n×1 com elementos di = − 12σ2tr{(XT X)−1Wi}, i = 1, . . . , n

e, como antes, Wi = (∂2µi/∂βr∂βs) e uma matriz quadrada de ordem p. Portanto,o vies b e, simplesmente, a estimativa de mınimos quadrados da regressao normallinear de d sobre as colunas de X. Aqui, X depende do vetor β de parametrosverdadeiros. Entretanto, para se estimar b devem-se avaliar X e d em β.

Bates e Watts (1980) mostraram que o vies de Box esta relacionado coma medida de nao-linearidade γPE . Portanto, o vies pode ser reduzido por meiode alguma reparametrizacao do modelo e a expressao (13) pode indicar quaisparametros sao os maiores responsaveis por um valor apreciavel de nao-linearidade.

Cook et al. (1986) mostraram que d e igual a diferenca entre os valoresesperados das aproximacoes linear e quadratica para µ(β). Logo, o vies sera pequenose todos os elementos de d forem suficientemente proximos de zero, o que indicaque o modelo e, essencialmente, linear ou se d e ortogonal as colunas de X.

Em particular, para os modelos normais parcialmente nao-lineares definidospor

Yi = xTi β + δ fi(γ) + εi, (14 )

em que xi e a i-esima linha da matriz X, de dimensoes n × (p − 2), associadaaos parametros lineares em β = (β1, . . . , βp−2)

T e δ e γ sao parametros escalares,tem-se,

X = [X, f(γ), δf ′(γ)],

em que f(γ) = (f1(γ), . . . , fn(γ))T e f ′(γ) = (f ′

1(γ), . . . , f ′

n(γ))T . De (13), obtem-se

di = 2f ′

i(γ)Cov(δ, γ) + δf ′′

i (γ)Var(γ), i = 1, . . . , n.

Logo, o vies fica expresso por

b = −δ−1Cov(δ, γ)ρp − δ

2Var(γ)(XT X)−1XT τ ,

em que ρp e um vetor p × 1 de zeros com o valor um na ultima posicao e

(XT X)−1XT τ representa o vetor das estimativas dos coeficientes da regressao

normal linear de τ = (f ′′

1 (γ), . . . , f ′′

n (γ))T sobre X. A covariancia Cov(δ, γ)contribui somente para o vies de γ.


Box (1971), tambem, desenvolveu uma formula para avaliar o vies dosestimadores em uma nova parametrizacao, mostrando que o novo vies pode sercalculado por meio do vies da parametrizacao anterior.

Considere a reparametrizacao

φ = g(β),

sendo φ um escalar, g(.) uma funcao diferenciavel e β = (β1, . . . , βp)T . Seja bφ o

vies de ordem n−1 de φ. Box mostrou que

bφ = gT b +1

2tr{M Cov(β)},

em que g e um vetor p × 1 das derivadas parciais de primeira ordem de g(β)em relacao a β e M e uma matriz p × p de derivadas parciais de segunda ordem∂2g(β)/∂βr∂βs, r, s = 1, . . . , p. As duas quantidades g e M sao avaliadas em β.

A variancia de φ pode, tambem, ser expressa em funcao da covariancia de β

por

Var(φ) = tr{(ggT ) Cov(β)}.

Em particular, para p = 1

bφ = bdg(β)

dβ+

1

2Cov(β)

d2g(β)

dβ2

e

Var(φ) = Var(β){d2g(β)

dβ2

}2

,

sendo todas as derivadas avaliadas em β.

4.3 Ap erfeicoamento da estatıstica da raz ao de verossimilhancas

Cordeiro e Paula (1989a), utilizando a metodologia da correcao de Bartlett(1937) e as expansoes de Lawley (1956), corrigiram a estatıstica da razao deverossimilhancas w, ate ordem n−1, para a classe dos modelos exponenciais nao-lineares, que engloba como caso especial os modelos normais nao-lineares. Esse fatorde correcao, denotado por c, faz com que a estatıstica corrigida c−1w aproxime-semelhor da distribuicao qui-quadrado de referencia do que a estatıstica usual w.

Para ilustrar, suponha a particao β = (βTq , βT

p−q)T , p > q e a hipotese nula

H : βp−q = 0. Nesse caso, a estatıstica da razao de verossimilhancas e igual a

w = 2(ˆp − ˆq), em que ˆ

q e ˆp sao, respectivamente, o logaritmo da funcao de

verossimilhanca maximizada sob H e sob o modelo em pesquisa. A correcao deBartlett pode ser expressa como

c = 1 + (εp − εq)/(p − q), (15)


em que εp e um termo bastante complicado envolvendo valores esperadosde produtos de derivadas (cumulantes conjuntos) do logaritmo da funcao deverossimilhanca. Em especial, para os modelos normais nao-lineares, tem-se(Cordeiro e Paula, 1989a)

εp =σ2

4ψ =

σ2

4{2tr(Bd − BH) − 1T DMD1}, (16)

em que H = X(XT X)−1XT , B = {bij} com bij = tr{Wi(XT X)−1Wj(X

T X)−1},D = diag{d1, . . . , dn}, 1 e um vetor n × 1 de 1’s, M = I − H e o operador deprojecao ortogonal de vetores de Rn no subespaco gerado pelas colunas da matrizX e Bd e uma matriz diagonal de ordem n com os elementos da diagonal de B.

Mostra-se, utilizando (16), que

E

(SQ R e s

σ2

)∼= n − p − σ2

4ψ,

isto e, o valor esperado da soma de quadrados dos resıduos dividido por σ2 e,aproximadamente, igual a n−p, que e o valor esperado no modelo linear, mais umacontribuicao devida a nao-linearidade em µ(β), multiplicada por −σ2/4. Johansen(1983) relacionou ψ com a medida de nao-linearidade de Beale (1960).

Restringindo-se a subclasse dos modelos parcialmente nao-lineares (14), εp

reduz-se a

εp =σ2

4Var2(γ)1T ΓMΓ1,

em que Var(γ) e a variancia assintotica de γ, Γ = δ diag{f ′′

1 (γ), . . . , f ′′

n (γ)} e1T ΓMΓ1 e a soma de quadrados dos resıduos na regressao linear de Γ1 sobreas colunas de X. A correcao de Bartlett sera apreciavel quando Var(γ) for grande.

Na pratica, a correcao c deve ser estimada sob o menor modelo, isto e, tanto εp

quanto εq em (15) devem ser calculados sob H : βp−q = 0. Para ilustrar, suponha

que num modelo parcialmente nao-linear, o interesse e testar H : γ = γ(0 ). Logo, acorrecao de Bartlett e expressa por

c = 1 +σ2

4Var2(γ)1T ΓMΓ1,

em que as quantidades Var(γ), Γ e M devem ser avaliadas sob H. Aqui, em

particular, Var(γ) e o elemento que ocupa a posicao (p, p) da matriz (XT X)−1.

5 Aplicacoes

5 .1 Modelo para explicar a fracao media de cloro disponıvel num

produto

Os dados da Tabela 2 e representados na Figura 3 referem-se a um estudoda fracao media de cloro (y) disponıvel num produto manufaturado em funcao do


tempo (x = 2(3+`), ` = 1, 2, . . . , 18, em semanas) de fabricacao do mesmo (Drapere Smith, 1981).

Tabela 2 - Fracao media de cloro (y) disponıvel num produto manufaturado comofuncao do tempo de fabricacao

0,490 0,475 0,450 0,437 0,433 0,455 0,423 0,407 0,407

0,407 0,405 0,393 0,405 0,400 0,395 0,400 0,390 0,390

10 15 20 25 30 35 40

0.38

0.40

0.42

0.44

0.46

0.48

0.50

Tempo de fabricação (semanas)

Fraç

ão d

e cl

oro

disp

onív

el

10 15 20 25 30 35 40

0.38

0.40

0.42

0.44

0.46

0.48

0.50


Fraç

ão d

e cl

oro

disp

onív

el

Figura 3 - Grafico de dispersao das observacoes da fracao media de cloro disponıvelnum produto manufaturado como funcao do tempo de fabricacao (Tabela2) sem e com a curva ajustada.

Aos dados da Tabela 2 foi ajustado o modelo

y = α xβ + ε, ε ∼ N(0, σ2).

Os valores iniciais foram obtidos fazendo-se uma regressao linear de log(y) emrelacao a log(x), sendo α0 = exp(coef1) e β0 = coef2.

Na Tabela 3, estao as estimativas dos parametros do modelo. O valor de pmostra um efeito significativo do tempo. A representacao grafica da curva ajustadae dos valores observados e mostrada na Figura 3.

Apos o ajuste do modelo µi(β) = α xβi , foram calculados os resıduos ordinarios

ri = yi − µi(β), i = 1, . . . , 18, a matriz de projecao estimada H =X(

X

T X)−1

XT

,

sendo que X ={∂µi/∂α ∂µi/∂β

}=

{xβ

i αxβi log(xi)

}=

{µi/α µi log(xi)

},

os vetores x11 ={∂2µi/∂α2

}= 0, x12 = x21 =

{∂2µi/∂α∂β

}=

{xβ

i log(xi)}

={µi log(xi)/α

}e x22 =

{∂2µi/∂β2

}=

{αxβ

i log(xi)2}

={µi log(xi)

2}

e a matrizW∗ = (I − H)(x11 x12 x22).


Tabela 3 - Estimativas do modelo para a fracao media de cloro disponıvel numproduto manufaturado como funcao do tempo de fabricacao

Parametros Coeficientes Erro Valor pEstimados Padrao

α 0,63968 0,02045 8,89e-16

β -0,13569 0,01035 5,63e-10

Note que as duas primeiras colunas de W∗ pertencem a C(X). Logo,

(I− H)x11 = (I− H)x12 = 0 e o vetor H1r correspondera aos valores ajustados da

regressao linear de r sobre (I− H)x22, em que x22 e a terceira coluna da matriz W.

O vetor de resıduos projetados e definido por r − H1r. Como q = posto(H1) = 1,a perda de informacao ao passar do subespaco dos resıduos ordinarios para osubespaco dos resıduos projetados sera pequena, pois q e bem menor quandocomparado a n − p = 16. Os resıduos projetados estudentizados sao obtidosdiretamente de (9).

A Figura 4 exibe os graficos dos resıduos ordinarios ri, dos resıduosestudentizados ti e dos resıduos projetados si versus o tempo de fabricacao xi,i = 1, . . . , n. Particularmente, a observacao # 06 destaca-se como aberrante, o queparece estar em concordancia com o posicionamento desse ponto na Figura 3. Oprograma em R encontra-se no Apendice 1.

5.2 T ecnicas de diagnostico para avaliar a interacao entre duas drogas

Os dados da Tabela 4 referem-se a um estudo sobre misturas de insulina padraoe insulina modificada, Suberoyl A1-B29 (Darby e Ellis, 1976). A resposta medidafoi a conversao de glicose (3-3H ) para um certo tipo de lipıdio em celulas de ratos.

O interesse e testar a hipotese nula que o efeito das drogas e aditivo, isto e,as drogas tem acao similar versus a hipotese alternativa que a mistura das drogaspode ter uma acao melhor (sinergismo) ou pior (antagonismo). Sob a hipotese nula,considera-se que uma dose de x1 unidades de A e x2 unidades de B e equivalentea uma dose de x1 + ρ x2 de A. Dessa forma, um modelo parcialmente nao-linearproposto para analisar esse tipo de observacoes, e especificado como

y = µ + ε = α + δ log(x1 + ρ x2) + ε, Modelo 1

em que x1 e x2 rep resen ta m, resp ec tiv a men te, a s d o ses d a s d ua s d ro g a s A e B, δ

e a in c lin a c a o c o mum d a rela c a o lo g -d o se-resp o sta , ρ e a p o ten c ia d a d ro g a B emrela c a o a d ro g a A e ε ∼ N (0 , σ 2). O mo d elo a ltern a tiv o e ex p resso p o r

y = µ + ε = α + δ lo g [x1 + ρx2 + κ(ρx1x2)1/2] + ε, Modelo 2

Rev. B ra s. B io m ., S a o P a u lo , v .2 7 , n .3 , p .3 6 0 -3 9 3 , 2 0 0 9 3 7 9

10 15 20 25 30 35 40

−0.0

10.

000.

010.

02


Res

íduo

ord

inár

io

6

10 15 20 25 30 35 40

−2−1

01

2

Tempo de fabricação (semanas)R

esíd

uo e

stud

entiz

ado

6

10 15 20 25 30 35 40

−10

12

3


Res

íduo

pro

jeta

do e

stud

entiz

ado

6

Figura 4 - Grafico dos resıduos ordinarios, estudentizados ti e dos resıduosprojetados si versus os valores de xi para as ob servacoes de fracao mediade cloro disponıvel num produto manufaturado como funcao do tempode fab ricacao (T ab ela 2 ).


Tabela 4 - Resultados obtidos para um ensaio de mistura de insulina padrao einsulina modificada, S uberoy l A 1 -B 29

M istura Insulina P adrao: D ose total Respostas para 4 repeticoes

S uberoy l A 1 -B 29 (pmol l−1) C onversao de glicose(3 -3 H )

1 1 :0 20,9 1 4,0 1 4,4 1 4,3 1 5 ,2

41 ,9 24,6 22,4 22,4 26 ,7

2 1 :1 ,8 5 5 2,9 1 1 ,7 1 5 ,0 1 2,9 8 ,3

1 06 ,0 20,6 1 8 ,0 1 9 ,6 20,5

3 1 :5 ,5 6 1 01 ,0 1 0,6 1 3 ,9 1 1 ,5 1 5 ,5

202,0 23 ,4 1 9 ,6 20,0 1 7 ,8

4 1 :1 6 ,7 1 8 1 ,0 1 3 ,8 1 2,6 1 2,3 1 4,0

3 6 2,0 1 5 ,8 1 7 ,4 1 8 ,0 1 7 ,0

5 1 :5 0,0 26 1 ,0 8 ,5 9 ,0 1 3 ,4 1 3 ,5

5 22,0 20,6 1 7 ,5 1 7 ,9 1 6 ,8

6 1 :1 5 0,0 3 09 ,0 1 2,7 9 ,5 1 2,1 8 ,9

6 1 7 ,0 1 8 ,6 20,0 1 9 ,0 21 ,1

7 0:1 3 40,0 1 2,3 1 5 ,0 1 0,1 8 ,8

6 8 1 ,0 20,9 1 7 ,1 1 7 ,2 1 7 ,4

em que κ representa a interacao entre as drogas. S e κ = 0 significa que h aacao similar entre as duas drogas, κ > 0 representa sinergismo e κ < 0 significaantagonismo.

P ara encontrar os valores iniciais foram fixados ρ = 0, 05 e κ = 0, 5 , com osdemais valores iniciais sendo iguais as estimativas de mınimos quadrados do modelo

y = α + δt + ε,

em que t = log[x1 + ρ0x2 + κ0(ρ0x1x2)1/2]. A Tabela 5 apresenta a analise de

variancia para os dois modelos. Nota-se que existe evidencia significativa para arejeicao da h ipotese H : κ = 0. Na Tabela 6 , estao as estimativas dos parametrosdo modelo com interacao.

C omo κ < 0, constata-se uma indicacao de antagonismo entre as duas drogas,isto e, que a mistura de ambas produz um efeito menor do que a soma dos efeitosindividuais das duas drogas. O valor de ρ = 0, 046 indica que a insulina padrao e,aproximadamente, 22 vezes mais eficaz do que a insulina na forma suberoy l A 1 -B 29 .

U tilizando tres tecnicas diferentes de diagnostico, apresentadas na S ecao 3 ,sera mostrado que a significancia de H : κ = 0 deve-se, essencialmente, a algumas“ misturas” extremas que incluem apenas uma das drogas. A primeira tecnica,


Tabela 5 - Analise de variancia

Modelo G.L . do Soma de quad. G.L . Soma V alor F V alor p

resıduo do resıduo de quad.

1 53 243,971

2 52 220,182 1 23,789 5,6182 0,02152

Tabela 6 - E stimativas do modelo para avaliar a interacao entre duas drogas

Parametros Coeficientes E rro V alor p

E stimados Padrao

α -17,342113 2,737100 5,65e-08

δ 10,540148 0,792484 < 2e − 16

ρ 0,046492 0,003263 < 2e − 16

κ -0,300843 0,120408 0,0157

consiste na analise do grafico dos resıduos ordinarios estudentizados t′is versus

a ordem das observacoes, que e exibido pela Figura 5. Pode-se notar nessegrafico, duas observacoes (8 e 12) mal ajustadas, as quais poderao ser consideradasaberrantes se a distribuicao dos resıduos for, aproximadamente, normal. A primeiradessas observacoes refere-se a uma “mistura” extrema.

Pode-se verificar pontos infl uentes atraves do grafico da distancia de Cook ,apresentado na Figura 5. Observou-se que quatro pontos (8, 12, 50 e 52) apresentamvalores muito maiores do que os demais, sendo um indicativo de pontos infl uentes.

A segunda tecnica e uma analise da variacao de primeiro passo da estimativaκ, quando uma observacao e excluıda. O grafico de ∆ 1

i κ versus a ordem dasobservacoes e apresentado na Figura 6. Nota-se que tres observacoes 8, 50 e 52sao mais infl uentes, sendo que cada uma individualmente, quando excluıda, causauma variacao de, aproximadamente, 15% na estimativa de κ, confirmando o quefoi mostrado, analisando-se as distancias de Cook . E ssas tres observacoes sao de“misturas” extremas.

Na Figura 6, tem-se o grafico de y − µ versus (I − H)xκ (grafico da variaveladicionada, terceira tecnica) para avaliar a infl uencia conjunta das observacoes naestimativa κ, em que xκ e um vetor n × 1 com as derivadas de µ(β , κ) em relacaoa κ, computado sob a hipotese nula H : κ = 0.

Observando a Figura 6, nota-se que tres observacoes (8, 50 e 5) tem umainfl uencia desproporcional na inclinacao da reta com o coeficiente κ. De fato, aretirada desses pontos conduz a estimativa κ = −0, 160 (D P (κ) = 0, 127), quesignifica um aumento de, aproximadamente, 47% apontando para uma interacao


0 10 20 30 40 50

−2−1

01

2

Índice

Res

íduo

est

uden

tizad

o8

12

0 10 20 30 40 50

0.00

0.05

0.10

0.15

Índice

Dis

tânc

ia d

e C

ook

8

12

50

52

Figura 5 - Grafico dos resıduos estudentizados e das distancias de Cook versus aordem das observacoes do modelo para avaliar a interacao entre duasdrogas.

0 10 20 30 40 50

−0.4

−0.2

0.0

0.2

Índice

∆ iκ

8

50

52

−3 −2 −1 0 1 2

−4−2

02

4

(I − H)x~κ

Res

íduo

ord

inár

io −

mod

elo

redu

zido

5

8

50

Figura 6 - Grafico da variacao em κ versus a ordem das observacoes e da variaveladicionada correspondente a estimativa κ do modelo para avaliar ainteracao entre duas drogas.


nula entre as duas drogas. E provavel que esse tipo de influencia seja devida adoses extremas, exageradas. Provavelmente, a relacao linear entre o log(dose) e aresposta, que e considerada quando as drogas sao analisadas separadamente, naodeve continuar valendo nesses casos. O programa em R encontra-se no Apendice 2.

5.3 Reparametrizacao do modelo MMF usando as medidas de Bates e

W atts e de Box

Considere o modelo MMF expresso por

µ = (βγ + αxδ)/(γ + xδ),

e o conjunto de dados da Tabela 7 (Ratkow sky, 1983, pag. 88), em que y representaa producao e x o tempo de cultivo de pastagens.

Tabela 7 - Dados de producao (y) de pastagens e tempos de cultivo (x)

x 9 14 21 28 42 57 63 70 79y 8,93 10,80 18,59 22,33 39,35 56,11 61,73 64,62 67,08

Na convergencia, obtem-se as estimativas α = 80, 96, β = 8, 895, γ = 49.577e δ = 2, 828. As medidas de nao-linearidade sao estimadas por γI N = 0, 180 eγP E = 90, 970 e o valor crıtico para um nıvel de significancia de 0, 05 e igual a1/2

√F = 0, 219, sendo F o quantil de 0, 95 de uma distribuicao F com 4 e 5 graus

de liberdade. Portanto, a nao-linearidade devida a parametrizacao do modelo ealtamente significativa, enquanto a nao-linearidade intrınseca e nao-significativa.

Para determinar quais os parametros que, possivelmente, estao causando essanao-linearidade acentuada, utiliza-se a medida de Box para o vies. E usual expressaro vies como uma porcentagem da estimativa correspondente. Para as estimativas α,β, γ e δ, essas porcentagens igualam, respectivamente, 1, 525%, −1, 643%, 119, 478%e 0, 921%. Nota-se, portanto, um valor muito grande para o vies de γ, indicandoque, possivelmente, uma reparametrizacao de γ possa reduzir o vies.

Ratkow sky (1983) sugere a simulacao das distribuicoes das estimativas comvies acentuado, para se ter uma ideia da reparametrizacao a ser aplicada. Nesseexemplo, a transformacao φ = g(γ) = log(γ) e a mais recomendada, obtendo-se

φ = log γ = 10, 81 e bφ = 0, 1087. Logo, a porcentagem do vies agora e igual

a 0, 1087 × 100/10, 81 ∼= 1%, uma reducao substancial em relacao a porcentageminicial. O programa em R encontra-se no Apendice 3.


Agradecimentos

Agradecemos ao apoio financeiro do CNPq e do Projeto PROCAD da CAPES.

CORDEIRO, G. M.; PRUDENTE, A. A.; DEMETRIO, C. G. B. A review of normalnonlinear models. Rev. Bras. Biom., Sao Paulo, v.27, n.3, p.360-393, 2009.

ABSTRACT: The normal nonlinear models continue to have a special treatment in

theoretical and applied statistics even after the introduction of the g eneralized linear

models. The normal nonlinear models are applied to several fi elds such as econometrics,

ag riculture, pharmacolog y , b iolog y , eng ineering , education, chemistry , etc. Their main

properties are derived from theoretical assumptions and their parameters are in g eneral

interpretab le. U nlik e linear models, the adeq uacy of the fi tting of the normal nonlinear

models are evaluated not only b y means of the diag nostic techniq ues, b ut also b y

ex tending the nonlinear b ehavior. M odels w ith nonlinear b ehavior far aw ay from the

linear structure may lead to asy mptotic results w hich do not hold in situations w here

small samples are used. In this article w e review the estimation and the main diag nostic

techniq ues and measures of nonlinearity . Applications are performed to the analy sis of

real data sets.

K E Y W O RD S: Bias correction; diag nostic analy sis; measures of curvature; normal

nonlinear models; projected residual.

R eferencias

BARTLETT, M. S. Properties of suffi ciency and statistical tests. P roc. R. S tat. S oc.

S er. A , London, v.160, n.1, p.268-282, 1937.

BATES, D. M.; W ATTS, D. G. Relative curvature measures of nonlinearity. J . R.

S tat. S oc. S er. B, London, v.42, p.1-25, 1980.

BEALE, E. M. L. Confidence region in nonlinear estimation. J . R. S tat. S oc. S er.

B, London, v.22, n.1, p.41-76, 1960.

BOX , M. J . Bias in non-linear estimation (with discussion). J . R. S tat. S oc. S er.

B, London, v.33, n.2, p.171-201, 1971.

BUSINGER, P.; GOLUB, G. H. Linear least squares solutions by householdertransformations. N um. M ath ., Heidelberg, v.7, n.3, p.269-276, 1965.

COOK , R. D.; TSAI, C. L. Residual in nonlinear regression. Biometrika, London,v.72, p.23-29, 1985.

COOK , R. D.; TSAI, C. L.; W ei, B. C. Bias in non-linear regression. Biometrika,London, v.73, n.3, p.615-623, 1986.

COOK , R. D.; W EISBERG, S. Resid uals an d in fl uen ce in regression . London:Chapman and Hall, 1982.


CORDEIRO, G. M. Improved likelihood ratio statistics for generalized linearmodels. J. R. Stat. Soc. Ser. B, London, v.45, n.2, p.401-413, 1983.

CORDEIRO, G. M. On the corrections to the likelihood ratio statistics. Biometrika,London, v.74, p.265-274, 1987.

CORDEIRO, G. M.; LIMA NETO, E. A. Modelos parametricos. In: SIMPOSIONACIONAL DE PROBABILIDADE E ESTATISTICA, 16., 2004, Caxambu. Anais

... Caxambu: ABE, 2004.

CORDEIRO, G. M.; PAULA, G. A. Improved likelihood ratio statistics forexponential family nonlinear models. Biometrika, London, v.76, p.93-100, 1989a.

CORDEIRO, G. M.; PAULA, G. A. Modelos de regressao para analise de dadosunivariados. In: COLOQ UIO BRASILEIRO DE MATEMATICA, 17., 1989, Riode Janeiro. Resumos ... Rio de Janeiro: IMPA, 1989b.

COX, D. R.; SNELL, E. J. A general definition of residuals (with discussion). J. R.

Stat. Soc. Ser. B, London, v.30, n.2, p.248-275, 1968.

DARBY , S. C.; ELLIS, M. J. A test for synergism between two drugs. Ap p l. Stat.,London, v.25, n.3, p.296-299, 1976.

DAVIDIAN, M.; GILTINAN, D.M. Nonlinear Models for repeated measurementdata: an overview and update. J. Agric., Biol. E nviron. Stat., Alexandria, v.8, n.4,p.387-419, 2003.

DRAPER, N.; SMITH, H. Ap p lied regression analy sis. New Y ork: John Wiley, 1981.709p.

EAVES, D.M. On Bayesian nonlinear regression with an enzyme example.Biometrika, London, v.70, n.2, p.373-379, 1983.

GELMAN, A.; CARLIN, J.B., STERN, H.S.; RUBIN, D.B. Bayesian data analy sis.London: Chapman and Hall/ CRC, 2004. 698p.

GILTINAN, D.; CAPIZ Z I, T.; MALANI, H. Diagnostic tests for similar action oftwo compunds. Ap p l. Stat., London, v.37, p.39-50, 1988.

GUTTMAN, I.; MEETER, D. On Beale’s measures of nonlinearity. T echnometrics,Washington, v.7, n.4, p.623-637, 1965.

JENNRICH, R. I. Asymptotic properties of nonlinear least-squares estimation. Ann.

Math. Stat., Ann Harbor, v.20, p.633-643, 1969.

JOHANSEN, S. Some topics in regression [with discussion and replay]. Scand. J.

Stat., Oxford, v.10, n.3, p.161-194, 1983.

LAWLEY , D. N. A general method for approximating to the distribution oflikelihood ratio criteria. Biometrika, London, v.43, n.3-4, p.295-303, 1956.

NELDER, J. A.; WEDDERBURN, R. W. M. Generalized linear models. J. R. Stat.

Soc. Ser. B, London, v.135, n.3, p.370– 384, 1972.

PINHEIRO, J.C.; BATES, D.M. Mixed-eff ects models in S and S-Plus. New Y ork:Springer, 2000. 528p.


RATKOWSKY, D. A. Nonlinear regression modelling. New York: Marcel Dekker,1983.

RATKOWSKY, D. A. H andbook of nonlinear regression models. New York: MarcelDekker, 1990. 241p.

RICHARDS, F. J. A. Flexible growth function for empirical use. J. Exper. Bot.,v.10, n.29, p.290-300, 1959.

SANDLAND, R. L.; McGILCHRIST, C. A. Stochastic growth curve analysis.Biometrics, London, v.35, n.1, p.255-271, 1979.

SOUSA, D. Algumas consideracoes sobre regressao nao-linear. Dissertacao(Mestrado e Estatıstica) - Instituto de Matematica e Estatıstica - USP, Sao Paulo,1986.

SOUZA, G. S. Introducao aos modelos de regressao linear e nao-ninear. Brasılia:EMBRAPA, 1998.

Recebido em 10.03.2009.

Aprovado apos revisao em 03.10.2009.


Apendice 1

## Lendo os dados

cloro= read.table(file="C:/dados_cloro.txt",header=T)

attach(cloro)

## Ajustando um modelo de regress~ao linear para obter os valores iniciais

fit.lm<-lm(log(y)~log(x))

alpha0<-fit.lm$coef[1]

beta0<-fit.lm$coef[2]

## Ajustando um modelo de regress~ao n~ao linear

fit.nls<- nls(y ~ alpha*(x^beta), start=list(alpha=exp(alpha0), beta=beta0))

summary(fit.nls)

## Grafico dos valores observados

plot(x, y, ylim=c(0.38,0.50), pch=8, xlab="Tempo de fabricac~ao (semanas)",

ylab="Frac~ao de cloro disponıvel")

## Grafico dos valores observados e ajustados

plot(x, y, ylim=c(0.38,0.50), pch=8, xlab="Tempo de fabricac~ao (semanas)",

ylab="Frac~ao de cloro disponıvel")

lines(x, fitted(fit.nls), lty=1)

## Parametros

alpha=coef(fit.nls)[1]

beta=coef(fit.nls)[2]

p<-length(coef(fit.nls))

n<-length(y)

Beta=matrix(coef(fit.nls),p)

## Primeira e segunda derivada de mu em relac~ao ao vetor de parametros Beta

g<- expression(alpha*(x^beta)) # inserir o modelo

g<-deriv(g, c("alpha","beta"), hessian = TRUE) # especificar os parametros para derivar

X<-attr(eval(g),"gradient") # X til

Xt=t(X)

XtX<- crossprod(X)

W<-attr(eval(g),"hessian") # p matrizes n x p

Wn<-matrix(W,n*p,p) # matriz W com dimens~ao (n*p) x p

## Matriz de projec~ao

H<-X%*%solve(XtX)%*%Xt

h<-diag(H)

r<-summary(fit.nls)$resid ## Resıduo ordinario

s<-summary(fit.nls)$sigma ## Desvio padr~ao

t<-r/(s*sqrt(1-h)) ## Resıduo estudentizado

## Grafico dos resıduos ordinarios

plot(x, r, pch=8, xlab="Tempo de fabricac~ao (semanas)",

ylab="Resıduo ordinario")

abline(0,0,lty=1, col="black")

identify(x,r)

## Grafico dos resıduos estudentizados

plot(x, t, pch=8, ylim=c(-2.5,2.5), xlab="Tempo de fabricac~ao (semanas)",


ylab="Resıduo estudentizado")

abline(-2,0,lty=2, col="black")


identify(x,t)

## Resıduo projetado

I=diag(n)

xij<-matrix(0,n,1)

T<-matrix(0,n,p*p)

for (j in 1:(p*p))

{

l=(1+((j-1)*n)):(n+((j-1)*n))

xij<-W[l]

T[,j]= as.vector(xij)

}

## Vetor xtil_{22}

Tf=cbind(T[,4])

S<-(I-H)%*%Tf

St=t(S)

StS<- crossprod(S)

H1<-S%*%solve(StS)%*%St

V<-cbind(X,S)

Vt=t(V)

VtV<- crossprod(V)

H2<-V%*%solve(VtV)%*%Vt

trIH2<-sum(diag(I-H2))

sigma2=(t(r)%*%(I-H2)%*%r)/trIH2

sigma=sqrt(sigma2)

rp<-((I-H2)%*%r)/(sigma*diag(sqrt(I-H2)))

## Grafico dos resıduos projetado estudentizados

plot(x, rp, pch=8, xlab="Tempo de fabricac~ao (semanas)",

ylab="Resıduo projetado estudentizado")

identify(x,rp)


Apendice 2

## Lendo os dados

drogas= read.table(file="C:/dados_drogas.txt",header=T)

attach(drogas)

## Ajustando um modelo de regress~ao linear para obter os valores iniciais

rho0<-0.05

kappa0<-0.5

t<-log(x1+rho0*x2+kappa0*(rho0*x1*x2)^0.5)

fit.lm<-lm(y ~ t)

alpha0<-fit.lm$coef[1]

delta0<-fit.lm$coef[2]


fit.nls<- nls(y ~ alpha+delta*log(x1+rho*x2+kappa*(rho*x1*x2)^0.5),

start=list(alpha=alpha0,delta=delta0,rho=rho0,kappa=kappa0))

summary(fit.nls)

## Ajustando um modelo de regress~ao n~ao linear com kappa=0 - modelo reduzido

fit.nls1<- nls(y ~ alpha+delta*log(x1+rho*x2),

start=list(alpha=alpha0,delta=delta0,rho=rho0))

summary(fit.nls1)

## Analise de Variancia

anova(fit.nls1, fit.nls)

## Parametros

alpha<-coef(fit.nls)[1]

delta<-coef(fit.nls)[2]

rho<-coef(fit.nls)[3]

kappa<-coef(fit.nls)[4]


n<-length(y)

dalpha<-rep(1,n)

ddelta<-log(x1+rho*x2+kappa(rho*x1*x2)^0.5)

drho<-delta*(x2+0.5*kappa*((x1*x2)^0.5)*(rho^(-0.5)))/(x1+rho*x2+kappa*(rho*x1*x2)^0.5)

dkappa<-delta*((rho*x1*x2)^0.5)/(x1+rho*x2+kappa*(rho*x1*x2)^0.5)

X<-matrix(c(dalpha, ddelta, drho, dkappa), n)

Xt=t(X)

XtX<- crossprod(X)

H<-X%*%solve(XtX)%*%Xt # Matriz de projec~ao

h<-diag(H)

r<-summary(fit.nls)$resid # Resıduo ordinario

s<-summary(fit.nls)$sigma # Desvio padr~ao

t<-r/(s*sqrt(1-h)) # Resıduo estudentizado

D<-(t^2/p)*(h/(1-h)) # Distancia de Cook

## Grafico dos resıduos estudentizados

plot(t, pch=8, xlab="Indice", ylab="Resıduo estudentizado")

abline(-2,0,lty=2, col="black")


identify(t)


## Grafico da distancia de Cook

plot(D, pch=8, xlab="Indice", ylab="Distancia de Cook")

identify(D)

## Variac~ao de primeiro passo da estimativa Kappa

beta<-as.vector(coef(fit.nls))

beta1i<-matrix(0,n,p)

for(i in 1:n)

{

beta1i[i,]<-beta-((solve(XtX)%*%(X[i,]*r[i]))/(1-h[i]))

}

DP_kappa<-coef(summary(fit.nls))[4,2]

k<-as.vector(kappa-beta1i[,4])/DP_kappa

# Grafico Variac~ao em Kappa

plot(k, pch=8, xlab="Indice", ylab=expression(Delta[i]*hat(kappa)))

identify(k)

## Modelo Reduzido

alpha.1<-coef(fit.nls1)[1]

delta.1<-coef(fit.nls1)[2]

rho.1<-coef(fit.nls1)[3]

p.1<-length(coef(fit.nls1))

n<-length(y)

dalpha.1<-rep(1,n)

ddelta.1<-log(x1+rho.1*x2)

drho.1<-delta.1*(x2)/(x1+rho.1*x2)

X.1<-matrix(c(dalpha.1, ddelta.1, drho.1), n)

Xt.1=t(X.1)

XtX.1<- crossprod(X.1)

H.1<-X.1%*%solve(XtX.1)%*%Xt.1 # Matriz de projec~ao

h.1<-diag(H.1)

r.1<-summary(fit.nls1)$resid # Resıduo ordinario

# Xkappa modelo completo sob a hipotese H: kappa=0

dkappa.0<-delta*((rho*x1*x2)^0.5)/(x1+rho*x2)

I=diag(nrow(X.1))

vec<-(I-H.1)%*%dkappa.0

# Grafico da variavel adicionada

plot(vec, r.1, pch=8, xlab=expression((I-hat(H))*tilde(x)[hat(kappa)]),

ylab="Resıduo ordinario - modelo reduzido")

identify(vec,r.1)

## Ajustando um modelo de regress~ao n~ao linear sem as observac~oes 5,8,50

fit<- nls(y ~ alpha+delta*log(x1+rho*x2+kappa*(rho*x1*x2)^0.5),

start=list(alpha=alpha0,delta=delta0,rho=rho0,kappa=kappa0), subset=-c(5,8,50))

summary(fit)


Apendice 3

## Lendo os dados

prod_past= read.table(file="C:/dados_pastagens.txt",header=T)

attach(prod_past)


fit.nls<- nls(y ~ (beta*gamma+alpha*x^delta)/(gamma+x^delta),

start=list(alpha=70, beta=9, gamma=50000, delta=3))

summary(fit.nls)

## Parametros

alpha=coef(fit.nls)[1]

beta=coef(fit.nls)[2]

gamma=coef(fit.nls)[3]

delta=coef(fit.nls)[4]


n<-length(y)

Beta=matrix(coef(fit.nls),p)

## Primeira e segunda derivada de mu em relac~ao ao vetor de parametros Beta

g<- expression((beta*gamma+alpha*x^delta)/(gamma+x^delta)) # inserir o modelo

g<-deriv(g, c("alpha","beta","gamma","delta"), hessian = TRUE) # especificar os parametros para derivar

X<-attr(eval(g),"gradient") # X til

Xt=t(X)

XtX<- crossprod(X)

W<-attr(eval(g),"hessian") # p matrizes n x p

Wn<-matrix(W,n*p,p) # matriz W com dimens~ao (n*p) x p

I=diag(n)

## Matriz de projec~ao

H<-X%*%solve(XtX)%*%Xt

h<-diag(H)

r<-summary(fit.nls)$resid ## Resıduo ordinario

s<-summary(fit.nls)$sigma ## Desvio padr~ao

t<-r/(s*sqrt(1-h)) ## Resıduo estudentizado

## Vies de ordem n^{-1} de Beta

wi<-matrix(0,p,p)

d=c(0)

w=matrix(0,n,p^2)

for (i in 1:n)

{

wi<-rbind(Wn[i,],Wn[n+i,],Wn[(2*n)+i,],Wn[(3*n)+i,]) # Altera de acordo com a qtd de parametros

w[i,]= as.vector(wi) # Armazena wi como um vetor na i-esima linha de w

d[i]=-0.5*s^2*(sum(diag(solve(XtX)%*%wi)))

}

b=solve(XtX)%*%Xt%*%d

abs(b)/abs(Beta) # Se maior que 0,2 a estimativa possui vies apreciavel

(b*100)/Beta # Vies como uma porcentagem da estimativa


## Medidas de n~ao-linearidade

rho<-s*sqrt(p)

Xrho=X/rho

qrstr <- qr(Xrho) #decomposic~ao QR de Xrho

Q <- qr.Q(qrstr,complete=TRUE)

R <- qr.R(qrstr,complete=TRUE)

Rtil <- qr.R(qrstr)

L=solve(Rtil)

ui=matrix(0,p,p)

u=matrix(0,n,p^2)

for (i in 1:n)

{

ui=t(L)%*%((matrix(w[i,],p,p))/rho)%*%L

u[i,]= as.vector(ui) # Armazena ui como um vetor na i-esima linha de u

}

a<-t(Q)%*%u

# Valor crıtico, sendo F o quantil de 0,95 de uma distribuic~ao F

F=qf(0.95, 1, 1)

vc=1/(2*sqrt(F))


Documents

Gauss Moutinho CORDEIRO A ndrea A ndrade P RU DEN T E ...jaguar.fcav.unesp.br/RME/fasciculos/v27/v27_n3/A4_Andreia.pdf · exibem varia»c~ao heterog^enea e muitas vezes do tipo heterosced¶astica