Archive

Posts Tagged ‘lm’

Interpretação da matriz de covariância das estimativas

Quero falar de um erro comum de interpretação de resultados em análise de regressão através de um exemplo. Considere que você tem pessoas com diferentes pesos e diferentes alturas. Facilmente você aceita que uma pessoa mais alta tem maior peso, certo? Ou seja, existe uma correlação positiva entre peso e altura. Pois bem, vamos simular observações desse experimento e em seguida ajustar uma regressão linear simples para relação peso-altura.

require(MASS)

set.seed(12345)
da <- as.data.frame(mvrnorm(80, c(altura=30, peso=60),
                            matrix(c(2,1.9,1.9,3), 2,2)))
plot(peso~altura, da)

m0 <- lm(peso~altura, data=da) # ajusta a reta
summary(m0)              # estimativa dos parâmetros
abline(m0)               # adiciona uma reta ao gráfico
vcov(m0)                 # covariância das estimativas

Aqui está o ponto que eu quero comentar: a interpretação da matriz de covariância. Perceba que a covariância foi negativa. Já vi gente interpretando isso da seguinte forma: espera-se que pessoas mais altas tenham menor peso. Duplamente errado. Primeiro nos sabemos por experiência que a correlação de peso e altura é positiva. Segundo, a matriz de covariância se refere as estimativas dos parâmetros e não as variáveis envolvidas. Nunca a interprete dessa forma.

Então como interpretar? Bem, a matriz de covariância das estimativas, é um reflexo da função objetivo ao redor da solução. A função objetivo nesse caso é minimizar a soma de quadrados (SQ). Então se eu aumento o valor do parâmetro b0, o parâmetro b1 tem que diminuir para compensar, ou seja, para diminuir a SQ. Ocorre um efeito compensatório aqui. Além do mais, esse efeito pode ser eliminado facilmente aplicando uma reparametrização, como por exemplo, centrar a covariável na média.

da$alturac <- with(da, altura-mean(altura)) # centrando a covariável
par(mfrow=c(1,2))
plot(peso~altura, da)
abline(m0)
plot(peso~alturac, da)

m1 <- lm(peso~alturac, data=da) # ajusta a reta
summary(m1)              # estimativa dos parâmetros
abline(m1)               # adiciona uma reta ao gráfico
vcov(m1)                 # covariância das estimativas

require(ellipse)
par(mfrow=c(1,2)) # regiões de confiança conjunta
plot(ellipse(vcov(m0)), type="l")
plot(ellipse(vcov(m1)), type="l")

cov2cor(vcov(m0))
cov2cor(vcov(m1))

Perceba que os resultados são os mesmos em termos de estatísticas de teste e medidas de ajuste. Isso é esperado por eu só fiz uma traslação da altura. Mas o importante é que agora a matriz de covariância tem covariâncias praticamente nulas, que é resultado da translação. O intercepto é então o valor esperado para alturac igual a zero (centro dos dados). Veja porque essa covariância é zero: se eu aumentar b0, não adiante eu alterar b1 que não vai diminuir a SQ, e vice-versa, porque agora eles são ortogonais. Ortogonalidade entre parâmetros é uma característica desejável pois permite você inferir sobre um deles sem considerar o outro. Além disso, tem vantagens do ponto de vista de estimação por métodos numéricos de otimização. Em outras palavras, pegando conceitos de probabilidade, se a distribuição amostral de duas variáveis aleatórias é normal bivariada com covariância nula, a distribuição condicional de A|B é igual a distribuição marginal de A pela independência. Reforçando, eu não preciso conhecer valores de B para descrever A. Considerando tudo que foi argumentado, é sempre preferível que você adote o modelo que apresente menor covariância entre os parâmetros. Até a próxima ridícula.

Categorias:inferência Tags:, ,

Semântica para descrever modelos

