Arquivo

Posts Tagged ‘predição’

Ajuste simultâneo de regressão/anova para várias respostas

Valores observados e ajustados acompanhados do intervalo de 95% em função da dose de potássio e nível de água para as variáveis resposta estudadas.

Na realização de experimento é comum a observação de várias variáveis respostas em uma mesma unidade amostral. Em termos estatísticos, o que observamos é um vetor aleatório. Por exemplo, em experimento de com soja para fins de ensaio de competição de cultivares o pesquisador pode observar: altura da planta, altura da primeira vagem, massa seca da planta, número de vagens por planta, teores foliares de N, P, K, Mg, Ca, peso de 100 grãos, proporção de sementes viáveis, dias para fechar o ciclo, todas elas observadas em cada unidade experimental. Nesse caso então, observamos um vetor aleatório de 12 variáveis aleatórias. O ideal seria empregarmos um modelo multivariado, porém como as respostas podem ser contínuas não limitadas, contagens, contínuas limitadas, contínuas com censura, proporções, etc, ainda não dispomos de um modelo de distribuição multivariado que contemple essa situação. O mais comum é fazer a análise univariada e como todas essas variáveis estão sob o efeito das mesmas fontes de variação, a análise estatística vai empregar o mesmo modelo (pelo menos no tocante a parte determinística). Sendo assim, tem como fazer o estudo de todas as variáveis reposta de uma vez só? Sim. Para saber como continue lendo.

Existem diversas formas de aplicar tarefas em grupos ou colunas. A família apply do stats oferece recurso para um número bem grande de situações. A família ply do plyr acrescenta mais opções. Além destes temos o doBy com mais algumas funcionalidades.

