## 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')