A função lm() do R resolve 90% dos casos de análise de experimentos. Não são para lm() modelos não lineares, de efeito aleatório, de distribuição diferente da gaussiana ou dados não independentes (correlacionados). Para uso eficiente da função é importante saber declarar os modelos. Para isso existe uma semântica (termo bem apropriado usado pelo Fernando Toledo) para declarar os modelos. Aqui eu vou listar alguns dos casos possíveis.

fórmula significado
A+B efeito simples aditivo entre fatores
A:B efeito do produto cartesiano de níveis dos fatores
A*B*C efeitos simples e todas interações possíveis
(A+B+C)^2 efeitos simples e interações de até 2 termos (duplas)
(A+B+C+C+E)^3 efeitos simples e interações de até 3 termos (triplas)
A/B efeito simples de A e de B aninhado em A, A+A:B
I(x^2) operador I() torna x^2 literal (remove o poder semântico de ^2)
-A:B remove o termo do modelo
A*B*C-A:B:C declara modelo triplo mas remove interação tripla, (A+B+C)^2
0+A remove intercepto (coluna de 1) do modelo
-1+A idem ao anterior
poly(x, degree=2) polinômio ortogonal de grau 2 em x
poly(x, degree=2, raw=TRUE) polinômio cru de grau 2 em x
A*(B+C) propriedade distributiva, A*B+A*C
A/(B+C) idem ao anterior, A/B+A/C
A/B/C propriedade recursiva, A+A:B+A:B:C
. representa todas as colunas, exceto a da resposta
(.)^2 inclui todas as interações duplas das variáveis

Além dessas opções de sintaxe, ainda temos a possibilidade de usar o operador semântico Error() para declarar termos de efeito aleatório. Esse operador só é interpretado pela função aov(). Com ele fazemos a análise de experimentos em parcelas subdivididas, por exemplo, de tal forma que a análise de variância apresenta os testes F considerando os quadrados médios corretos.
Além disso podermos usar a função terms(..., keep.order=TRUE) para fixar a ordem dos termos na análise de variância, pois os termos são ordenados pela ordem de magnitude, assim, sempre interações vem depois de efeitos simples, algumas vezes não desejamos isso, como é o caso da análise de experimentos em blocos incompletos. Os exemplos vão tornar esses apontamentos mais claros. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# observe as matrizes de delineamento geradas nessas situações

da <- replicate(5, gl(3,1), simplify=FALSE)
da <- data.frame(do.call(expand.grid, da))
names(da) <- LETTERS[1:ncol(da)]
da$x <- runif(nrow(da))

formulas <- list(~A+B, ~A:B, ~A*B, ~(A+B)^2, ~(A+B+C+D+E)^3,
                 ~A/B, ~A/B/C, ~A/(B+C), ~A*(B+C), ~A*(B+C)^2,
                 ~(A+B+C)^2-A:B, ~x, ~0+A, ~-1+A, ~0+x, ~A*x, ~A:x,
                 ~A/x, ~x+I(x^2), ~poly(x, degree=2),
                 ~poly(x, degree=2, raw=TRUE), ~A+A:B+C,
                 terms(~A+A:B+C, keep.order=TRUE))
names(formulas) <- apply(sapply(formulas, paste), 2, paste, collapse="")

X <- lapply(formulas, model.matrix, data=da)
lapply(X, head)

#-----------------------------------------------------------------------------
# para obter a equação de um polinômio para os níveis

soja <- read.table("http://www.leg.ufpr.br/~walmes/data/soja.txt",
                   header=TRUE, sep="\t", dec=",")
str(soja)

soja$A <- factor(soja$agua)
soja$K <- soja$potassio

m0 <- lm(rengrao~bloco+A*(K+I(K^2)), soja) # I() torna literal
m1 <- lm(rengrao~bloco+A*poly(K, degree=2, raw=TRUE), soja)  # original
m2 <- lm(rengrao~bloco+A*poly(K, degree=2, raw=FALSE), soja) # ortogonal

