© Ricardo Solar/UFMG - compartilhamento e utilização não-comercial livres. Não reproduzir sem autorização >http://doi.org/10.5281/zenodo.7392285
Agora, começamremos a fazer as análises de maneira mais direta no R. Vamos utilizar uma tabela que contém dados de altura e peso de indivíduos. Os dados foram coletados em turmas de R em vários anos, de turmas de Graduação e Pós-Graduação. Vamos primeiro importar os dados, detalhando como isso deve ser feito.
Ainda que saibamos que tudo que faremos neste curso é importante, essa parte merece uma atenção especial! Se você não conseguir importar seus dados de maneira correta e direta ao R, de muito pouca utilidade o programa será pra você! O carregamento de dados pode parecer demasiado técnico do ponto de vista computacional, e por isso farei algumas recomendações que podem facilitar esse processo. Aqui é importante destacar que tais recomendações não são as únicas formas de carregar seus dados, mas padronizar de alguma maneira essa tarefa pode diminuir suas dores de cabeça.
De maneira simplificada, importar seus dados no R envolve a seguinte sequência:
| ID | Var_1 | Var_2 | Var 3 |
|---|---|---|---|
| amostra_1 | valor_v1 | valor_v2 | valor_v3 |
| amostra_2 | valor_v1 | valor_v2 | valor_v3 |
| … | … | … | … |
| amostra_n | valor_v1 | valor_v2 | valor_v3 |
R em qual pasta o arquivo de trabalho está localizado - funções getwd e setwd que veremos abaixo.R - funções de importação de dados (read.table, read.csv, read.excel, etc)Environment e funções summary, str, plot, View, etc.Os dados podem ser importados para o R em vários formatos, basicamente todas as formas de texto existentes. Exatamente por esse motivo, recomendo que você escolha um formato preferencial que use sempre, evitando se confundir. É importante, todavia, ter uma noção dos principais formatos de dados que trabalha-se no R. Vamos ver alguns dos mais comuns:
.txt - o formato TXT é bastante versátil. É um formato de texto simples, onde os valores podem ser separados por espaços ou tabulações. É um arquivo de texto puro, isto é, sem formatação ou fórmulas. Isso é muito útil, pois esse tipo de arquivo não será importado para o R caso haja qualquer erro remanescente. Para uso no R, recomendo fortemente exportar o arquivo do Excel (ou seus equivalentes) como texto (separado por tabulações). A importação desse tipo de arquivo no R é bem direta, através da função read.table..csv - o formato é uma abreviação para comma separated values (de valores separados por vírgula). Também é um arquivo de texto puro, como o TXT. Há alguns detalhes, se arquivo CSV for exportado de uma versão que utiliza formatação em inglês do Excel, os valores em cada linha serão separados por vírgulas (,). Se o arquivo CSV for gerado usando uma versão que usa formatação dinamarquesa do Excel, a vírgula é usada para casas decimais e os valores são separados por ponto e vírgula (;). Essa é uma das questões técnicas que mencionei e que podem enrolar o uso de alguns arquivos. A importação desse tipo de arquivo no R é bem direta, através da função read.csv ou read.csv2.xls ou xlsx - é possível ainda importar os dados diretamente de um arquivo Excel. O R não faz essa importação direta e você precisará lançar mão de um pacote (veremos em breve) para importar esses arquivos (por exemplo o pacote readxl, com sua função read.excel). A importação direta do arquivo Excel é bem interessante, mas o arquivo vem carregado de formatações, fórmulas e adjuntos trazidos pelo Excel. É importante tomar cuidado na conferência dos dados.R?Há informações que você precisa saber antes de importar qualquer dado no R, que baseiam-se majoritariamente em como a tabela foi construída, vou listar alguns abaixo, mas vale a pena ir construindo seu formato de fazer a tabela para evitar problemas.
NA, tanto no formato, quanto maiúsculas?R. O padrão do programa é separar decimais por pontos, então se usou pontos, não precisa mais se preocupar.stringsAsFactors = T, uma vez que o padrão do programa é que eles sejam lidos como texto.Agora que sabemos tudo isso, e estamos certos do nosso arquivo, posso continuar a importar meus dados.
Para trabalhar com os dados no R, precisamos informar qual é o nosso diretório de trabalho. A forma mais simples de fazer isso, como já vimos, é criar uma pasta para cada análise que faremos e trabalhar dentro daquela pasta. Essa prática apresenta uma série de vantagens. A maior delas é saber que todos os arquivos que preciso importar estão na mesma pasta. Outra grande vantagem é que todo arquivo que eu exportar, será gravado naquela pasta também.
setwdNessa etapa, é importante ter alguma intimidade com seu computador. Isso vai fazer sua vida mais fácil.
Há alguns caminhos possíveis, mas todos eles devem em última instância lhe gerar um caminho no computador que possa ser colocado no script, para que em rodadas futuras daquela análise isso já esteja resolvido.
Estes dois comandos irão colocar num console um comando começado com setwd similar ao que apareceu aqui no endereço do meu computador:
setwd("~/Dropbox/UFMG/Disciplinas/R UFAC/Aulas/teste t e anova"). Veja que este endereço é exclusivo desta pasta no meu computador, haverá um diferente no seu. Copie este endereço do console e cole ele no seu script, para que em sua próxima rodada dessa análise não haja necessidade de fazer isso de novo.
Vamos agora importar nosso arquivo, que será utilizado nas análises a seguir. Conferindo:
*.txt;Com isso em mente, vamos lá.
Primeira parte, mudar o diretório:
setwd("~/Dropbox/UFMG/Disciplinas/R UFAC/Aulas/teste t e anova")
##caso eu queira conferir, posso usar os comandos
getwd() #para que o programa me mostre o diretório que ele está setado
## [1] "/Users/ricardosolar/Dropbox/UFMG/Disciplinas/R UFAC/Aulas/teste t e anova"
##e também
dir() ##para me mostrar os arquivos que estão na pasta que estou trabalhando.
## [1] "dados_qualif.txt"
## [2] "Herbivoria.txt"
## [3] "Regressão Linear simples e ANCOVA.Rmd"
## [4] "Regressão-Linear-simples-e-ANCOVA.html"
## [5] "Regressão-Linear-simples-e-ANCOVA.log"
## [6] "Regressão-Linear-simples-e-ANCOVA.tex"
## [7] "rsconnect"
## [8] "RTurmas.txt"
## [9] "Teste T e F - rotina usual.Rmd"
## [10] "Teste-T-e-F---rotina-urual.html"
## [11] "Teste-T-e-F---rotina-usual.html"
Agora posso importar os dados, que como vimos, estão com o nome RTurmas.txt. Como os dados estão em txt, vou usar o comando read.table, e detalhar os passos que usamos na função abaixo:
read.table( #comando base
"RTurmas.txt", #arquivo a ser importado
header = TRUE, #Informa que há um cabeçalho, ou seja, a primeira linha deve ser considerada como nomes dos dados, e não deve ser entendida como parte em si dos dados
stringsAsFactors = TRUE, #informo que, caso haja texto que se repita nas variáveis, devem ser considerados como fatores
dec = "." #como havia dito, no nosso caso esse comando é dispensável, mas deixo aqui registrado que, caso sua decimal seja ",", é aqui que deve ser informada.
)
O problema desse comando anterior é que ele só leu a planilha!! Para armazenar e utilizar no R, vale lembrar das lições anteriores … devemos designá-lo a um objeto, é bem simples, e fica assim no uso corriqueiro
dados <- read.table("RTurmas.txt", header = TRUE, stringsAsFactors = TRUE)
E pra poder conferir os dados, temos várias opções, vamos explorar algumas:
summary(x) - dá um resumo geral dos dados, focando nos valores e métricas das variáveissummary(dados)
## ID_Ind Altura Peso Sexo Turma
## Min. : 1.00 Min. :150.0 Min. : 44.50 F:45 Grad:45
## 1st Qu.:23.75 1st Qu.:163.0 1st Qu.: 56.77 M:47 Pos :47
## Median :46.50 Median :168.0 Median : 66.02
## Mean :46.50 Mean :170.1 Mean : 68.33
## 3rd Qu.:69.25 3rd Qu.:177.2 3rd Qu.: 76.03
## Max. :92.00 Max. :195.0 Max. :113.70
## Ano
## Min. :2015
## 1st Qu.:2016
## Median :2017
## Mean :2017
## 3rd Qu.:2017
## Max. :2019
str(x) - resume a estrutura das variáveisstr(dados)
## 'data.frame': 92 obs. of 6 variables:
## $ ID_Ind: int 1 2 3 4 5 6 7 8 9 10 ...
## $ Altura: int 162 160 165 165 164 165 158 160 166 165 ...
## $ Peso : num 44.5 46 52 53 54.2 56.1 57.6 58 62 79.1 ...
## $ Sexo : Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
## $ Turma : Factor w/ 2 levels "Grad","Pos": 1 1 1 1 2 2 2 2 2 2 ...
## $ Ano : int 2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
dim(x) - mostra os números de linhas e colunas dos dadosdim(dados)
## [1] 92 6
names(x) ou colnames(x) - mostra os nomes das colunas dos dados, muito útil para ver se isso foi importado corretamentenames(dados)
## [1] "ID_Ind" "Altura" "Peso" "Sexo" "Turma" "Ano"
colnames(dados)
## [1] "ID_Ind" "Altura" "Peso" "Sexo" "Turma" "Ano"
head(x) e tail(x) - mostra o início (head) ou o final (tail) dos dados, como padrão tendo as 6 primeiras ou 6 últimas linhas.#"cabeça" dos dados
head(dados)
#"cauda" dos dados
tail(dados)
plot(x) - por fim, posso fazer um plot geral dos dados, pra entender como estão distribuídos e, visualmente, identificar algum problemaplot(dados)
Com os dados conferidos e corretos, posso começar a trabalhar os mesmos.
Queremos inicialmente verificar se a altura dos estudantes que cursaram a disciplina diferem entre os sexos biológicos. Isso pode ser feito através de um teste T.
A primeira premissa é que a variável resposta (Altura) deve ter distibuição normal. Como verificar isso? Uma das formas é o teste de Shapiro, que é implementado no R através de shapiro.test(x). Quero testar a variável Altura, então preciso informar ao programa naquele formato dados$Altura. O teste de shapiro tem como hipótese nula que a variável em questão possui distribuição Normal.
shapiro.test(dados$Altura)
##
## Shapiro-Wilk normality test
##
## data: dados$Altura
## W = 0.97605, p-value = 0.08814
O teste de shapiro não foi significativo, ou seja, não rejeito a hipótese nula. Portanto, minha variável possui distribuição Normal. Posso prosseguir com o teste T. Como vimos antes, o teste T é implementado através do comando t.test(y~x). Que irei armazenar num objeto chamado testeT.
(testeT <- t.test(Altura~Sexo, data=dados)) ## colocar todo o comando entre parêntesis permite que o R crie o objeto E mostre seu resultado
##
## Welch Two Sample t-test
##
## data: Altura by Sexo
## t = -9.6706, df = 89.763, p-value = 1.418e-15
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -16.03837 -10.57155
## sample estimates:
## mean in group F mean in group M
## 163.3333 176.6383
Veja que terminei o comando com data=dados. Fiz isso para informar ao programa que os termos Altura e Sexo estão dentro do dataframe dados.
Quando fiz o teste, gravei seus resultados no objeto testeT para poder extair dele os resultados sem ter que copiar e colar do console. Isso será principalmente útil para os gráficos. Vamos começar bem simples, fazendo um gráfico contendo o valor das médias e o intervalo de confiança. Para isso, já utilizarei o pacote GGPLOT2. Ainda não explicarei a ideia de semântica de gráficos, que virá logo em breve.
Vamos ver o que está contido dentro do objeto - usaremos o comando str(x)
str(testeT)
## List of 10
## $ statistic : Named num -9.67
## ..- attr(*, "names")= chr "t"
## $ parameter : Named num 89.8
## ..- attr(*, "names")= chr "df"
## $ p.value : num 1.42e-15
## $ conf.int : num [1:2] -16 -10.6
## ..- attr(*, "conf.level")= num 0.95
## $ estimate : Named num [1:2] 163 177
## ..- attr(*, "names")= chr [1:2] "mean in group F" "mean in group M"
## $ null.value : Named num 0
## ..- attr(*, "names")= chr "difference in means between group F and group M"
## $ stderr : num 1.38
## $ alternative: chr "two.sided"
## $ method : chr "Welch Two Sample t-test"
## $ data.name : chr "Altura by Sexo"
## - attr(*, "class")= chr "htest"
Ok, temos que o que nos interessa nessa estrtura são os pontos $estimate, onde estão armazenadas as médias e $stderr, que intuitivamente sabemos que armazena o erro padrão. Para fazer o gráfico, preciso montar isso num dataframe que permita ao programa entender os valores da forma correta. É bem simples, vamos lá:
dados.plot <- data.frame(medias = testeT$estimate, intervalo=abs(testeT$stderr), sexo=c("F", "M")) ##este último eu deduzi pela ordem alfabética, que o R sempre segue como padrão, F vem antes de M.
Veja acima que pedi a criação de um objeto data.frame, com o nome de dados.plot. Nele coloquei os elementos contidos em testeT$estimate e os elementos em testeT$stderr, além de criar uma outra coluna que continha os nomes c("F", "M"). Para incluir os pontos reais que deram origem às médias, devo usar o geom_jitter, mas preciso informar quque os dados estão em outro dataframe, confira abaixo. Veja que para o intervalo de confiança, me interessam somente os valores em módulo, ou valores absolutos.
library(ggplot2)
ggplot(dados.plot, aes(x=sexo, y=medias))+
geom_jitter(data=dados, aes(x=Sexo, y=Altura), pch=18)+
geom_point(pch=21, size=4, fill="grey50")+
geom_errorbar(aes(ymin=medias-intervalo, ymax=medias+intervalo), position = position_dodge(), width=.3)+
xlab("Sexo biológico")+
ylab("Altura (cm)")
Assim então temos o teste T simples, sempre entre uma variável resposta e uma variável explicativa categórica (duas categorias).
Para o teste T pareado, nosso interesse é diferente. Queremos saber **como a mesma amostra se comporta em dois momentos distintos, ou seja, quero acompanhar essas variáveis em dois momentos.
Para exemplificar o uso desse tipo de teste, vamos trabalhar com dados de um exemplo em que as informações foram coletadas antes e depois de plantas serem expostas à herbivoria. A hipótese é que a quantidade de compostos secundários das folhas aumenta após a exposição à herbivoria.
dadosherb <- read.table("Herbivoria.txt", h=T, stringsAsFactors = T)
summary(dadosherb)
## Planta Tempo Concentracao
## Min. : 1.0 after :10 Min. : 8.70
## 1st Qu.: 3.0 before:10 1st Qu.:11.97
## Median : 5.5 Median :13.45
## Mean : 5.5 Mean :13.06
## 3rd Qu.: 8.0 3rd Qu.:14.30
## Max. :10.0 Max. :16.50
A nossa pergunta então é bem importante para definir o modelo. Meu interesse aqui é saber se a concentração de compostos aumenta depois do tratamento. Portanto, nosso modelo é a \(Concentração ~ Tempo\), levando em conta que eu já espero que os valores sejam maiores depois e que há valores pareados para a mesma planta antes e depois.
shapiro.test(dadosherb$Concentracao)
##
## Shapiro-Wilk normality test
##
## data: dadosherb$Concentracao
## W = 0.96403, p-value = 0.6271
(testeTpar <- t.test(Concentracao~Tempo, data=dadosherb, paired=TRUE, alternative="greater"))
##
## Paired t-test
##
## data: Concentracao by Tempo
## t = 2.272, df = 9, p-value = 0.0246
## alternative hypothesis: true mean difference is greater than 0
## 95 percent confidence interval:
## 0.1043169 Inf
## sample estimates:
## mean difference
## 0.54
Nossos resultados, primeiro do teste de Shapiro, mostram que posso usar o teste T, pois a distribuição da variável é Normal. O restultado do teste T pareado, com a hipótese alternativa de que depois do tratamento os valores são maiores, é significativa (p=0.025).
Os gráficos de testes T pareados são um pouco diferentes, pois nos interessa ver o comportamento de cada ponto. Para isso, vamos usar um pacote que permite essa plotagem, tendo como base o pacote ggplot2.
Vou aproveitar esse momento para mostrar como filtar dados, usando o pacote library(dplyr). Observem que para fazer os gráfico, preciso separar os dados, no caso entre before e after. O comando filter permite filtrar um dataframe de acordo com uma condição determinada. No nosso caso, observe que filtrei before e depois after. Como só era do meu interesse a coluna `$Concentracao’, adicionei essa indexação ao final do comando.
# Plot paired data
library(PairedData)
before <- dplyr::filter(dadosherb, Tempo=="before")$Concentracao
after <- dplyr::filter(dadosherb, Tempo=="after")$Concentracao
pd <- paired(before, after)
plot(pd, type = "prof")+
geom_line(col="coral", linewidth=1)+
#geom_jitter(width=.1)
geom_point(col="#801FA4")
E se nosso interesse fosse testar mais de uma variável, ou variáveis com mais de dois níveis? Nessa situação, podemos lançar mão da Análise de Variância(ANOVA), com o teste F. De fato, veremos que é sempre possível usar esse teste, pois quando entrarmos em GLM, esse tipo de análise permite acomodar modelos não-normais, o que nos ajuda muito!
Uma coisa importante aqui é que inicialmente também temos a premissa de normalidade das populações de onde tiramos os dados. Poderíamos também seguir com o teste de Shapiro para checar a normalidade, mas temos algumas observações a fazer. Em geral, uma ANOVA de um teste bastante robusto contra violações da premissa de normalidade, desde que os tamanhos de amostra sejam suficientemente grandes.
Além disso, se você tiver tamanhos de amostra extremamente grandes, testes estatísticos como o teste de Shapiro-Wilk quase sempre dirão que seus dados não são normais. Por esse motivo, geralmente é melhor inspecionar seus dados visualmente, através de análise dos resíduos. Simplesmente olhando para os gráficos, você pode ter uma boa ideia se os dados são ou não normalmente distribuídos. Faremos isso consistentemente de agora em diante.
Agora que já sabemos isso tudo, podemos testar nossa nova pergunta de interesse com aquele primeiro conjunto de dados, que contém as informações das turmas de R.
Nossa pergunta agora é se o Peso da turma é função do Sexo biológico e do fato de estas pessoas estarem na Graduação ou na Pós. Nosso modelo agora passa a ser então Peso~Sexo+Turma. Para isso, usaremos o comando aov, e colocaremos esse modelo sob prova. Para testar as premissas do modelo, usaremos a análise de resíduos, implementada através do pacote library(DHARMa).
modelo <- aov(Peso~Sexo+Turma, data=dados)
library(DHARMa)
simulateResiduals(modelo, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.204 0.252 0.48 0.524 0.132 0.152 0.2 0.212 0.356 0.904 0.34 0.448 0.488 0.62 0.68 0.244 0.224 0.316 0.884 0.912 ...
anova(modelo, test="F")
Portanto, ambas as variáveis são singificativas. Como ambas possuem apenas dois níveis, podemos parar por aqui e fazer o gráfico. Novamente, para fazer o gráfico, vou criar um dataframe com as informações de interesse, médias e erro padrão, além da identidade das variáveis. Aqui vou lançar mão de uma função muito interessante contida no pacote library(Rmisc), chamada summarySE. Nessa função preciso sempre informar quem são os dados, quem é a variável a ser medida (y) e quem são as variáveis agrupadoras (x).
library(Rmisc)
## Loading required package: plyr
(dados.plot <- summarySE(dados, measurevar = "Peso", groupvars = c("Sexo", "Turma")))
Vejam que tenho os dados agora, e posso prosseguir com o ggplot2. A explicação de como incorporar os dados será feita em sala.
ggplot(dados, aes(y=Peso, x=Turma, group=Sexo))+
geom_jitter(aes(y=Peso, x=Turma, fill=Sexo), position = position_dodge(.9), pch=21)+
geom_point(data= dados.plot, position = position_dodge(.9), size=3, pch=21, col="black", fill="grey50")+
geom_errorbar(data= dados.plot, aes(ymax=Peso+se, ymin=Peso-se), position = position_dodge(.9), width=.3, col="grey10")
Pronto, finalizamos esse módulo com essas duas análises bastante utilizadas. Daqui por diante no curso, seguiremos compreendendo melhor o uso dos Modelos Lineares Generalizados.