Análise de Regressão

## Objetivos:
## a) Apresentar os recursos básicos do R para análise de regressão
## b) Fundamentação em regressão linear simples
## c) Fundamentação em álgebra de matrizes
## Exemplos via algebra de matrizes

N   <- c(10, 20, 30, 40, 50, 60, 70)

Pro <- c(1000, 2300, 2600, 3900, 5400, 5800, 6600)

(X <- cbind(1, N))        # Ajustando um modelo linear de 1 grau 
#(X <- cbind(N))           # Ajustando um modelo linear de 1 grau forçando para a origem
#(X <- cbind(1, N, N^2))   # Ajustando um modelo linear de 2 grau

(Y <- cbind(Pro))

(XlX     <- t(X) %*% X)

(XlX_inv <- solve(XlX))

(XlY     <- t(X) %*% Y)

(b_est   <- XlX_inv %*% XlY) # vetor das estimativas dos parâmetros do modelo ajustado: inv(X'X) * (X'Y)

(Y_est   <- X %*% b_est)     # vetor dos valores Y estimados pelo modelo ajustado

(res     <- Y - Y_est)       # vetor dos erros

plot(Pro ~ N,
     xlim=c(0, 70),
     ylim=c(0, 7000),
     pch=19,
     col='darkgreen')

lines(Y_est ~ N,
      col='red')

points(Y_est ~ N,
       col='red',
       pch=19)


##===============================================================================
## ÊNFASE NAS FUNÇÕES DO R - USO ROTINEIRO
##===============================================================================
## Gerando dados de um experimento virtual
## x: variável independente (regressor/tratamento)
## y: variável dependente (fator resposta)

set.seed(10082006)
x <- sort(rep(1:5,
              6))             # 5 tratamentos e 6 repetições

r <- rep(1:6,
         5)                   # repetições

y <- 2 + 1 * x + rnorm(30, 
                       mean=0, 
                       sd=1)  # definindo a variável de resposta, segundo um modelo linear

## pré-concebido (y = 2 + 1*x)
## com uma pertubação (erro experimental) intencional adicional
(dad  <- data.frame(x, 
                    r, 
                    y))

## Visualizando os dados
with(dad,
     plot(y ~ factor(x),
          ylim=c(0, 10)))

## Observar que x visualmente influencia y

## ANOVA preliminar segundo delineamento inteiramente casualizado - DIC
av <- aov(y ~ factor(x),
          data=dad)

summary(av)              

## Pode concluir que o efeito dos tratamentos (níveis de x) é significativo.
## Como os tratamentos são quantitativos, a análise deve ser, necessariamente,
## conduzida via análise de regressão.

## Ou seja, não se pode realizar testes de comparação de médias múltiplas,
## nem análise de contrastes, para comparar os níveis de x.

## É necessário estudar se é possível estabelecer uma relação
## funcional de y em função de x. Em outras palavras, ajustar um modelo
## matemático que permita a predição de y a medida que x varia no intervalo
## considerado, além de testar a validade do modelo.

## Calculando as médias de tratamentos
dad_m <- aggregate(y ~ x,
                   data=dad,
                   mean)
dad_m

## Ajustando modelo linear e visualizando os resultados
with(dad_m,
     plot(y ~ x,
          ylab='m(y)',
          ylim=c(0, 10),
          pch=19,
          col='blue'))

## Polinomio de 1 grau
reg <- lm(y ~ x,
          data=dad_m)

## Polinomio de 2 grau
# reg <- lm(y ~ x + I(x^2),
#           data=dad_m)

summary(reg)

ANOVAreg <- anova(reg)
ANOVAreg

## Conclusão: o r2 é elevado e 98,54% da variação na variável de resposta y,
## decorrente da variação na variável independente x, é explicada
## pelo modelo linear (y = 1.62941 + 1.09960x) ajustado.

## Em outras palavras, o modelo linear ajustado explica satisfatoriamente as
## as variações da variável dependente y em decorrência da variação na variável
## independente (ou regressora, x).

## Predizer valores y a partir de valores x
newdad <- data.frame(x=c(1.5, 2.5, 3.5, 4.5))

y_pred <- predict(reg,
                  newdata=newdad)

y_pred

## Visualizar médias e modelo linear ajustado
lines(fitted(reg) ~ dad_m$x,
      col='red')

## Visualizar valores preditos
points(x=c(1.5, 2.5, 3.5, 4.5),
       y=y_pred,
       pch=19,
       col='red')