cbind(coef(m0), coef(m1), coef(m2)) # 2 primeiras colunas são iguais

m3 <- lm(rengrao~0+A/(K+I(K^2))+bloco, soja) # para ter os coeficientes
coef(m3) # os interceptos, termos lineares e termos quadráticos
c0 <- coef(m3)[grep("bloco", names(coef(m3)), invert=TRUE)]
matrix(c0, ncol=3, dimnames=list(levels(soja$A), c("b0","b1","b2")))

#-----------------------------------------------------------------------------
# tranformações dentro da fórmula

m0 <- lm(rengrao~bloco+A*(K+sqrt(K)), soja)  # não requer I()
m0 <- lm(rengrao~bloco+A*I(K-mean(K)), soja) # centrar na média

#-----------------------------------------------------------------------------
# fazendo de conta que o experimento é em parcelas subdividas para ilustrar
# o uso do operador semântico Error()

soja$K <- factor(soja$K)

# combinar os níveis de bloco e A nos identificamos as parcelas
m4 <- aov(rengrao~bloco+A*K+Error(bloco:A), soja) # dá mensagem de aviso
summary(m4) # anova separa pelos estratos

soja <- transform(soja, parcela=interaction(bloco, A))

# é o mesmo modelo
m5 <- aov(rengrao~bloco+A*K+Error(parcela), soja) # não dá mensagem de aviso
summary(m5) # anova separa pelos estratos

#-----------------------------------------------------------------------------
# fazendo uso da função terms() para fixar a ordem dos termos na anova()

bib1 <- read.table("http://www.leg.ufpr.br/~walmes/data/bib1.txt",
                   header=TRUE, sep="\t")
str(bib1)

bib1 <- transform(bib1,
                  repetição=factor(repetição),
                  variedade=factor(variedade),
                  bloco=factor(bloco))
str(bib1)

# repetição + bloco dentro de repetição + variedade
m0 <- lm(y~repetição/bloco+variedade, bib1)
anova(m0) # por ser de maior ordem o termo repetição:bloco é o último

m1 <- lm(terms(y~repetição/bloco+variedade, keep.order=TRUE), bib1)
anova(m1) # ordem obedeceu a ordem dentro da fórmula

#-----------------------------------------------------------------------------
Categorias:delineamento Tags:, ,

Experimento fatorial duplo com um tratamento adicional

Produção de cultivares em função do nível de fertilizante (95% de confiança). As linhas horizontais em cinza representam o desempenho da cultivar de referência. Linhas pretas representam o desempenho das demais como função do nível de fertilizante. A seta aponta o nível de fertilizante que confere a máxima produção da cultura. Análise baseada em dados simulados.

Em algumas áreas do conhecimento é quase regra fazer experimentos incluindo uma tratamento adicional. Um experimento fatorial é aquele que possui duas ou mais fontes de variação. Em geral, têm-se presente no experimento todos os níveis possíveis das combinações entre os níveis desses fatores, ou em outras palavras, todos os níveis possíveis do produto cartesiano. Os fatores podem ser qualitativos (ou categóricos) ou quantitativos (ou métricos). Um tratamento adicional é um nível que não corresponde à estrutura fatorial, ou seja, é um nível que não veio da combinação dos níveis dos fatores. Para ficar mais claro, segue alguns exemplos:

  • Para um determinado município, sabe-se que a melhor cultivar de soja é a A com produção máxima ao receber a aplicação de 100 kg de fertilizante. Uma empresa está estudando 5 novas cultivares para o município e precisa definir a dose de fertilizante para a máxima produção. Além de definir a dose ótima, deseja-se comparar a produção com a cultivar A. Para isso faz-se um experimento com os 5 níveis de cultivar de soja (B, C, D, E, F) e um número de níveis suficientes de fertilizante que permita determinar o ponto de máximo, diga-se, 5 níveis (50, 75, 100, 125, 150 kg). O experimento terá 25 níveis (5×5) vindos da porção fatorial e um nível de referência, chamado de testemunha ou controle, que é cultivar A com 100 kg de fertilizante. Esse experimento é um fatorial 5×5+1.
  • Em um experimento com herbicidas, deseja-se definir a dose mínima de aplicação para controlar plantas daninhas. São 4 níveis de novos herbicidas (A, B, C, D) e 6 níveis de concentração (0,1; 0,5; 1; 2; 5; 10%). A variável resposta de interesse é produtividade da cultura. Para verificar se os produtos causam toxidade à cultura, tem-se um nível que é capina manual (testemunha positiva), que é o controle de plantas daninhas sem produto químico, onde se espera a maior produtividade. Para ver a eficiência desses herbicidas, no estudo se incluí um nível que é sem a capina (testemunha negativa), onde se espera a menor produtividade. Este experimento é um fatorial 4×6+2.