Bem, vamos ao que interessa. No CMR eu faço ajustes de regressão com polinômios para dados de um experimento com a cultura da soja (vasos em casa de vegetação). Foram estudados o nível de água (umidade do solo em relação à capacidade de campo) e potássio aplicado ao solo. O objetivo desse experimento foi verificar se o fornecimento de potássio pode minimizar impactos do deficit hídrico. Para acessar o artigo completo clique aqui: Umidade do solo e doses de potássio na cultura da soja. Não deixe de conferir a lista de links abaixo. Em próximos artigos do Rídiculas vamos aprofundar o estudo nesses dados comparando curvas. Até lá.

  • A brief introduction to “apply” in R;
  • plyr – The split-apply-combine strategy for R;
  • plyr: Tools for splitting, applying and combining data;
  • R Grouping functions: sapply vs. lapply vs. apply. vs. tapply vs. by vs. aggregate vs.;
  • Say it in R with “by”, “apply” and friends;
  • The doBy package;
  • doBy: doBy – Groupwise summary statistics, general linear contrasts, population means (least-squares-means), and other utilities;
  • #-----------------------------------------------------------------------------
    # importa dados
    soja <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/soja.txt",
                       header=TRUE, sep="\t", dec=",")
    str(soja)
    soja <- transform(soja, a=agua, A=factor(agua),
                      k=potassio, K=factor(potassio))
    resp <- 4:7 # índice das variáveis resposta
    str(soja)
    
    #-----------------------------------------------------------------------------
    # ver
    
    soja2 <- cbind(stack(soja[,resp]), soja[,-resp])
    
    require(lattice)
    
    xyplot(values~k|ind, groups=A, data=soja2,
           scales="free", jitter.x=TRUE, type=c("p","a"))
    
    #-----------------------------------------------------------------------------
    # ajuste do modelo saturado de efeitos categóricos (modelo de médias)
    
    m0 <- apply(soja[,resp], 2,
                function(r){ aov(r~bloco+A*K, soja) })
    lapply(m0, anova)
    
    #-----------------------------------------------------------------------------
    # ajuste de polinômio para potassio por nível de água
    
    m1 <- apply(soja[,resp], 2,
                function(r){ aov(r~bloco+A/poly(k, 4), soja) })
    lapply(m1, anova)
    names(coef(m1[[1]]))
    
    inter <- m1[[1]]$assign==3        # indices das interações
    eff <- names(coef(m1[[1]]))[inter]
    togrep <- outer(c("A37.5","A50","A62.5"), c("1$","2$","(3|4)$"),
                    paste, sep=":.*") # nomes de busca
    or <- order(togrep)
    tolist <- outer(c("A37.5","A50.0","A62.5"), c("lin","qua","lof"),
                    paste, sep="|")   # nomes de rótulo
    list.desd <- sapply(togrep, grep, eff)
    names(list.desd) <- tolist
    list.desd <- list.desd[or]
    list.desd                         # lista com os índices dos desdobramentos
    
    # anovas com desdobramento de efeito de potássio em cada nível de água
    lapply(m1, summary, split=list("A:poly(k, 4)"=list.desd))
    
    #-----------------------------------------------------------------------------
    # com base nas anovas acima selecionamos os modelos
    
    form <- list( rengrao~bloco+A/poly(k, 2, raw=TRUE),
                 pesograo~bloco+A/poly(k, 2, raw=TRUE),
                    kgrao~bloco+A/poly(k, 2, raw=TRUE),
                    pgrao~bloco+A/poly(k, 3, raw=TRUE))
    names(form) <- names(soja)[resp]
    
    m2 <- lapply(form, aov, data=soja, contrasts=list(bloco=contr.sum))
    lapply(m2, coef)
    
    #-----------------------------------------------------------------------------
    # predição a partir do modelo
    
    str(soja)
    pred <- expand.grid(bloco="I", A=levels(soja$A), k=seq(0,180,5))
    
    #aux <- lapply(m2, predict, newdata=pred, interval="confidence")
    aux <- lapply(m2, function(r){ # remove o efeito do bloco na predição
      predict(r, newdata=pred, interval="confidence")-coef(r)["bloco1"]})
    aux <- lapply(aux, as.data.frame)
    str(aux)
    
    require(plyr)
    pred <- cbind(ldply(aux), pred)
    pred$.id <- factor(pred$.id, levels=levels(soja2$ind))
    
    #-----------------------------------------------------------------------------
    # gráfico final
    
    fl <- c("K nos grãos (mg)", "Peso 100 grãos (g)",
            "P nos grãos (mg)", "Rendimento de grãos (g)")
    cbind(levels(pred$.id), levels(soja2$ind), fl)
    
    require(latticeExtra)
    display.brewer.all()
    mycol <- brewer.pal(6, "Set1")
    
    panel.ciH <- function(x, y, ly, uy, subscripts, ...){
      y <- as.numeric(y)
      x <- as.numeric(x)
      or <- order(x)
      ly <- as.numeric(ly[subscripts])
      uy <- as.numeric(uy[subscripts])
      panel.polygon(c(x[or], x[rev(or)]), c(ly[or], uy[rev(or)]),
                    col=1, alpha=0.05, border=NA)
      panel.lines(x[or], ly[or], lty=3, lwd=0.5, col=1)
      panel.lines(x[or], uy[or], lty=3, lwd=0.5, col=1)
      panel.xyplot(x, y, subscripts=subscripts, ...)
    }
    
    #png(file="f031.png", width=500, height=500)
    trellis.par.set(superpose.symbol=list(col=mycol),
                    superpose.line=list(col=mycol),
                    strip.background=list(col="gray90"))
    xyplot(values~k|ind, groups=A, data=soja2,
           scales="free", jitter.x=TRUE,
           xlab="Potássio no solo", ylab=NULL,
           strip=strip.custom(factor.levels=fl),
           auto.key=list(title="Níveis de água (%)", cex.title=1.2,
             columns=3, points=FALSE, lines=TRUE))+
      as.layer(xyplot(fit~k|.id, groups=A, pred, type="l", lwd=1.5,
                      ly=pred$lwr, uy=pred$upr,
                      panel.groups=panel.ciH,
                      panel=panel.superpose))
    #dev.off()
    
    #-----------------------------------------------------------------------------
    # pode-se tentar comparar as curvas ajustando modelos apropriados sob H0
    # exemplo, será que as curvas vermelha e verde para kgrao e pesograo não
    # podem ser representados pela mesma curva?
    # as curvas azul e ver para pgrao diferem de uma reta horizontal?
    
    #-----------------------------------------------------------------------------
    
    Categorias:regressão Tags:, , ,

    Como fazer regressão polinomial com diferentes graus

    Ajuste de modelos de regressão polinomial com diferentes graus para o polinômio. Para as cultivares da esquerda o polinômio é quadrático, para a da direita o polinômio é cúbico.

    No post Predição em modelos de regressão com blocos foram apresentados procedimentos para estimação e predição de parâmetros em modelos de regressão com blocos. Nesse post aqui eu vou estender o exemplo permitindo o ajuste com diferentes graus para o polinômio que descreve a relação da resposta com a preditora contínua. Nos códigos abaixo vou apresentar aspectos dos polinômios ortogonais, do cálculo do R² e de como fazer o ajuste com diferentes graus de polinômio em um mesmo modelo, ou em outras palavras, em uma mesma lm(). No final, como de costume, confecciono o gráfico resultado do modelo ajustado com as bandas de confiança. Mais uma vez deixei o código comentado passo a passo para orientar o leitor. Até a próxima ridícula!

    #-----------------------------------------------------------------------------
    # dados de índice agronômico de milho de um experimento
    
    da <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/anovareg.txt",
                     header=TRUE, sep="\t")
    str(da)
    
    #-----------------------------------------------------------------------------
    # objetivo: apresentar procedimentos para ajuste de polinômios e como ajustar
    # polinômios de diferentes graus em uma mesma análise
    
    #-----------------------------------------------------------------------------
    # vamos apresentar essas funções para o caso do ajuste de uma função e fazer
    # de conta que não temos blocos apenas por facilidade de exposição
    
    levels(da$cultivar)
    br <- subset(da, cultivar=="BR-300")
    str(br)
    
    #-----------------------------------------------------------------------------
    # ajustar diminuindo o nível do polinômio até que não tenhamos falta de ajuste
    # fazer isso de diferentes formas:
    # 1) fórmula explicita e anova seqüencial
    # 2) polinômios ortogonais e teste t
    
    # 1) fórmula explicita e anova seqüencial
    m1.5 <- lm(indice~dose+I(dose^2)+I(dose^3)+I(dose^4)+I(dose^5), data=br)
    anova(m1.5) # testes F apontam que os termos > 2 grau não contribuem
    m1.2 <- lm(indice~dose+I(dose^2), data=br)
    anova(m1.2, m1.5) # o abandono desses termos não resulta em falta de ajuste
    anova(m1.2)
    
    # 2) polinômios ortogonais e teste t
    m2.5 <- lm(indice~poly(dose, degree=5), data=br)
    anova(m2.5)   # anova não separa os termos
    summary(m2.5) # quando usamos polinômios orotoganis, o valor t^2 = F
    cbind("t^2"=summary(m2.5)$coeff[-1,3]^2, "F"=anova(m1.5)[-6,4]) # viu?
    
    m2.2 <- lm(indice~poly(dose, degree=2), data=br)
    summary(m2.2)
    
    #-----------------------------------------------------------------------------
    # antigamente experimentos agronômicos com fatores quantitativos eram
    # planejados para polinômios ortogonais. assim usava espaçamento regular entre
    # níveis. muitas pessoas acham que isso é regra mas era apenas algo imposto
    # pelas limitações de computação da época. com polinômios ortogonais a matriz
    # X'X é ortogonal, então facilmente invertível e obtém-se as estimativas dos
    # parâmetros. com essas estimativas e erros padrões aplicava-se o teste t (ao
    # invés do teste F) para testar a significância dos termos. os temos não
    # significativos eram abandonados mas a estimação não precisava ser refeita
    # por causa da ortogonalidade. observe os procedimentos abaixo.
    
    X <- model.matrix(m2.5) # matriz de delineamento ortogonal
    round(colSums(X), 4) # soma nas colunas dá zero
    round(t(X)%*%X, 4)   # ortogonal, portanto a inversa é a inversa da diagonal
    t(X)%*%br$indice     # as estimativas dos parâmetros, intercepto dividir por n
    
    cbind(X[,2:3], model.matrix(m2.2)[,2:3])[1:6,] # iguais colunas
    coef(m2.5); coef(m2.2) # as estimativas são iguais, não requer outro ajuste
    
    #-----------------------------------------------------------------------------
    # no entanto, esses coeficientes não tem interpretação pois a matriz não está
    # na escala das doses. aliás, polinômio não tem interpretação em nenhuma
    # escala. eles representam apenas uma aproximação local. interpretar a
    # magnitude das estimativas é perda de tempo. compare apenas o rank e a
    # direção (sinal). faça o teste de hipótese e interprete os valores preditos.
    # caso queira apresentar a equação ajustada, terá que obter as estimativas
    # sem polinômios ortogonais.
    
    #-----------------------------------------------------------------------------
    # obter as estimativas sem polinômios ortogonais com a função poly()
    
    m2.2 <- lm(indice~poly(dose, degree=2, raw=TRUE), data=br)
    summary(m2.2) # estimativas na escala do fator, são as mesmas do modelo m1.2
    
    #-----------------------------------------------------------------------------
    # muitos tutoriais/aplicativos obtém o R² do ajuste fazendo ajuste às médias.
    # *eu* não considero essa prática porque esse R² não reflete a razão
    # ruído/sinal dos dados para com o modelo, mas da média dos dados para com o
    # modelo. esse tipo de procedimento só é possível quando temos repetição de
    # níveis. veja os três exemplos simulados a seguir
    
    desv <- rnorm(12,0,1)      # desvios aleatórios
    a <- 0; b <- 1             # valores dos parâmetros
    x1 <- seq(1,12,1)          # sem repetição de nível
    x2 <- rep(seq(3,12,3),e=3) # com repetição de nível
    y1 <- a+x1*b+desv          # valores observados de y1
    y2 <- a+x2*b+desv          # valores observados de y2
    
    m1 <- lm(y1~x1); summary(m1) # R² correto
    m2 <- lm(y2~x2); summary(m2) # R² correto
    
    # R² calculado com relação ao desvio dos ajustados para as médias
    md <- tapply(y2, x2, mean)            # tiramos a média das repetições
    aj <- unique(round(fitted(m2),5))     # valores ajustados pelo modelo
    r2 <- 1-sum((md-aj)^2)/sum((md-mean(md))^2); r2 # R² como desvio das médias
    
    #-----------------------------------------------------------------------------
    # não sei dizer qual a justificativa para o uso desse procedimento para obter
    # o R². é um procedimento atraente porque o R² será *sempre* maior e talvez
    # essa seja a razão da sua adoção (que na *minha opinião* é errada). esse R²
    # não associa desvio das observações. então o valor desse R² não é a
    # porcentagem de explicação que o modelo dá para as observações, e sim para a
    # média delas! como assim se estou modelando as observações e não a média
    # amostral delas? veja esse exemplo mais abrupto
    
    #-----------------------------------------------------------------------------
    # criando um novo par x, y
    
    x3 <- x2+c(-0.001,0,0.001) # deslocamento mínimo faz perder as repetições
    y3 <- a+x3*b+desv
    
    cbind(y2,y3,x2,x3)
    
    # ambos pares x, y são praticamente a mesma coisa e o R² será muito próximo
    # para ambos. porém, para x2 pode-se usar o R² relativo a desvio da média.
    # se eu calcular esse R² ele vai diferir muito pois representa outra razão
    # ruído/sinal. *eu* não gosto disso.
    
    m3 <- lm(y3~x3)
    c(summary(m1)$r.sq, r2, summary(m2)$r.sq, summary(m3)$r.sq) # compara R²
    
    #-----------------------------------------------------------------------------
    # para mais críticas sobre R² procure sobre a análise do quarteto de Anscombe.
    
    #-----------------------------------------------------------------------------
    # ajustar modelos de regressão para os 3 níveis de cultivar e decompor a SQ.
    # aqui usaremos o fator de níveis ordenados que tem como padrão o contraste
    # de polinômios ortogonais. na poly() conseguimos especificar o grau e no
    # fator ordenado não conseguimos, mas as somas de quadrados podem ser
    # facilmente desdobradas se usarmos aov().
    
    m0 <- aov(indice~bloco+cultivar*ordered(dose), data=da)
    summary(m0)    # pode-se usar anova(m0) também
    summary.lm(m0) # quadro de estimativas dos parâmetros
    
    #-----------------------------------------------------------------------------
    # para obter o quadro de análise de variância com desdobramentos dos termos
    # do polinômio dentro dos níveis de cultivar precisamos mudar a fórmula do
    # modelo, é apenas uma mudança a nível de colunas mas o espaço do modelo ainda
    # é o mesmo. compare as colunas das matrizes obtidas com os comando abaixo.
    
    aux <- expand.grid(cultivar=gl(2,1), dose=1:3) # estrutura experimental
    model.matrix(~cultivar*ordered(dose), aux)     # cultivar cruzado com dose
    model.matrix(~cultivar/ordered(dose), aux)     # dose aninhado em cultivar
    
    #-----------------------------------------------------------------------------
    # ajustando modelo com fórmula corresponde ao desdobramento desejado
    # o procedimento abaixo é para criar a lista que especifica a partição nas
    # somas de quadrados, envolve algumas funções que vale a pena consultar o help
    
    da$D <- ordered(da$dose) # cria uma nova coluna apenas por conforto
    m1 <- aov(indice~bloco+cultivar/D, data=da) # modelo com a fórmula conveniente
    
    m1$assign # valores iguais a 3 representam os estimativas do 3 termo
    tm <- names(coef(m1))[m1$assign==3]                  # nomes dos termos
    tm <- gsub("^cultivar(.*):.*(.{2})$", "\\1-\\2", tm) # edita os nomes
    ord <- c(sapply(levels(da$cultivar), grep, x=tm, simplify=TRUE)) # ordena
    split <- as.list(ord)   # cria a lista com a posição dos termos
    names(split) <- tm[ord] # atribui nome aos elementos da lista
    
    #-----------------------------------------------------------------------------
    # quadro de análise de variância com somas que quadrados particionadas.
    # com ela escolhemos o grau do polinômio adequado para cada cultivar.
    
    summary(m1, split=list("cultivar:D"=split))
    
    #-----------------------------------------------------------------------------
    # uma das cultivares requer 3 grau, as outras ficam com o 2 grau. juntar os
    # termos não significativos do polinômio para criar o termo falta de ajuste.
    
    do.call(c, split) # ver os índices e fazer o agrupamento gerando nova split
    split <- list(Ag.L=1, Ag.Q=4,         Ag.lof=c(7,10,13),
                  BR.L=2, BR.Q=5,         BR.lof=c(8,11,14),
                  Pi.L=3, Pi.Q=6, Pi.C=9, Pi.lof=c(12,15))
    summary(m1, split=list("cultivar:D"=split)) # lof não significativos, ok!
    
    #-----------------------------------------------------------------------------
    # não iremos apresentar as estimativas do modelo m1 nem seus valores ajustados
    # porque esses se referem ao modelo completo. temos que ajustar o modelo
    # reduzido removendo os temos apontado pelo teste F. este modelo é o modelo
    # que apresentaremos as estimativas e faremos predição. note que o grau do
    # polinômio agora muda com o nível de cultivar. teremos que arrumar um
    # artifício ajustar um modelo considere isso. então vou criar uma indicadora
    # para multiplicar o termo cúbico e sumir com ele quando não for necessário.
    
    levels(da$cultivar)
    da$ind <- 0; da[da$cultivar=="Pioneer-B815",]$ind <- 1 # a indicadora do nível
    da$ind # o valor 1 é associado ao nível de cultivar que tem o termo cúbico
    
    #-----------------------------------------------------------------------------
    # o termo cúbico só vai existir para o nível Pionner por causa da indicadora.
    # vai aparecer um NA para os outros níveis
    
    m2 <- lm(indice~bloco+cultivar*(dose+I(dose^2)+I(ind*dose^3)), data=da)
    summary(m2)
    
    #-----------------------------------------------------------------------------
    # os NA implicam em algumas conseqüência como no funcionamento da predict().
    # para evitá-los podemos criar uma coluna para representar esse termo cúbico e
    # e declara uma fórmula de modelo que não retorne NA.
    
    da$Pi.C <- with(da, ind*dose^3) # é a coluna do cubo da dose para Pioneer
    m3 <- lm(indice~bloco+cultivar*(dose+I(dose^2))+Pi.C, data=da,
             contrast=list(bloco=contr.sum))
    summary(m3)
    
    #-----------------------------------------------------------------------------
    # faz o gráfico final de predição pelo método matricial
    
    # data.frame da predição
    pred <- expand.grid(cultivar=levels(da$cultivar), dose=seq(0,300,10))
    pred$Pi.C <- with(pred, ifelse(cultivar=="Pioneer-B815", dose^3, 0))
    
    # matriz de predição sem colunas para os efeitos de blocos
    X <- model.matrix(~cultivar*(dose+I(dose^2))+Pi.C, data=pred)
    head(X)
    
    # vetor de estimativas sem os efeitos de blocos
    b <- coef(m3)[-grep("bloco", names(coef(m3)))]
    
    # vetor de erros padrões das estimativas
    se <- predict(m3, newdata=cbind(bloco=levels(da$bloco)[1], pred),
                  se.fit=TRUE)$se
    
    # quantil superior da distribuição t para um IC de 95%
    tc <- qt(0.975, df=df.residual(m3))
    
    # margem de erro
    me <- se*tc
    
    # data.frame com valores preditos e IC
    pred$fit <- X%*%b
    pred$lwr <- pred$fit-me; pred$upr <- pred$fit+me
    
    require(reshape)
    
    pred <- melt(pred[,-3], id=c("cultivar","dose"))
    names(pred)[ncol(pred)] <- "indice"
    da$variable <- "obs"
    pred <- rbind(da[,names(pred)], pred)
    
    #-----------------------------------------------------------------------------
    # gráfico com observados, preditos e IC pelo método matricial
    
    require(lattice)
    
    #png("f019.png", w=500, h=250)
    xyplot(indice~dose|cultivar, groups=variable, data=pred,
           distribute.type=TRUE, layout=c(3,1),
           type=c("l","l","p","l"), col=c(2,3,1,3), lty=c(1,2,0,2),
           strip=strip.custom(bg="gray90"),
           xlab="Dose de nitrogênio (kg/ha)", ylab="Índice agronômico",
           sub=list("predição livre do efeito de bloco", cex=0.8, font=3),
           key=list(space="top", columns=3, type="o", divide=1,
             lines=list(pch=c(1,NA,NA), lty=c(0,1,2), col=c(1,2,3)),
             text=list(c("valores observados","valores preditos","IC 95%"))))
    #dev.off()
    
    #-----------------------------------------------------------------------------
    

    Expressões, bandas de confiança e legenda

    Nesse post vou apresentar como confeccionar um gráfico para representar o ajuste de um modelo de regressão. Esse gráfico vai conter elementos básicos de um ajuste, ou seja, valores observados, valores preditos, bandas de confiança, legenda e equação estimada (veja figura abaixo), que são elementos que normalmente aparecem em gráficos de artigos científicos e gráficos gerados por outros aplicativos estatísticos (menos interessantes que o R, lógico!).

    Valores observados e preditos pelo modelo de regressão com o intervalo de confiança (95%).

    Vamos usar os dados contidos no data.frame “anscombe”. Esse conjunto de dados é usado em livros de regressão para exemplificar a utilidade das medidas de diagnóstico de resíduos. Em outra oportunidade vou postar algo sobre isso.

    No script abaixo, na linha [10] foi ajustado o modelo de equação da reta

    y_{1i} = \beta_0 + \beta_1 \cdot x_{1i} + \epsilon_i.

    aos dados. O resumo desse ajuste é apresentado em [11]. Foi feito a predição da resposta segundo a equação acima e as estimativas dos parâmetros no passo [16:17]. Usamos um intervalo de 50 pontos para obtermos linhas suaves. As estimativas dos parâmetros e R² foram arredondadas e armazenadas em um vetor [19].

    No último bloco do script, trocamos a opção para que os decimais sejam representados por vírgula tanto no console quanto nos gráficos [24], pois na língua portuguesa o separador decimal é a virgula. Em [25] fazemos o diagrama de dispersão dos dados observados, adicionando as retas do ajuste e das bandas de confiança em [26]. Em [27:30] é colocada a legenda ao gráfico no canto inferior direito, identificando com ponto de linhas os valores observados, valores preditos e as bandas de confiança de 95% para o valor predito. Finalmente, em [31:33] é adicionado à equação ajustada ao gráfico. Usamos a função locator() para que o usuário clique no gráfico para adicionar a equação. A função substitute() foi usada para automaticamente substituir os valores dos parâmetros por suas estimativas na equação.

    #-----------------------------------------------------------------------------
    # dados
    
    data(anscombe)
    str(anscombe)
    
    #-----------------------------------------------------------------------------
    # ajuste do modelo
    
    m0 <- lm(y1~x1, data=anscombe)
    summary(m0)
    
    #-----------------------------------------------------------------------------
    # fazendo a predição intervalar num grid regular mais fino
    
    pred <- data.frame(x1=seq(min(anscombe$x1), max(anscombe$x1), length=50))
    pred$y1 <- predict(m0, newdata=pred, interval="confidence")
    str(pred)
    cf <- format(c(coef(m0), summary(m0)$r.squared), digits=4)
    
    #-----------------------------------------------------------------------------
    # fazendo o gráfico, descomente para salvar a figura em png
    
    #png("f001.png", w=500, h=300); par(mar=c(5.1,4.1,2.1,2.1))
    options(OutDec=",")
    plot(y1~x1, data=anscombe, xlab=expression(x[1]), ylab=expression(y[1]))
    with(pred, matlines(x1, y1, col=c(1,2,2), lty=c(1,2,2)))
    legend("bottomright",
           legend=c("valores observados", "valores preditos",
             "intervalo de confiança (95%)"),
           lty=c(NA,1,2), col=c(1,1,2), pch=c(1,NA,NA), bty="n")
    text(locator(n=1), # ao salvar a figura use x=7.5, y=10 no lugar da locator()
         label=substitute(hat(y)[1]==a+b%.%x[1]~~~(R^2==c),
           list(a=cf[1], b=cf[2], c=cf[3]))) # clicar no gráfico
    #dev.off()
    
    #-----------------------------------------------------------------------------