Arquivo

Posts Tagged ‘mtext’

Como fazer legendas em gráficos

Gráficos com legendas e assinatura.

 

Para os gráficos do pacote graphics a adição de legendas é feita por meio da função legend(). Essa função possuí diversas opções para confecção de diversos tipos de legenda (veja os exemplos da documentação).

Nessa ridícula eu fiz alguns exemplos de como colocar legendas em gráficos de funções, gráfico com pontos e retas, gráfico do ajuste de um modelo de regressão, gráfico de barras, legenda fora da área gráfica, escolha da posição da legenda com o mouse e como colocar assinatura nos gráficos. A documentação da função esclarece o significado dos argumentos usados. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# legenda para uma função

curve(x+x^2, 0, 3, lty=1, col=2, lwd=2)
legend("topleft", legend="f(x)", lty=1, col=2, lwd=2, bty="n")

#-----------------------------------------------------------------------------
# legenda para gráfico de pontos ligados

x <- 1:10
y <- x+rnorm(x,0.5)
plot(y~x, type="o", pch=19)
legend("bottomright", legend="valores observados",
       lty=1, col=1, bty="n", pch=19)

#-----------------------------------------------------------------------------
# legenda para pontos e reta de regressão

m0 <- lm(y~x)
plot(y~x)
abline(m0, col=2, lty=2, lwd=2)
legend("top", legend=c("valores observados", "valores ajustados"),
       lty=c(NA,2), col=c(1,2), lwd=1:2, bty="n", pch=c(1,NA))

#-----------------------------------------------------------------------------
# legenda para gráficos de barras

z <- matrix(1:4, 2, 2)
colnames(z) <- c("A","B")
barplot(z, beside=TRUE, col=c("forestgreen", "palegreen"))
legend("topleft", legend=c("tipo 1", "tipo 2"),
       fill=c("forestgreen", "palegreen"), bty="n")

#-----------------------------------------------------------------------------
# legenda fora da área gráfica disposto em colunas

barplot(z, beside=TRUE, col=c("forestgreen", "palegreen"))
legend(x=1, y=4.5, xpd=TRUE, ncol=2, legend=c("tipo 1", "tipo 2"),
       fill=c("forestgreen", "palegreen"), bty="n")

#-----------------------------------------------------------------------------
# escolhendo a posição da legenda com o mouse

barplot(z, beside=TRUE, col=c("forestgreen", "palegreen"))
legend(locator(1), xpd=TRUE, ncol=2, legend=c("tipo 1", "tipo 2"),
       fill=c("forestgreen", "palegreen"), bty="n")

#-----------------------------------------------------------------------------
# nota de rodapé no gráfico com assinatura, data e versão

barplot(z, beside=TRUE, col=c("forestgreen", "palegreen"))
legend("topleft", xpd=TRUE, ncol=2, legend=c("tipo 1", "tipo 2"),
       fill=c("forestgreen", "palegreen"), bty="n")
mtext(side=4, at=par("usr")[3], adj=0,
      cex=0.8, col="gray40", line=-1,
      text=paste("Walmes Zeviani --",
        format(Sys.time(), "%d/%m/%Y %H:%M:%S --"),
        R.version.string))

#-----------------------------------------------------------------------------
Anúncios
Categorias:gráficos Tags:, , , ,

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)

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