Não foi possível enviar o arquivo. Será algum problema com as permissões?
Diferenças
Aqui você vê as diferenças entre duas revisões dessa página.
|
cursos:rbelem:dia3a [2010/05/27 11:59] paulojus criada |
cursos:rbelem:dia3a [2010/05/29 11:20] (atual) paulojus |
||
|---|---|---|---|
| Linha 1: | Linha 1: | ||
| + | ====== Dia 4 - a ====== | ||
| + | |||
| + | 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. | ||
| + | |||
| <code R> | <code R> | ||
| + | ## definindo número de dados | ||
| idados<- 50 | idados<- 50 | ||
| - | adados<-20000 | + | ##adados<-20000 |
| - | #Equação verdadeira de relacionamento (usada para simular Fase I e gerar Modelo de Regressão) | + | # Equação verdadeira de relacionamento |
| - | coef0 <- c(10, 2, -3, 2.3) | + | ## 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 | #Equação alterada | ||
| - | coefal<-c(10, 6, -3, 2.3) | + | #coefal<-c(10, 6, -3, 2.3) |
| #Fase I | #Fase I | ||
| + | ## simulando o valor das covariáveis X1 e X2 | ||
| x1 <-runif(idados, min=-1, max=1) | x1 <-runif(idados, min=-1, max=1) | ||
| x2 <-runif(idados, min=-1, max=1) | x2 <-runif(idados, min=-1, max=1) | ||
| - | | + | ## simulando os erros: |
| e <-rnorm(idados, mean=0, sd=4) | e <-rnorm(idados, mean=0, sd=4) | ||
| + | ## criando a matriz X | ||
| x <-cbind(1, x1, x2, x1*x2) | 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 <- model.matrix(~ x1*x2) | ||
| X | X | ||
| - | model.matrix(~ I(x1*x2)) | + | ## 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))) | ||
| Linha 26: | Linha 47: | ||
| da <-data.frame (x1, x2, x1*x2) | da <-data.frame (x1, x2, x1*x2) | ||
| + | ## colocando as covariáveis num data-frame | ||
| da <-data.frame (x1=x1, x1=x2) | 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) | y <- vector(length=idados) | ||
| ii<-1:idados | ii<-1:idados | ||
| Linha 34: | Linha 61: | ||
| } | } | ||
| - | yy <- as.vector(X %*% coef0 + e) | + | ## Com: note que não precisa fazer o loop for() pois o R é linguagem vetorial |
| - | all.equal(y, yy) | + | ## alem disto poderia criar o y direto no data-frame |
| + | da$y <- X %*% coef0 + e | ||
| da$y <-y | da$y <-y | ||
| - | da$y <- as.vector(X %*% coef0 + e) | ||
| - | rm(x1, x2, y, yy) | + | 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 + 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) | ||
| - | mod <-lm(y~ x1 + x2 + I(x1*x2), data=da) | ||
| - | mod <-lm(y~ x1 + x2 + x1:x2, data=da) | ||
| - | mod <-lm(y~ x1*x2, data=da) | ||
| sdmod <-sqrt(sum(residuals(mod)^2)/df.residual(mod)) | sdmod <-sqrt(sum(residuals(mod)^2)/df.residual(mod)) | ||
| + | ## .. ou direto de summary() | ||
| sdmod <- summary(mod)$sigma | sdmod <- summary(mod)$sigma | ||
| | | ||
| Linha 75: | Linha 107: | ||
| - | ## | + | ## 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 | solve(t(X) %*% X) %*% t(X) %*% y | ||
| + | ## operação feita de maneira mais adequada: | ||
| solve(crossprod(X), crossprod(X,y)) | solve(crossprod(X), crossprod(X,y)) | ||
| + | </code> | ||
| - | </code> | ||