Perceba que pelas descrições acima ficaram evidentes às hipóteses de cada experimento. No primeiro é definir a dose ótima de cada cultivar e depois comparar o desempenho das cultivares na melhor dose com a cultivar A. No segundo é verificar se os herbicidas causam toxides à cultura (comparar com a testemunha positiva), qual a perda de produção quando não se faz o controle de plantas daninhas (comparar com a testemunha negativa) e definir a dose mínima para controlar plantas daninhas para cada herbicida.

O modelo estatístico que representa o primeiro experimento é

Y_{ijk} = M+C_i + F_j + CF_{ij} + \epsilon_{ijk} \quad \text{e} \quad Y_{k} = M + A + \epsilon_{k},

em que Y são os valores observados, M um valor inerente à todas as observações, C representa o efeito dos níveis de cultivar, F representa o efeito dos níveis de fertilizante, CF representa o efeito da combinação dos níveis de cultivar e fertilizante, A é o efeito do tratamento adicional e \epsilon é um erro aleatório.

Veja que até aqui foi tudo muito comum, com exceção ao fato de ser necessário escrever o modelo com duas expressões, uma para descrever a parte fatorial e outra para a parte adicional. No R não conseguimos declarar o modelo com duas fórmulas e é justamente sobre isso que vou tratar nesse post.

Podemos representar o modelo acima como um fatorial incompleto. Pense assim: são 6 nível de cultivar e apenas 5 deles combinados com 5 níveis de fertilizante. É assim que vamos declarar o modelo no R. Com algumas manipulações, iremos repartir a soma de quadrados da porção fatorial da porção adicional. A repartição será ortogonal que nem sempre representa as hipóteses que nós estabelecemos. Então por que fazemos essa repartição? Só porque ela é ortogonal.

Vou usar dados simulados para exemplificar o procedimento. O efeito dos níveis de fertilizante serão representados por um polinômio de segundo grau, apenas para que apresente ponto de máximo. No script abaixo eu gero os dados, faço os ajuste dos modelos, apresento os quadros de análise de variância, faço a repartição da soma de quadrados, a estimação da dose ótima de fertilizante para a resposta máxima e a comparação com a testemunha. O código foi comentado para orientar o leitor. Se você quer saber como é a análise do fatorial com dois tratamentos adicionais, deixe um comentário nesse post. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# gerando dados conforme o modelo estatístico 5x5+1, depois explico o fert=101

adi <- expand.grid(cult=gl(1,5,la=LETTERS[1]), fert=101)
fat <- expand.grid(cult=gl(5,5,la=LETTERS[2:6]), fert=seq(50,150,25))
da <- rbind(adi, fat)

theta <- c(c(-194.29, -197.26, -197.85, -203.03, -190.20, -190.45),
           c(9.1797, 8.2686, 8.6437, 9.3438, 8.8773, 8.1872),
           c(-0.03382, -0.03479, -0.03632, -0.03341, -0.03597, -0.03675))
