Objetos do R podem ser separados entre objetos de dados (vetores, matrizes, arrays, data-frames e listas) e outros tipos de objetos. As fórmulas constituem um tipo especial de objeto no Rque representam simbolicamente relação entre variáveis e/ou objetos de dados. Fórmulas são em geral usadas em funções gráficas e funções que analisam dados a partir de algum modelo definido pela fórmula.
Nesta seção vamos fazer uma breve introdução ao uso de fórmulas através de alguns exemplos de análises de dados. Para isto iremos utilizar o conjunto de dados mtcars disponível com o R. Este conjunto contém características técnicas de diversos modelos de automóvel. Para carregar os dados e listar os nomes das variáveis utilize os comandos a seguir. Lembre-se ainda que help(mtcars) irá fornecer mais detalhes sobre estes dados.
> data(mtcars) > names(mtcars)
[1] "mpg" "cyl" "disp" "hp" "drat" "wt" "qsec" "vs" "am" "gear" [11] "carb"
Algumas (mas não todas!) funções gráficas do R aceitam uma fórmula como argumento. Em geral tais funções exibem gráficos para explorar a relação entre variáveis. O R possui dois tipos de sistemas gráficos: (i) gráficos base (base graphics) e (ii) gráficos lattice. Os exemplos mostrados aqui se referem apenas ao primeiro sistema. Gráficos lattice são disponibilizados pelo pacote lattice. no qual as fórmulas são ainda mais importantes e largamente utilizadas.
A Figura 52 mostra dois tipos de gráficos que são definidos a partir de fórmulas. No primeiro a variável de rendimento (mpg: milhas por galão) é relacionada com uma variável categórica (cyl: número de cilindros). No segundo caso o rendimento é relacionado com o peso do veículo. A fórmula do tipo \(y \sim {} x\) pode ser lida como: a variável \(y\) é explicada por \(x\).
> with(mtcars, boxplot(mpg ~ cyl)) > with(mtcars, plot(mpg ~ cyl)) > with(mtcars, plot(mpg ~ wt))
A Figura 53 mostra agora um exemplo onde o gráfico de rendimento explicado pelo peso é feito para cada número de cilindros separadamente. Neste caso a formula usa o símbolo | para indicar condicionamento e é do tipo \(y \sim {} x|A\) podendo ser lida como a relação entre \(y\) e \(x\) para cada nível da variável \(A\).
> coplot(mpg ~ wt | as.factor(cyl), data = mtcars, + panel = panel.smooth, rows = 1)
Assim como no caso de gráficos, algumas funções de análise de dados também aceitam fórmulas em seus argumentos. Considere o exemplo do teste-\(t\) para comparação de duas amostras na comparação do rendimento de veículos com câmbio automático e manual. No exemplo a seguir mostramos o uso da função de duas formas que produzem resultados idênticos, uma sem usar fórmula e outra usando fórmula.
> with(mtcars, t.test(mpg[am==0], mpg[am==1], var.eq=T))
Two Sample t-test
data: mpg[am == 0] and mpg[am == 1]
t = -4.1061, df = 30, p-value = 0.000285
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-10.84837 -3.64151
sample estimates:
mean of x mean of y
17.14737 24.39231
> with(mtcars, t.test(mpg ~ am, var.eq=T))
Two Sample t-test
data: mpg by am
t = -4.1061, df = 30, p-value = 0.000285
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-10.84837 -3.64151
sample estimates:
mean in group 0 mean in group 1
17.14737 24.39231
Portanto em \(mpg \sim {} am\) pode-se ler: rendimento (mpg) explicado por tipo de câmbio (am). De forma similar a função para comparação de variâncias também pode utilizar fórmulas.
> with(mtcars, bartlett.test(mpg ~ am))
Bartlett test of homogeneity of variances
data: mpg by am
Bartlett's K-squared = 3.2259, df = 1, p-value = 0.07248
A fórmula é um objeto do R e possui a classe formula. Desta forma, funções que tenham métodos para esta classe tratam o objeto adequadamente. Por exemplo, no caso de t.test recebendo uma formula como argumento o método formula para t.test é disponível, como indica a documentação da função.
## S3 method for class 'formula': t.test(formula, data, subset, na.action, ...)
A seguir reforçamos estas ideias e vemos alguns comandos aplicados à manipulação de fórmulas. As funções all.vars() e terms() são particularmente úteis para manipulação de fórmulas e objetos dentro de funções.
> class(mpg ~wt)
[1] "formula"
> form1 <- mpg ~ wt > class(form1)
[1] "formula"
> all.vars(form1)
[1] "mpg" "wt"
> terms(form1)
mpg ~ wt attr(,"variables") list(mpg, wt) attr(,"factors") wt mpg 0 wt 1 attr(,"term.labels") [1] "wt" attr(,"order") [1] 1 attr(,"intercept") [1] 1 attr(,"response") [1] 1 attr(,".Environment") <environment: R_GlobalEnv>
Entre os diversos usos de fórmulas, o mais importante deles é sem dúvida o fato que fórmulas são utilizadas na declaração de modelos estatísticos. Um aspecto particularmente importante da linguagem S, e portanto no programa R, é que adota-se uma abordagem unificada para modelagem, o que inclui a sintaxe para especificação de modelos. Variáveis respostas e covariáveis (variáveis explanatórias) são sempre especificadas de mesma forma básica, ou seja, na forma \(resposta \sim {} covariavel\), onde:
No restante deste texto vamos, por simplicidade, considerar que há apenas uma variável resposta que poderá ser explicada por uma ou mais covariáveis.
Considere, para o conjunto de dados mtcars, ajustar um modelo que explique o rendimento (Y:mpg) pelo peso do veículo (X:wt). O modelo linear é dado por: \[Y = \beta _0 + \beta _1 X + \epsilon \; ,\] e pode ser ajustado no R usando lm() (lm : linear model). Na sintaxe da chamada função \(mpg \sim {} wt\) lê-se: mpg é modelado por wt, através de um modelo linear lm(), o que implica no modelo acima. A Figura 54 mostra os dados e a linha sólida mostra a equação do modelo ajustado.
> reg1 <- lm(mpg ~ wt, data=mtcars) > reg1
Call: lm(formula = mpg ~ wt, data = mtcars) Coefficients: (Intercept) wt 37.285 -5.344
Note que a fórmula apenas especifica a relação entre as variáveis resposta e explanatórias e não implica que o modelo seja necessariamente linear. A linearidade é dada pela função lm(). Portanto a mesma fórmula pode ser usada para outros tipos de ajuste como o mostrado na linha tracejada do gráfico, resultante de regressão local polinomial obtida por loess().
Nem todas as funções que relacionam variáveis aceitam formulas, como por exemplo o caso da regressão por núcleo (kernel) dada por ksmooth() cujo ajuste é mostrado na linha pontilhada. Outras funções estendem a notação de funções como é o caso do ajuste por modelos aditivos generalizados gam() mostrado na linha sólida grossa, onde o termo s() indica que a variável resposta deve ser descrita por uma função suave da covariável incluída neste termo.
> with(mtcars, plot(mpg ~ wt)) > abline(reg1) > reg2 <- loess(mpg ~ wt, data=mtcars) > wts <- with(mtcars, seq(min(wt), max(wt), len=201)) > lines(wts, predict(reg2, data.frame(wt=wts)), lty=2) > lines(with(mtcars, ksmooth(wt,mpg,band=1)), lty=3) > require(mgcv) > reg3 <- gam(mpg ~ s(wt), data=mtcars) > lines(wts, predict(reg3, data.frame(wt=wts)), lwd=2) > legend("topright", c("lm()","loess()","ksmooth()", "gam()"), lty=c(1:3,1), lwd=c(1,1,1,2))
Nos exemplos acima é interessante notar o uso de predict() que é utilizada para predizer o valor da resposta para um conjunto arbitrário de valores da covariável, baseando-se no modelo ajustado. No exemplo utilizamos este recurso para produzir o gráfico com a "curva"do modelo ajustado para uma sequência de valores da covariável. Para a função lm() utilizamos apenas abline() devido ao fato que esta função retorna a equação de uma reta que é interpretada a traçada por um método abline. Entretanto predict() também poderia ser usada e a reta traçada com o comando a seguir. Esta forma é mais flexível para traçar funções (modelos) ajustados que sejam mais complexos que uma equação de uma reta.
> lines(wts, predict(reg1, data.frame(wt=wts)))
As formulas admitem operadores aritméticos em seus termos. Por exemplo considere a relação entre o rendimento (mgp) e a potência (hp). A linha sólida no gráfico da esquerda da Figura 55 sugere que o modelo linear não descreve bem a relação entre estas variáveis, enquanto no gráfico da direita sugere a relação é melhor descrita por um modelo linear entre o rendimento e o logaritmo de potência. Na chamada das funções utilizamos a operação aritmética log() diretamente na fórmula, sem a necessidade de transformar os dados originais.
> with(mtcars, plot(mpg ~ hp)) > abline(lm(mpg ~ hp, data=mtcars)) > with(mtcars, plot(mpg ~ log(hp))) > abline(lm(mpg ~ log(hp), data=mtcars))
Uma outra possibilidade para os dados originais é o ajuste de um modelo dado por uma função
polinomial, conforme mostrado na linha tracejada no gráfico da esquerda da Figura 55 e que é
ajustado com os comandos a seguir. Neste ajuste é importante notar que a variável quadrática deve
ser especificada com I(hp^2) e o uso de I() é obrigatório para garantir que os valores de hp sejam de
fato elevados ao quadrado. O uso de hp^2 possui um significado diferente que veremos na próxima
sessão.
> polA <- lm(mpg ~ hp + I(hp^2), data=mtcars) > hps <- seq(50,350,len=200) > lines(hps, predict(polA, data.frame(hp=hps)), lty=2)
Uma outra forma de especificar regressões polinomiais é com o uso de poly(), onde o grau do desejado do polinômio é um argumento desta função conforme ilustrado nos comandos a seguir. No exemplo é importante notar que a interpretação dos parâmetros é diferente devido ao fato de que polinômios ortogonais são utilizados. Entretanto os valores preditos e as estatísticas de ajuste são iguais. O ajuste por polinômios ortogonais é numericamente mais estável e portanto deve ser preferido quando possível. Quando se usa as opções default a função poly() vai sempre contruir polinômios ortogonais. Caso queira-se usar potências usuais deve-se adicionar à chamada desta função o argumento raw=T.
> polA
Call: lm(formula = mpg ~ hp + I(hp^2), data = mtcars) Coefficients: (Intercept) hp I(hp^2) 40.4091172 -0.2133083 0.0004208
> polB <- lm(mpg ~ poly(hp,2), data=mtcars) > polB
Call: lm(formula = mpg ~ poly(hp, 2), data = mtcars) Coefficients: (Intercept) poly(hp, 2)1 poly(hp, 2)2 20.09 -26.05 13.15
> hps <- seq(50,350,by=50) > predict(polA, data.frame(hp=hps))
1 2 3 4 5 6 7 30.79574 23.28645 17.88123 14.58009 13.38303 14.29005 17.30114
> predict(polB, data.frame(hp=hps))
1 2 3 4 5 6 7 30.79574 23.28645 17.88123 14.58009 13.38303 14.29005 17.30114
Vamos considerar agora um outro exemplo de ajuste de modelo linear, agora para o conjunto de dados women que fornece peso (weight) em libras (lbs) e altura (height) em polegadas (in) de 15 mulheres americanas de 30-39 anos. Os comandos a seguir mostram os quatro ajustes indicados na Figura 56. O primeiro (linha fina sólida) é uma regressão linear, o segundo (linha fina tracejada) é uma regressão linear com intercepto igual a zero, isto é, a reta passa pela origem. O terceiro (linha sólida grossa) é uma regressão quadrática e o quarto (linha sólida grossa) é uma regressão quadrática passando pela origem. Neste exemplo fica então ilustrado que a adição do termo + 0 na fórmula faz com que o intercepto do modelo seja nulo e apenas o parâmetro referente ao coeficiente angular da reta seja estimado.
> data(women) > wm1 <- lm(weight ~ height, data=women) > wm2 <- lm(weight ~ height+0, data=women) > wm3 <- lm(weight ~ height + I(height^2), data=women) > wm4 <- lm(weight ~ height + I(height^2) + 0, data=women) > with(women, plot(weight ~ height)) > hgs <- seq(58,72,l=200) > lines(hgs, predict(wm1, data.frame(height=hgs))) > lines(hgs, predict(wm2, data.frame(height=hgs)), lty=2) > lines(hgs, predict(wm3, data.frame(height=hgs)), lwd=2) > lines(hgs, predict(wm4, data.frame(height=hgs)), lty=2, lwd=2) > legend("topleft", c("linear", "linear s/ intercepto", "quadrático", "quadrático s/ intercepto"), lty=c(1,2,1,2), lwd=c(1,1,2,2), cex=0.85)
Nos exemplos anteriores a variável resposta era explicada por apenas uma variável explanatória. Isto pode ser expandido considerando-se a presença de duas ou mais variáveis explicativas. A Tabela 5 resume as principais operações possíveis para definir modelos com uma ou duas variáveis e que podem ser estendidas para o caso de mais variáveis.
| Termos | Especificação |
A + B | Efeitos principais \(A\) e \(B\) |
A : B | Termo de interação entre \(A\) e \(B\) |
A * B | Efeitos principais e interação, corresponde a \(A + B + A:B\) |
B %in% A | \(B\) dentro (aninhado) de \(A\) |
A/B | Efeito principal e aninhado, corresponde a \(A + B\%in\%A\) |
A-B | tudo de \(A\) exceto o que está em \(B\) |
A^k | Todos os termos de \(A\) e interação de ordem \(k\) |
A + 0 | exclui o intercepto de modelo |
I() | operador de identidade aritmética, ver explicação no texto |
Esta notação é uma implementação das ideias propostas por Wilkinson e Rogers para especificação de modelos estatísticos. (G. N. Wilkinson. C. E. Rogers. Symbolic Description of Factorial Models for Analysis of Variance. Applied Statistics, Vol. 22, No. 3, 392-399. 1973).
Para ilustrar algumas destas opções vamos considerar novamente o conjunto de dados mtcars ajustando modelos para o rendimento (mpg) explicado pelo peso (wt) e potência (hp) do veículo. Nos comandos a seguir mostramos os coeficientes estimados a partir de cinco formas de especificação de modelos.
> coef(lm(mpg ~ I(wt + hp), data=mtcars))
(Intercept) I(wt + hp) 30.2877307 -0.0680239
> coef(lm(mpg ~ wt + hp, data=mtcars))
(Intercept) wt hp 37.22727012 -3.87783074 -0.03177295
> coef(lm(mpg ~ I(wt * hp), data=mtcars))
(Intercept) I(wt * hp) 27.74564216 -0.01487156
> coef(lm(mpg ~ wt * hp, data=mtcars))
(Intercept) wt hp wt:hp 49.80842343 -8.21662430 -0.12010209 0.02784815
> coef(lm(mpg ~ (wt + hp)^2, data=mtcars))
(Intercept) wt hp wt:hp 49.80842343 -8.21662430 -0.12010209 0.02784815
> coef(lm(mpg ~ I((wt + hp)^2), data=mtcars))
(Intercept) I((wt + hp)^2) 24.4985252043 -0.0001625815
Os resultados sugerem que as fórmulas definem modelos diferentes, exceto pelos termos wt * hp e (wt * hp)^2 onde o mesmo modelo é especificado de duas formas alternativas. Os modelos ajustados para explicar o rendimento mpg denotado por \(Y\) são:
Chama-se atenção ao fato que a notação de "potência" em (wt+hp)^2 não indica uma operação
aritmética mas sim a inclusão de todos os efeitos principais e interações até as de ordem indicada pela
potência. Para incluir a operação aritmética de potência é necessário utilizar I() no termo a ser
exponenciado.
De forma geral, a mensagem é de que os operadores soma (+), produto (*), divisão (/) e potência
(
) têm nas fórmulas o papel de definir numa notação simbólica quais e como os termos devem
ser incluídos no modelo. Em fórmulas, tais operadores só indicam operações aritméticas
com os termos envolvidos quando utilizados dentro de I(). A função I() garante que a
expressão nela contida seja avaliada como uma função aritmética, tal e qual está escrita ("as
is").
Na tabela 6 são ilustradas mais algumas especificações de modelos. No caso marcado com (*) os modelos são equivalentes porém os coeficientes resultantes são diferentes como comentado sobre polinômios ortogonais na Sessão 23.5.
| Declaração | Modelo equivalente | Descrição |
A+B*C | A+B+C+B:C | todos efeitos principais e interação dupla apenas |
| entre \(B\) e \(C\) | ||
A+B*(C+D) | A+B+C+D+B:C+B:D | todos efeitos principais e interações duplas |
| de \(B\) com \(C\) e \(B\) com \(D\) | ||
A*B*C | A+B+C+A:B+A:C+B:C+A:B:C | todos efeitos principais e interações possíveis |
(A+B+C)^3 | A+B+C+A:B+A:C+B:C+A:B:C | três covariáveis e interações de ordem 2 e 3 |
| (igual ao anterior) | ||
(A+B+C)^2 | A+B+C+A:B+A:C+B:C | três covariáveis e interações de ordem 2 |
(A+B+C)^3 - A:B:C | (A+B+C)^2 | três covariáveis e interações de ordem 2 |
(A+B+C)^2 - A:C | A+B+C+A:B+B:C | três covariáveis e interações de ordem 2, |
| exceto por \(A:C\) | ||
A+I(A^2)+I(A^3) | poly(A,3) | regressão polinomial cúbica em \(A\) (*) |
A+I(A^2)+I(A^3) | poly(A,3,raw=TRUE) | regressão polinomial cúbica em \(A\) |
A+I(A^2)+B | poly(A,2)+B | termos lineares em \(A\) e \(B\) e quadrático em \(A\) (*) |
A+I(A^2)+B | poly(A,2,raw=TRUE)+B | termos lineares em \(A\) e \(B\) e quadrático em \(A\) |
y ~ . | A+B+... | inclui como covariáveis todas as variáveis no |
| objeto de dados, exceto a resposta | ||
y ~ . - A | B+... | inclui como covariáveis todas as variáveis no |
| objeto de dados, exceto a resposta e a variável \(A\) | ||
Uma vez que um objeto contenha uma fórmula, é possível obter uma nova fórmula que seja uma modificação da original utilizando update.formula().
> form1 <- y~x1+x2+x3 > form1
y ~ x1 + x2 + x3
> form2 <- update.formula(form1, . ~ . - x2) > form2
y ~ x1 + x3
A lógica da sintaxe é que o primeiro argumento recebe uma fórmula inicial e o segundo indica a
modificação. O caracter ponto (\(\cdot \)) indica tudo. Ou seja, em . ~ . - x2 entende-se: a nova fórmula
deverá possuir tudo que estava do lado esquerdo, e tudo do lado direito, excluindo a variável
x2.
Este mecanismos é útil para evitar que fórmulas precisem ser totalmente redigitadas a cada redefinição do modelo, o que é útil ao se investigar vários modelos que são obtidos uns a partir de outros. O mecanismo também reduz a chance de erros nas especificações uma vez que garante a igualdade daquilo que é indicado pela notação de ponto (\(\cdot \)).