#===============================================================================
#. Título   : Análise quantitativa de dados - AQE
#. Curso    : CET076 - Metodologia e Estatística Experimental
# Autor    : José Cláudio Faria/UESC/DCET
# Data     : 30/08/2021 21:51:18
# 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 álgebra 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')