X <- model.matrix(~-1+cult/(fert+I(fert^2)), data=da)
da$eta <- X%*%matrix(theta)
set.seed(2)
da$y <- da$eta+rnorm(nrow(da),0,30)

require(lattice)
xyplot(y~fert|cult, data=da, type=c("p","a"))

#-----------------------------------------------------------------------------
# usar o confundimento

da$Fert <- factor(da$fert) # esse passo explica o fert=101
levels(da$Fert) # o nível 101 ocorre nas mesmas observações onde
levels(da$cult) # ocorre o nível A, portanto estão confundidos

#-----------------------------------------------------------------------------
# ajuste do modelo considerando os fatores todos de níveis categórios
# poderia ir direto para a análise com polinômio mas esse passo ilustra
# o procedimento para 2 fatores de níveis categóricos

m0 <- aov(y~cult*Fert, data=da)
anova(m0) # compare os graus de liberdade com o número de níveis

# 6-1 gl para cult
# 5-1 gl para Fert
# (6-1)*(5-1) para cult:Fert? não, é (5-1)*(5-1)
# pois existe níveis de fert em apenas 5 níveis de cultivar

#-----------------------------------------------------------------------------
# os contrates entre níveis são os seguintes

contrasts(da$cult) # todo nível de cult contra o nível A
contrasts(da$Fert) # todo nível de Fert contra o nível 50

#-----------------------------------------------------------------------------
# para decompor a soma de quadrados em parte fatorial e parte adicional
# devemos mudar esses contrastes. Um deles deve ser o nível adicional vs todos
# os demais níveis daquele fator. Para isso podemos trocar a ordem dos níveis
# e o tipo de contraste para Helmert

contrasts(C(da$cult, contr=helmert)) # nível F contra os demais, trocar para A
contrasts(C(da$Fert, contr=helmert)) # nível 150 trocar para 101

da$cult <- factor(da$cult, levels=rev(levels(da$cult)))
da$Fert <- factor(da$Fert, levels=c(50,75,100,125,150,101))

levels(da$cult) # A é o último nível
levels(da$Fert) # 101 é o último nível

contrasts(C(da$cult, contr=helmert)) # A vs demais
contrasts(C(da$Fert, contr=helmert)) # 101 vs demias

m0 <- aov(y~cult*Fert, data=da, # passando o tipo de contraste
          contrasts=list(cult=contr.helmert, Fert=contr.helmert))
anova(m0) # mesmo quadro de anova

# agora vamos repartir a soma de quadrados. Note que existe uma hipótese em
# teste quando fazemos esse procedimento. Essa hipótese surge porque estamos
# fazendo um partição ortogonal específica.

summary(m0, expand=FALSE, split=list("cult"=list(fat=1:4, adi=5)))

# de 1:4 são estimativas dos contrastes entre níveis da porção fatorial
# 5 é a estimativa do contraste entre o nível adicional contra os outros
# e essa é a hipótese que a partição ortogonal testa
# mas que nem sempre é a hipótese que temos interesse
# agora é só aplicar os testes para as hipótese de interesse, sejam elas
# comparações múltiplas, contrastes planejados, etc.

#-----------------------------------------------------------------------------
# ajustando modelo polinomio de 2 grau para o efeito de fertilizante

m1 <- lm(y~-1+cult/(fert+I(fert^2)), data=da)
summary(m1)

# usamos essa fórmula para termos interpretação direta das estimativas
# perceba que para o nível A obtivemos NA, lógico, não temos como estimar
# efeito de fertilizante porque não temos níveis para essa cultivar

#-----------------------------------------------------------------------------
# estimando as doses que conferem máxima resposta em cada cultivar
# solução é derivar o modelo quadrático e igualar a zero

D(expression(b0+b1*x+b2*x^2), "x") # igualar a zero e isolar x
fert.max <- -0.5*coef(m1)[6+1:6]/coef(m1)[12+1:6]
fert.max <- fert.max[-6]
fert.max # as doses de fertilizante que conferem o máximo para cada cultivar

