Página Inicial > gráficos > Equação junto ao gráfico na lattice

Equação junto ao gráfico na lattice

Reposta à salinidade para quatro espécies de craca. A regulação osmótica está ligada à capacidade adaptativa do organismo ao ambiente. O gráfico mostra os valores observados e a curva de regressão estimada acompanhada de bandas de confiança (95%). As equações em cada gráfico representam o modelo ajustado que descreve a curva.

Os polinômios são os modelos de regressão mais empregados em algumas áreas do conhecimento. Isso não quer dizer que eles sejam modelos bons. Na verdade, isso depende do cada um define como um bom modelo. Os polinômios são uma aproximação local. Isso significa que eles não devem ser usados para predições fora da amplitude observada das covariáveis (região experimental). Em outras palavras, não se deve fazer extrapolações, mas pode-se, com certos cuidados, fazer interpolações. Os polinômios são desprovidos de interpretação a partir do segundo grau. A partir do segundo grau, apenas o intercepto tem interpretação pois é o valor da função quando a covariável (x) é zero. Por isso, interpretações em geral são perdas de tempo. O que pode-se fazer são julgamentos quando ao sinal do coeficiente de maior grau que por exemplo, num modelo quadrático, indica se a curvatura é para cima (ou não). Outro ponto é que os polinômios apresentam problema de colinearidade. Para superar essa deficiência pode-se centrar a covariável ou usar polinômios ortogonais. Apesar do polinômio já ser sem interpretação, o uso de polinômios ortogonais é menos interpretável ainda pois os coeficientes não são taxas em relação à x, x², x³…, mas para uma transformação linear de x, essa que induz a ortogonalidade. Resumindo, os polinômios ortogonais são úteis pois minimizam o problema de colineridade. Então, do ponto de vista de reportar uma equação de predição, aquela que descreve o valor esperado de y (resposta) em função de x, é mais útil usar o polinômio cru. Isso porque, mesmo sem ter interpretação, toda pessoa sabe traçar a curva de um polinômio, é algo familiar, visto desde o ensino médio, e novamente, não que isso seja um ponto positivo para um modelo. De qualquer forma, o post de hoje ensina a ajustar polinômios e representar no gráfico os valores observados, os valores ajustados com intervalo de confiança e a equação dentro do gráfico. Os procedimentos abaixo eu desenvolvi para atender demandas de consultoria. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# definições da sessão

require(lattice)
require(latticeExtra)
require(grid)

#-----------------------------------------------------------------------------
# leitura

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

xyplot(Osm~Salinidade|Espécie, data=osmo, type=c("p","smooth"))

#-----------------------------------------------------------------------------
# vou desconsiderar o fato dos dados serem heterocedasticos e ajustar um lm
# em outra matéria trataremos os desvios de pressupostos
# ajuste usando polinômios ortogonais

m0 <- aov(Osm~Espécie*poly(Salinidade,2), data=osmo)
par(mfrow=c(2,2)); plot(m0); layout(1)

qqmath(~residuals(m0)|osmo$Espécie)
xyplot(sqrt(abs(residuals(m0)))~fitted(m0)|osmo$Espécie, type=c("p","smooth"))
# desvio de pressupostos, trataremos em outra matéria do blog

summary(m0) # inferências sem validade pois existe desvio de pressupostos

#-----------------------------------------------------------------------------
# partição das somas de quadrados, declara Salinidade dentro de Espécie

m1 <- aov(Osm~Espécie/poly(Salinidade,2), data=osmo)

names(coef(m1))
attr(terms(m1), "term.labels")
ncoef <- 2 # número de efeitos por expécie=2 (linear e quadrático)
ngrup <- 4 # número de níveis de espécie=4
split.list <- split(matrix(1:(ncoef*ngrup), ngrup, ncoef), f=1:ngrup)
names(split.list) <- levels(osmo$Espécie)
split.list <- lapply(split.list,
                     function(x){ names(x) <- c("lin","qua"); x})
split.list <- list(a=split.list)
names(split.list) <- attr(terms(m1), "term.labels")[2]

summary(m1, split=split.list)

split.list <- as.list(1:(ncoef*ngrup))
names(split.list) <- c(outer(levels(osmo$Espécie), c("lin","qua"),
                             paste, sep=": "))
split.list <- list(a=split.list)
names(split.list) <- attr(terms(m1), "term.labels")[2]

summary(m1, split=split.list)

#-----------------------------------------------------------------------------
# prepara as equações para colocar no gráfico, função para polinômiais

coefs <- coef(lm(Osm~-1+Espécie/(Salinidade+I(Salinidade^2)), data=osmo))
coefs <- as.numeric(formatC(coefs, format="f", digits=3, flag="-"))

get.expression <- function(coef, digits){
  coef <- as.numeric(formatC(coef, format="f", digits=digits, flag="-"))
  postext <- c("","*x", paste("*x^", 2:(length(coef)-1), sep=""))
  sig <- c(ifelse(coef[1]<0, "-", ""), ifelse(coef[-1]<0, "-", "+"))
  txt <- paste(paste(sig, abs(coef), postext, sep=""), collapse=" ")
  eval(parse(text=paste("expression(y==", txt, ")", sep="")))
}

coefs <- split(matrix(coefs, ncol=3), f=1:4)
eqn <- lapply(coefs, get.expression, digits=3)
eqn

#-----------------------------------------------------------------------------
# nomes científicos em itálico nas tarjas

fl <- as.expression(lapply(levels(osmo$Espécie),
    function(x){ bquote(italic(.(x))) }))
fl

#-----------------------------------------------------------------------------
# predição

pred <- expand.grid(Espécie=levels(osmo$Espécie), Salinidade=seq(0,40,by=2))
pred <- cbind(pred, predict(m1, newdata=pred, interval="confidence"))

#-----------------------------------------------------------------------------
# carrega funções que serão necessárias para fazer as bandas de confiança

source("http://www.leg.ufpr.br/~walmes/ridiculas/ridiculas_functions.R")
grep("ciH", ls(), value=TRUE) # nome das funções (painéis) que vamos usar

#-----------------------------------------------------------------------------
# gráfico final

#png("f035.png", width=600, height=500)
xyplot(Osm~Salinidade|Espécie, data=osmo, col=1,
       xlab=expression(Salinidade~(ppm)),
       ylab=expression(Osmolalidade~(mOsm~kg^{-1}~de~H[2]*O)),
       strip=strip.custom(bg="gray90", factor.levels=fl))+
  as.layer(xyplot(fit~Salinidade|Espécie, data=pred, type="l", col=1,
                  ly=pred$lwr, uy=pred$upr,
                  prepanel=prepanel.ciH, panel=panel.ciH))
dimLay <- trellis.currentLayout()
for(col in 1:ncol(dimLay)){
  for(row in 1:nrow(dimLay)){
    trellis.focus("panel", col, row, highlight=FALSE)
    grid.text(label=eqn[[dimLay[row,col]]],
              x=0.03, y=1-0.03, just=c("left","top"))
  }
}
trellis.unfocus()
#dev.off()

#-----------------------------------------------------------------------------
About these ads
Categoriasgráficos Tags:, ,
  1. Nenhum comentário ainda.
  1. No trackbacks yet.

Deixe um comentário

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s

Seguir

Obtenha todo post novo entregue na sua caixa de entrada.

Junte-se a 52 outros seguidores

%d blogueiros gostam disto: