Vamos revisar, comentar e alterar um código enviado na lista brasileira do do R para ilustrar alguns aspectos do uso adequado da linguagem.

Aproveitamos o exemplo para discutir a especificação de modelos através do uso de fórmulas no R.

## definindo número de dados  
idados<- 50
##adados<-20000
 
# Equação verdadeira de relacionamento 
## deseja simular de y = beta_0 + beta_1 X1 + beta_2 X2 + beta_3 X1.X2 
# (usada para simular Fase I e gerar Modelo de Regressão)
coef0 <- c(10, 2, -3, 2.3)    ## parâmetros para simular o modelo
 
#Equação alterada
#coefal<-c(10, 6, -3, 2.3)
 
#Fase I
## simulando o valor das covariáveis X1 e X2
x1 <-runif(idados, min=-1, max=1)
x2 <-runif(idados, min=-1, max=1)
## simulando os erros:  
e <-rnorm(idados, mean=0, sd=4)
 
## criando a matriz X 
x <-cbind(1, x1, x2, x1*x2)
## Comentários: 
##     não precisa criar com a coluna x1*x2 
##     pode usar a função model.matrix() para montar a matriz a partir
##         da espeficicação do modelo
X <- model.matrix(~ x1*x2)
X
 
## Sobre formulas no R:
## Notar que os operadores +, * , ^ etc, tem um sentido diferente em formulas
## I()  "força" a operação algébrica entre as variáveis
head(model.matrix(~ x1*x2))
head(model.matrix(~ I(x1*x2)))
head(model.matrix(~ x1+x2))
head(model.matrix(~ I(x1+x2)))
 
 
#invx<-solve(crossprod(x)) #operação (X'X)^-1
da <-data.frame (x1, x2, x1*x2)
 
## colocando as covariáveis num data-frame
da <-data.frame (x1=x1, x1=x2)
## Com.:note que poderi já ter criado aqui para evitar objetos desnecessários:
da <-data.frame(x1=runif(idados, min=-1, max=1), 
                 x2=runif(idados, min=-1, max=1)
               )
 
## simulando os dados
y  <- vector(length=idados)
ii<-1:idados
for (ii in 1:idados) {
  y[ii]<- coef0%*%x[ii,]+e[ii]
}
 
## Com: note que não precisa fazer o loop for() pois o R é linguagem vetorial
##       alem disto poderia criar o y direto no data-frame
da$y <- X %*% coef0 + e
da$y  <-y
 
all.equal(y, da$y)
 
## Com: procure sempre apagar o que não é mais necessário
rm(x1, x2, y)
 
## note a diferença (e semelhancas) entre os modelos
mod   <-lm(y~ x1 + x2 + x1*x2, data=da)
mod
lm(y~ x1 + x2 + I(x1*x2), data=da)
lm(y~ x1 + x2 + x1:x2, data=da)
lm(y~ x1*x2, data=da)
lm(y~ (x1 + x2)^2, data=da)
 
 
sdmod <-sqrt(sum(residuals(mod)^2)/df.residual(mod))
## .. ou direto de summary()
sdmod <- summary(mod)$sigma
 
# Modelo, com a retirada de uma rodada de outliers...
mod  <-lm(y~ x1 + x2 + x1*x2, data=da[abs(residuals(mod))<3*sdmod,])
sdmod<-sqrt(sum(residuals(mod)^2)/df.residual(mod))
 
#Fase II
# Amostra monitorada
xs  <-vector(length=adados)
xs1 <-runif(adados, min=-1, max=1)
xs2 <-runif(adados, min=-1, max=1)
 
es  <-rnorm(adados, mean=0, sd=4)
xs  <-cbind(1, xs1, xs2, xs1*xs2)
 
res <-vector(length=adados)
j   <- 1:adados
res[j] <- 0
ypred <- predict(mod, newdata=data.frame(x1=xs1, x2=xs2))
 
for (j in 1:adados) {
 res[j]<- (coefal%*%xs[j,]+es[j])-ypred[j] 
}
 
 
## Obtendo as estimativas de beta sem usar lm() usando contas matriciais
## beta = (X'X)^(-1) X'y
## operação correta porém  ineficiente e inadequada:
solve(t(X) %*% X) %*% t(X) %*% y
## operação feita de maneira mais adequada:
solve(crossprod(X), crossprod(X,y))