#-----------------------------------------------------------------------------
# vamos ajustar com a aov() que omite os NA para então usar a estimable()

m2 <- aov(y~-1+cult/(fert+I(fert^2)), data=da)
summary.lm(m2) # não contém NA

#-----------------------------------------------------------------------------
# agora é só estimar a produção máxima para cada nível de cultivar passando a
# para estimable() a matrix que identifica as funções lineares envolvidas

con <- cbind(diag(5), 0, diag(5)*fert.max, diag(5)*fert.max^2)
con <- rbind(con, 0); con[6,6] <- 1
rownames(con) <- paste(levels(da$cult),
                       c(round(fert.max,1),"100.0"), sep="-")
con

require(gmodels)
e0 <- estimable(m2, cm=con, conf.int=0.95)
e0 # quadro de estimativas da produção nas doses de máximo

#-----------------------------------------------------------------------------
# fazer um gráfico da produção máxima com os intervalos

xlim <- extendrange(e0[,6:7])
barchart(rownames(e0)~Estimate, data=e0, xlim=xlim,
         xlab="Produção", ylab="Cultivar-dose",
         panel=function(x,y,...){
           panel.barchart(x,y,...)
           panel.arrows(e0$Upper.CI, y, e0$Lower.CI, y,
                        code=3, angle=90, length=0.1)
         })

#-----------------------------------------------------------------------------
# fazendo as comparações das produções máximas com a testemunha

contr <- t(sapply(1:5, function(l) con[l,]-con[6,]))
rownames(contr) <- paste(rownames(con)[-6], "vs A")
e1 <- estimable(m2, cm=contr, conf.int=0.95)
e1

# segundo os contrastes, apenas o cultivar D é superior ao A quando recebe
# a dose 158 de fertilizante

#-----------------------------------------------------------------------------
# confeccionando o gráfico final

pred <- expand.grid(cult=levels(da$cult)[-6], fert=seq(50,170,1))
aux <- predict(m2, newdata=pred, interval="confidence")
aux <- stack(as.data.frame(aux))
pred <- cbind(pred, aux)
str(pred)

require(latticeExtra)

#png("f020.png", w=600, h=300)
xyplot(y~fert|cult, subset(da, cult!="A"), col="green4",
       xlim=extendrange(pred$fert), ylim=c(0,525), layout=c(5,1),
       xlab="Nível de fertilizante", ylab="Produtividade da cultura",
       strip=strip.custom(bg="gray90"),
       panel=function(...){
         panel.xyplot(...)
         panel.abline(h=e0[6,c(1,6,7)], lty=c(1,2,2), col="gray50")
         panel.arrows(fert.max[which.packet()], e0[which.packet(),1],
                      fert.max[which.packet()], 0, length=0.1)
       })+
  as.layer(xyplot(values~fert|cult, groups=ind, data=pred,
                  type="l", col=1, lty=c(1,2,2)))
trellis.focus("panel", 1, 1)
panel.text(50, e0[6,1], label="Cultivar A", pos=4, font=3, cex=0.8)
trellis.unfocus()
#dev.off()

#-----------------------------------------------------------------------------

Análise de resíduos em regressão não linear

Conjunto de gráficos para análise dos resíduos de um modelo de regressão não linear.

A análise dos resíduos de um modelo é feita para verificar a plausividade das pressuposições envolvidas. Os modelos linerares de regressão clássico, ou seja, aqueles em que as observações são realizações independentes (independência) e apresentam a mesma dispersão (homogeneidade de variância), podem ser ajustados no R com a função lm() e aov(). Para objetos dessas classes existe um método plot.lm() que apresenta os gráficos de análise de resíduo. Porém, modelos de regressão não linear podem ser ajustados com a função nls() que não possui um método para análise de resíduos.

Para aplicação de inferência (testes de hipótese, intervalos de confiança, etc), o modelo não linear é aproximado linearmente. Com isso eu quero dizer que podemos usar as mesmas técnicas de análise de resíduos para modelos de regressão não linear. A matriz gradiente do modelo (de derivadas da função em relação ao vetor de parâmetros) pode ser usada dentro da função lm() e com isso podemos obter facilmente os gráficos de análise de resíduos.

Esse procedimento ainda permite obter a análise de variância para o modelo de regressão não linear, fazer a partição ortogonal da soma de quadrados total em devido ao modelo de regressão e devido aos desvios de regressão. No entanto, esse não deveria ser o quadro apresentado se compararmos o que obtemos com modelo de regressão linear. No modelo linear particionamos a soma de quadrados total corrigida para o intercepto. O nosso modelo de regressão não linear se torna um modelo de intercepto se V e D forem zero (ver modelo do CMR). Portanto, o modelo nulo seria aquele em que estimamos apenas A e que podemos comparar com o modelo completo por meio da função/método anova.nls(). Esse é o teste de hipótese para o modelo de regressão não linear. Vale lembrar que alguns modelos não lineares não podem ser reduzidos a modelos de intercepto via restrição paramétrica. Para esses modelos deve-se ter cuidados ao usar o R², inclusive.

Abaixo apresento o procedimento para ajuste de um modelo não linear, extração da matriz gradiente e obtenção dos gráficos de resíduos. No final eu fiz uma função (R2()) para facilitar a obtenção do quadro de análise de variância e R². Não é porque você tem uma função para isso que vai sair aplicando ela em qualquer modelo não linear, ok? Até a próxima ridícula.

#-----------------------------------------------------------------------------
# dados de postássio liberado em função do tempo

klib <- data.frame(k=c(51.03, 57.76, 26.60, 60.65, 87.07, 64.67,
                     91.28, 105.22, 72.74, 81.88, 97.62, 90.14,
                     89.88, 113.22, 90.91, 115.39, 112.63, 87.51,
                     104.69, 120.58, 114.32, 130.07, 117.65, 111.69,
                     128.54, 126.88, 127.00, 134.17, 149.66, 118.25,
                     132.67, 154.48, 129.11, 151.83, 147.66, 127.30),
                   t=rep(c(15, 30, 45, 60, 75, 90,
                     120, 150, 180, 210, 240, 270), each=3))

#-----------------------------------------------------------------------------
# ajustando o modelo de regressão não linear aos dados

n0 <- nls(k~A*t/(V+t)+D*t, data=klib, start=list(A=90, V=15, D=0.21))
summary(n0)

#-----------------------------------------------------------------------------
# extraindo a matriz gradiente avaliada nas estimativas dos parâmetros

F <- attr(n0$m$fitted(), "gradient")
F

#-----------------------------------------------------------------------------
# passando a matriz gradiente para a lm(), importante remover intercepto

m0 <- lm(k~-1+F, data=klib)

#-----------------------------------------------------------------------------
# gráfico de análise dos resíduos

#png("f007.png", w=500, h=500);
par(mfrow=c(2,2), mar=c(5.1,4.1,4.1,2.1))
plot(m0)
mtext("Análise de resíduos para modelo de regressão não linear",
      outer=TRUE, line=-2, cex=1.4) 
layout(1)
#dev.off()

#-----------------------------------------------------------------------------
# veja que o ajuste é o mesmo pelas medidas abaixo

cbind(fitted(m0), fitted(n0))       # valores ajustados
cbind(residuals(m0), residuals(n0)) # valores preditos
c(summary(m0)$sig, summary(n0)$sig) # estimativa do desvio padrão residual
vcov(m0); vcov(n0)                  # matriz de covariância das estimativas

#-----------------------------------------------------------------------------
# quadro de anova com SQ de regressão e SQ de resíduo

anova(m0) # partição da soma de quadrados total
anova(n0, lm(k~1, klib)) # SQ do modelo não linear corrigido para intercepto

#-----------------------------------------------------------------------------
# função que retorna a anova e R2 para modelos de regressão não linear

R2 <- function(nls.obj){
  da <- eval(nls.obj$data)
  resp.name <- all.vars(summary(nls.obj)$formula)[1]
  form <- paste(resp.name, "~1", sep="")
  m0 <- lm(form, da)
  an <- anova(nls.obj, m0)
  sqn <- deviance(nls.obj)
  sqe <- deviance(m0)
  r2 <- 1-(sqn/sqe)
  aov <- data.frame(fv=c("regression","residuals"),
                    gl=c(-an$Df[2],an$Res.Df[1]),
                    sq=c(-an$Sum[2],an$Res.Sum[1]))
  aov$qm <- aov$sq/aov$gl
  aov$F <- c(aov$qm[1]/aov$qm[2], NA)
  aov$"Pr(>F)" <- c(1-pf(aov$F[1], df1=aov$gl[1], df2=aov$gl[2]), NA)
  names(aov) <- c(" ","Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)")
  return(list(anova=aov, R2=r2))
}

R2(n0)

#-----------------------------------------------------------------------------

Ajuste de muitos modelos de regressão linear

Pessoas que estão migrando de outro aplicativo para o R normalmente vivem perguntando “o R faz isso?, o R faz aquilo?…”. A resposta é “Sim, o R faz!”. Circula entre usuários o seguinte ditado: “não se pergunta o que o R faz, e sim como ele faz”.

Acontece que o R não é um software específico de ajuste de modelos de regressão, análise multivariada, análise de experimentos ou qualquer outra técnica. O R é geral. Os pacotes do R, por sua vez, são capazes de concentrar funções específicas. Um aplicativo que só faz ajuste de modelos de regressão tem espaço de sobra para incluir diversas opções, contem lista modelos para ajustar aos dados, enfim, mas tudo isso porque é destinado apenas para esse fim.

Por outro lado, no R é possível desenvolver procedimentos para qualquer tipo de tarefa. Se você quer ajustar uma lista de modelos e depois escolher o de maior R², ótimo! No você também consegue isso. Veja no CMR abaixo.

# gera dados
da <- data.frame(x=runif(100), z=5*rpois(100, lambda=7), w=runif(100, 50, 100))
da$y <- with(da, 12+0.1*x+0.05*z+0.34*w+0.2*sqrt(z)+0.1*x*w)+rnorm(100,0,0.1)
 
# vetor com as fórmulas específicando diferentes modelos lineares
form <- c(mod1=y~x, mod2=y~x+z, mod3=y~x+I(x^2), mod4=y~x+z+w)
 
# ajuste dos modelos
ajustes <- lapply(form, function(f){ m0 <- lm(f, data=da); m0 })
 
lapply(ajustes, summary) # quadro geral de estimativas e qualidade
lapply(ajustes, anova)   # quadro de anova sequencial
lapply(ajustes, coef)    # vetor de estimativas
sapply(ajustes, function(a){ summary(a)$r.squared})     # R²
sapply(ajustes, function(a){ summary(a)$adj.r.squared}) # R² ajustado
sapply(ajustes, function(a){ summary(a)$sigma})         # QMR
sapply(ajustes, deviance)                               # SQR
sapply(ajustes, df.residual)                            # GLR
lapply(ajustes, function(a){ summary(a)$coeff})         # tabela de estimativas
do.call(rbind, lapply(ajustes, function(a){ summary(a)$coeff})) # junta das tabelas
sapply(ajustes, fitted)    # valores ajustados
sapply(ajustes, residuals) # resíduos da análise
sapply(ajustes, vcov)      # matriz de covariância das estimativas
apply(sapply(ajustes, residuals), 2, shapiro.test) # normalidade dos resíduos

Sempre faça a avaliação das pressuposições do modelo antes de aplicar inferências. A escolha de modelo apenas pelo valor do R² não é aconselhada.