apostila_R2_André Oliveira Souza

Sumario

1 Leitura de dados usando copiar colar. 2

2 Recursos graficos 3

3 Filtro de banco de dados. 12

4 Editor LaTeX + R + Sweave no linux 134.1 Sweave linux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144.2 knitr no linux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144.3 Salvar imagens linux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144.4 Sites importantes de apoio ao R na internet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144.5 A mudanca de diretorio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144.6 Configuracoes basicas para Linux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154.7 Ajuda e demostracoes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154.8 Interface amigavel para R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154.9 Adicionando repositorios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154.10 Gerador de tabela em LaTeX e HTML . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154.11 Manipulacao de diretorios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164.12 Box plot com media destacadas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

5 Funcao densidade probabilidade fdp 18

6 Matematica 19

7 Tabela da ANOVA 197.1 Distribuicao Binomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 217.2 BoxPlot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

8 Multivariada 22

9 Modelos lineares generalizados - GLM 239.1 Redefinicao de nveis de referencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 259.2 Contrates em GLM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 259.3 Funcao de maxima verossimilhanca penalizada . . . . . . . . . . . . . . . . . . . . . . . . . . . . 259.4 Teste estatsticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 259.5 Estatstica Descritivas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 279.6 Biblioteca ExpDes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

10 Ler funcao externa 4110.1 Dados Multinacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4210.2 Grafico Lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4710.3 ScottKnott, o TukeyC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6010.4 TukeyC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6410.5 UsingR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6610.6 Regressao linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7510.7 Experimentacao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7910.8 Estatstica nao parametrica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8310.9 Regressao logistica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8410.10Sturges . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

1

Notas sobre o R Linux

Andre Oliveira Souza

27 de dezembro de 2013

1 Leitura de dados usando copiar colar.

Neste topico ilustrei tres formas de um arquivo para dentro do R apenas com copiar e colar. Depois deselecionar seus dados Ctrl+C em seu editor use read.table dentro do R.

1. dados

2 3 3 39.45

dados2

> rf h hist(rf, prob = TRUE, main="",xlab="Classes",ylab="Proporc~ao",col="limegreen",labels = sprintf("%s%%",round(prop.table(table(cut(rf, h$breaks)))*100,2)))

Classes

Prop

or

o

0 200 400 600 800 1000

0.00

000.

0005

0.00

100.

0015

13.33% 26.67%

40%

6.67%

13.33%

Figura 1: Porcentagem no histograma

> f=function(x){ 17.37011+ 40.11928*(1-exp(-0.02134*x))}

> f1=function(x){ 23.06253+ 38.21512*(1-exp(-0.01981*x))}

> curve(f, 0,96, xlab='Tempo (horas)', ylab='Degradabilidade', xlim=c(0,120), ylim=c(0,100), bty='n', col=1, lty=1)> curve(f1, 0,96, add=T, lty=2, col=2)

> legend(20,90, c('Tratamento 1', 'Tratamento 2'), col=c(1,2), lty=c(1,2), bty='n')

0 20 40 60 80 100 120

020

4060

8010

0

Tempo (horas)

Deg

rada

bilid

ade

Tratamento 1Tratamento 2

Figura 2: Adicionar curvas

4

> plot(density(rf),xlab="x",ylab="Densidade",main="",col=2,lty=2)

0 500 1000

0.00

000.

0005

0.00

100.

0015

x

Den

sida

de

Figura 3: Plotar densidade

> hist(rnorm(500,mean=2,sd=sqrt(4)),freq=FALSE,main="X~N(2,4)",xlab="",ylab="Densidade")

> curve(dnorm(x,mean=2,sd=sqrt(4)),col=2,lty=2,lwd=2,add=TRUE)

X~N(2,4)

Den

sida

de

4 2 0 2 4 6 8

0.00

0.05

0.10

0.15

0.20

Figura 4: Curva normal

5

> par(mfrow=c(2,1))

> x hist(x,prob=T, ylim = c(0, 0.2),main="",ylab="Freque^ncia",xlab="Classes",col="lightgreen")

> curve(dchisq(x, df = 4), col = 2, lty = 2, lwd = 2, add = TRUE)

> BMI hist(BMI,prob=T, xlab="Classes", main="", col="lightgreen", xlim=c(15,35), ylim=c(0, .20),ylab="Frequencia")

> curve(dnorm(x, mean=mean(BMI), sd=sd(BMI)), add=TRUE, lty=2, col = 2)

Classes

Freq

unc

ia

0 2 4 6 8 10 12 14

0.00

0.05

0.10

0.15

0.20

Classes

Freq

uenc

ia

15 20 25 30 35

0.00

0.05

0.10

0.15

0.20

Figura 5: Curvas

6

> par(mfrow=c(2,3))

> BMI histinfo hist(BMI, breaks=20, main="Breaks=20")

> hist(BMI, breaks=5, main="Breaks=5")

> hist(BMI,breaks = "Sturges", main="Sturges")

> hist(BMI, freq=FALSE, main="Density plot")

> hist(BMI, freq=FALSE, xlab="Body Mass Index", main="Distribution of Body Mass Index", col="lightgreen", xlim=c(15,35), ylim=c(0,0.20))

> curve(dnorm(x, mean=mean(BMI), sd=sd(BMI)), add=TRUE, col="darkblue", lwd=2)

Histogram of BMI

BMI

Freq

uenc

y

20 25 30

050

100

150

Breaks=20

BMI

Freq

uenc

y

20 25 30

050

100

150

Breaks=5

BMI

Freq

uenc

y

15 20 25 30 35

010

020

030

040

050

060

0

Sturges

BMI

Freq

uenc

y

20 25 30

050

100

150

Density plot

BMI

Den

sity

20 25 30

0.00

0.05

0.10

0.15

Distribution of Body Mass Index

Body Mass Index

Den

sity

15 20 25 30 35

0.00

0.05

0.10

0.15

0.20

Figura 6: curvas

7

> x=rnorm(999900,70,6)

> y=rnorm(999900,12,2)

> z=rnorm(999900,30,1)

> hist(x,xlim=range(x,y,z),main="")

> hist(y,add=TRUE, border='blue')> hist(z,add=TRUE, border='red')

x

Freq

uenc

y

0 20 40 60 80 100

050

000

1000

0015

0000

2000

0025

0000

3000

00

Figura 7: Adicionar novo histograma num ja existente

8

> par(mfrow=c(1,2))

> curve(dnorm(x,mean=1,sd=sqrt(2)),lwd=2,from=-6,to=17,ylab="Densidade")

> curve(dnorm(x,mean=10,sd=sqrt(2)),col=2,lwd=2,add=T)

> legend('topright',legend=c(expression(mu==1),expression(mu==10)),text.col=c(1,2),cex=1.5,bty="n")> title("Comparando a distribuic~ao normal com medias diferentes")

> curve(dnorm(x,mean=4,sd=sqrt(3)),lwd=2,from=-15,to=25,ylab="Densidade")

> curve(dnorm(x,mean=4,sd=sqrt(20)),col=2,lwd=2,add=T)

> legend('topright',legend=c(expression(sigma^2==3),expression(sigma^2==20)),text.col=c(1,2),cex=1.5,bty="n")> title("Comparando a distribuic~ao normal com varia^ncias diferentes")

5 0 5 10 15

0.00

0.05

0.10

0.15

0.20

0.25

x

Den

sida

de

= 1 = 10

Comparando a distribuio normal com mdias diferentes

10 0 10 20

0.00

0.05

0.10

0.15

0.20

x

Den

sida

de

2 = 32 = 20

Comparando a distribuio normal com varincias diferentes

Figura 8: Curva normal

9

> require(agricolae)

> data(growth)

> attach(growth)

> h1 normal.freq(h1,col="blue")

Classes

Freq

unc

ia

6 8 10 12 14 16

02

46

8

Figura 9: Curva ao histograma

10

> TCF=seq(68,80)

> RG=seq(32,37)

> f1=function(x,y)3696.046+62.183*x-33.414*y

> PCS=outer(TCF,RG,f1)

> persp(TCF,RG,PCS,theta=-50,phi=-0,shade=2,ticktype = "simple",box=T,border=T,r = sqrt(20), d = 1,col="blue")

TCFRG

PCS

Figura 10: 3D

11

> require(lattice)

> b0

[11,] "c" "5.09184556095365"

[12,] "c" "5.00340549838538"

[13,] "c" "4.73680938607341"

> colnames(test) Trat Trat

[1] "a" "b" "a" "b" "a" "b" "a" "b" "a" "b" "c" "c" "c"

> Resposta Resposta

[1] "6.73639601703966" "4.08931250781133" "4.6219843334702" "5.10199927530805"

[5] "4.28784778823526" "5.3458819717746" "3.55664508814877" "5.90286217355144"

[9] "6.71059636804312" "6.37526460075227" "5.09184556095365" "5.00340549838538"

[13] "4.73680938607341"

> dad=subset(test,Trat %in% c("a","c")) # Selecionar apenas trat a e c

> dad

trat resposta

[1,] "a" "6.73639601703966"

[2,] "a" "4.6219843334702"

[3,] "a" "4.28784778823526"

[4,] "a" "3.55664508814877"

[5,] "a" "6.71059636804312"

[6,] "c" "5.09184556095365"

[7,] "c" "5.00340549838538"

[8,] "c" "4.73680938607341"

> dad2=subset(test,Trat=="a") # Selecionar apenas trat a

> dad2

trat resposta

[1,] "a" "6.73639601703966"

[2,] "a" "4.6219843334702"

[3,] "a" "4.28784778823526"

[4,] "a" "3.55664508814877"

[5,] "a" "6.71059636804312"

4 Editor LaTeX + R + Sweave no linux

Ao usar R ambiente Linux para sair do manual R basta teclar Q e para terminar a execucao de determinadoprocesso Ctrl+C alem de quando precisar do caminho de algum arquivo no linux basta um Ctrl+L e acompilacao de arquivos no Linux e feita por R CMD Sweave meu arquivo.rnw e para extrair apenaso codigo R R CMD Stangle meu arquivo.rnw . Dois plugins especiais para o gedit com finalidade deintegracao entre com R.sao:

1. Gedit-latex-plugin que pode ser instalado com;sudo aptitude purge gedit-latex-plugin via terminal do Linux.

2. Rgedit que integra R com Gedit.

http://www.stattler.com/article/using-gedit-or-rgedit-r

http://sourceforge.net/projects/rgedit/

A melhor altenativa de IDE/livre e o Rstudio interface multiplataforma que tem perfeita integracao com R eem 4 janelas pre definidas, perimite atuliazar pacotes e instalar pacotes de forma interativa. Pode ser baixadaem http://www.rstudio.org/

1. R+LATEX+ SWEAVE

2. R+LATEX+Knitr

3. R+HTML+Markdown

13

4.1 Sweave linux

Para compilar o arquivo .rnw basta usar dentro do R:

1. setwd(/home/andre/pasta)

2. Sweave(nomedoarquivo.rnw)

3. system(pdflatex nomedoarquivo.tex)

No RStudio a compilacao e feita direta no menu compilacaocompilePdf(/media/KINGSTON/apostilaR2/apostilaR.rnw)

4.2 knitr no linux

Geracao de relatorios dinamicos elegante, flexvel e rapida com R. O pacote foi projetado para ser um me-canismo de transparencia para a geracao de relatorios dinamicos com R, resolver alguns problemas de longadata em Sweave. knitr tem mesmo proposito do sweave. Instale o pacote knitr e crie seu arquivo .rnw ins-tall.packages(knitr, dep=T) nao e necesario carregar nenhum pacote no .rnw nem mesmo a bibliotecaknitr precisa ser carregada. A complilacao no RStudio e feita direta no menu compilacao desde que modifiqueo tipo de arquivo .rnw na configuracao do RStudio.

4.3 Salvar imagens linux

Salvar uma imagem deve-se proceder da seguinte forma.

> pdf("arquivo.pdf",wid=640, hei=480) # Salvar em pdf

> plot(1:10, 1:10)

> dev.off()

null device

1

> png("arquivo.png",wid=640, hei=480) # Salvar em png

> plot(1:10, 1:10)

> dev.off()

null device

1

4.4 Sites importantes de apoio ao R na internet

http://ridiculas.wordpress.com

http://www.leg.ufpr.br/~paulojus/embrapa/Rembrapa

https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br

www.leg.ufpr.br/ce223

http://cran.r-project.org/web/views/

http://www.leg.ufpr.br/doku.php/pessoais:walmes

4.5 A mudanca de diretorio

A mudanca de diretorio pode ser feita da seguinte forma.

> getwd() # Area de trabalho atual

[1] "/media/Meus arquivos/Pendrive/KINGSTON/sweave/apostilaR2"

> setwd("/media/Meus arquivos/Pendrive/KINGSTON/sweave/apostilaR2") # Mudanca de diretorio

14

4.6 Configuracoes basicas para Linux

Para ler diretamente um banco de dados direto no formato .xls uma alternativa no Linux e o pacote gdata. Quedeve ser instalado inicialmente com o comando direto no terminal R install.packages(gdata, dep=T) ouR CMD INSTALL nome do pacote .tar.gz para instalar um pacote ja baixado assim como a atualizacao deum pacote especfico pode ser feita por update.packages(ask=gdata) e update.packages(ask= ) ouupdate.packages(ask=F) para atualizar as bibliotecas ja instaladas. Para ver os banco de dados instaladose disponveis no R data(). A edicao de um banco de dados pode ser feita com o comando fix(dados) ouedit(dados) sendo que o ideal e configurar o editor de dados primeiramente options(editor=gedit).Instalar um arquivo tar.gz basta install.packages(caminho/PACOTE.tar.gz,repos=NULL) a versao decada pacote pode ser feita require(pacote) e depois veja as informacoes dessa sessao com sessionInfo(). Parasaber se esta instalado e onde esta find.package(NomeDoPacote) a descricao completa do pacote package-Description(NomeDoPacote) para ver apenas a versao packageDescription(NomeDoPacote)$Version.Para carregar um pacote de oculto (silencioso) require(xtable, quietly = TRUE) e instalar mais de umpacote ao mesmo tempo install.packages(c(knitr, xtable), dependencies = TRUE). Os pacotescarregados em uma secao do R podem ser vistos com search()

4.7 Ajuda e demostracoes

Para obter ajuda em HTML no Linux basta usar a rotina.

> require(xtable)

> help(print.xtable, help_type="html")

Para obter a demonstracao pode se usar a funcaodemo

> demo(plotmath)

4.8 Interface amigavel para R

O pacote Rcmdr foi desenvolvido por John Fox visando atender a esta demanda. Para utilizar este pacote bastainstala-lo e carregar com o comando require(Rcmdr) e o menu se abrira automaticamente. Atencao: Note queo Rcmdr nao prove acesso a toda funcionalidade do R mas simplesmente a alguns procedimentos estatsticosmais usuais.

> install.packages("Rcmdr",dep=T)

> requre(Rcmdr)

4.9 Adicionando repositorios

1. Ver Codename no terminal Linux lsb_release -a

2. sudo gedit /etc/apt/sources.list

3. deb http://cran-r.c3sl.ufpr.br/bin/linux/ubuntu Codename/

4. sudo apt-get update&&sudo apt-get upgrade

5. gpg keyserver subkeys.pgp.net recv-key E2A11821

6. gpg -a export E2A11821 | sudo apt-key add -

Ou coloca esta chave abaixo se a de acima nao for aceita.

7. gpg keyserver keyserver.ubuntu.com recv-key E084DAB9

8. gpg -a export E084DAB9 | sudo apt-key add -9. Pode adicionar chave via arquivo sudo apt-key add arquivo.txt

A verificacao do Codename tambem pode ser consultada aqui http://cran-r.c3sl.ufpr.br/bin/linux

4.10 Gerador de tabela em LaTeX e HTML

http://www.tablesgenerator.com/

15

4.11 Manipulacao de diretorios

> getwd()

[1] "/media/Meus arquivos/Pendrive/KINGSTON/sweave/apostilaR2"

> setwd("/media/Meus arquivos/Pendrive/KINGSTON/sweave/apostilaR2")

> require(xtable)

> require(gdata)

> dad str(dad) # Imprimindo caractersticas de banco de dados dad

'data.frame': 198 obs. of 6 variables:$ fruta : Factor w/ 3 levels "Abricc","Batata",..: 3 3 3 3 3 3 3 3 3 3 ...

$ situa : Factor w/ 2 levels "Agitac","estac": 2 2 2 2 2 2 2 2 2 2 ...

$ cond : Factor w/ 2 levels "ci","si": 2 2 2 2 2 2 2 2 2 2 ...

$ tempo : Factor w/ 11 levels "10D","12D","14D",..: 8 9 10 11 1 2 3 4 5 6 ...

$ qtde : Factor w/ 3 levels "e0","e100","e200": 1 1 1 1 1 1 1 1 1 1 ...

$ medida: num 1.99 0.9 0.89 1.92 1.93 1.93 1.93 1.93 2.04 2.04 ...

> attach(dad) # Reconhecendo as variaveis do banco de dados dad

> detach(dad) # Deixando de reconhecer as variaveis do banco de dados dad

> summary(dad) # Estatsticas descritivas do banco de dados dad

fruta situa cond tempo qtde medida

Abricc :66 Agitac:99 ci:132 10D :18 e0 :66 Min. :0.530

Batata :66 estac :99 si: 66 12D :18 e100:66 1st Qu.:1.440

Cupuacu:66 14D :18 e200:66 Median :1.950

16D :18 Mean :1.838

18D :18 3rd Qu.:2.098

20D :18 Max. :3.160

(Other):90

> head(dad) # seis primeiras linhas do banco de dados


1 Cupuacu estac si 2D e0 1.99






> dad2 head(dad2)








> dados4 dados4

x y z

1 2 3 3

2 3 6 4

3 4 7 5

4 5 8 6

5 6 9 8

6 5 0 9

> library(XLConnect)

> teste1 teste1

x y z

1 2 3 3

2 3 6 4

16

3 4 7 5

4 5 8 6

5 6 9 8

6 5 0 9

> teste2 head(teste2)








> require(xlsx)


> dados=read.xlsx("teste.xlsx",sheetIndex=1)

> str(dados)

'data.frame': 6 obs. of 3 variables:$ x: num 2 3 4 5 6 5

$ y: num 3 6 7 8 9 0

$ z: num 3 4 5 6 8 9

> summary(dados)

x y z

Min. :2.000 Min. :0.00 Min. :3.000

1st Qu.:3.250 1st Qu.:3.75 1st Qu.:4.250

Median :4.500 Median :6.50 Median :5.500

Mean :4.167 Mean :5.50 Mean :5.833

3rd Qu.:5.000 3rd Qu.:7.75 3rd Qu.:7.500

Max. :6.000 Max. :9.00 Max. :9.000

> dados3 dados3

Tra.Bloco.Y

1 a,A,83

2 a,B,63

3 a,C,55

4 b,A,86

5 b,B,69

6 b,C,61

7 c,A,103

8 c,B,79

9 c,C,79

10 d,A,116

11 d,B,81

12 d,C,79

13 e,A,132

14 e,B,98

15 e,C,91

> setwd("/media/Meus arquivos/Pendrive/KINGSTON/sweave/apostilaR2") # Nova area de trabalho

Analise descritiva de algumas variaveis do banco de dados dad.

> attach(dad)

> summary(medida) # Resumo descritivo da variavel medida

Min. 1st Qu. Median Mean 3rd Qu. Max.

0.530 1.440 1.950 1.838 2.098 3.160

> dados

4.12 Box plot com media destacadas

> par(mfrow=c(1,2))

> hist(medida,xlab="Medidas",ylab="Densidade",main="",col="limegreen")

> boxplot(medida~fruta,col="blue",main="")

> points(1:nlevels(fruta), tapply(medida,fruta, mean),pch=19,col="red")

Medidas

Den

sida

de

0.5 1.0 1.5 2.0 2.5 3.0 3.5

010

2030

4050

60

Abric Batata Cupuau

0.5

1.0

1.5

2.0

2.5

3.0

Figura 12: Estatsticas descritivas

5 Funcao densidade probabilidade fdp

Plotando uma fdp definida por mais de uma sentenca.

> x y=0 & x=1 & x=2 &x plot(x,y,type="l",ylab="f(x)")

1 0 1 2 3 4 5

0.0

0.1

0.2

0.3

0.4

0.5

x

f(x)

Figura 13: Funcao densidade probabilidade fdp

> f1

+ ret

+ }

> integrate(f1,-1,5) # Integrac~ao da func~ao f(x) entre -1 e 5

1 with absolute error < 7.4e-05

Para arredondar o valor da integral integrate(f1,-1,5)$value

6 Matematica

Determinacao dos divisores de um numero natural

> require(divisors)

> divisors(131)

$num

[1] 2

$divs

[1] 1 131

A resolucao de uma equacao pode ser feita assim. Considere f(x) = 3 4.x + x2

> polyroot(c(3, -4, 1))

[1] 1+0i 3+0i

O Resto da divisao de a por b e dado por

> 150 %% 13

[1] 7

C49,6 =49!

6!(496)!

> factorial(49)/(factorial(6)*factorial(49-6))

[1] 13983816

> y=c(1,2,3,4,4,5,6,7,8,0,9)

> y

[1] 1 2 3 4 4 5 6 7 8 0 9

> as.integer(y >= 3) # Maiores iguais a 3 coloque 1 menores que tre^s coloque 0

[1] 0 0 1 1 1 1 1 1 1 0 1

7 Tabela da ANOVA

As funcoes lm, aov e Anova (pacote car) sao as mais comuns para imprimir a tabela da ANOVA. Apenas parafins de ilustracao sera construido 3 modelos, cada um com uma das funcoes para fins de comparacao.


> require(gdata)

> dad2 attach(dad2)

> require(car)

> fit1=lm(Y~Tra+Bloco)

> options(scipen=5)

> summary(fit1)

Call:

lm(formula = Y ~ Tra + Bloco)

Residuals:

Min 1Q Median 3Q Max

-5.0 -3.0 -1.0 3.5 6.0

19

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 86.000 3.276 26.250 0.00000000477 ***

Trab 5.000 3.916 1.277 0.237456

Trac 20.000 3.916 5.108 0.000921 ***

Trad 25.000 3.916 6.384 0.000213 ***

Trae 40.000 3.916 10.215 0.00000723969 ***

BlocoB -26.000 3.033 -8.572 0.00002647667 ***

BlocoC -31.000 3.033 -10.220 0.00000721153 ***

---

Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

Residual standard error: 4.796 on 8 degrees of freedom

Multiple R-squared: 0.9696, Adjusted R-squared: 0.9467

F-statistic: 42.46 on 6 and 8 DF, p-value: 0.00001226

> options(show.signif.stars=FALSE) # Retira a indicac~ao de fatores significativos

> require(xtable)

> a=anova(fit1)

> x=xtable(a,caption ="Tabela ANOVA",align="llllll")

> print(x,caption.placement='top',table.placement="H")

Tabela 1: Tabela ANOVADf Sum Sq Mean Sq F value Pr(>F)

Tra 4 3090.00 772.50 33.59 0.0000Bloco 2 2770.00 1385.00 60.22 0.0000Residuals 8 184.00 23.00

> u=Anova(fit1)

> xtable(u)

Sum Sq Df F value Pr(>F)Tra 3090.00 4 33.59 0.0000Bloco 2770.00 2 60.22 0.0000Residuals 184.00 8

> p=Anova(fit1,type="III")

> xtable(p)

Sum Sq Df F value Pr(>F)(Intercept) 15848.57 1 689.07 0.0000Tra 3090.00 4 33.59 0.0000Bloco 2770.00 2 60.22 0.0000Residuals 184.00 8

20

7.1 Distribuicao Binomial

> par(mfrow=c(2,2))

> p n x px barplot(px,names=as.character(x))

> plot(0:8,choose(8,0:8),type="h",main="Coeficiente binomial")

> barplot(dbinom(0:4,4,0.1),names=as.character(0:4),xlab="x",ylab="f(x)")

> plot(dbinom(0:16, 16, 0.3), type = "h", xlab = "Numeros de fe^meas",ylab = "Probabilidade")

>

0 1 2 3 4 5 6 7 8 9 10

0.00

0.15

0 2 4 6 80

2050

Coeficiente binomial

0:8

choo

se(8,

0:8)

0 1 2 3 4

x

f(x)

0.0

0.3

0.6

5 10 15

0.00

0.10

0.20

Numeros de fmeas

Prob

abilid

ade

Figura 14: Distribuicao de probabilidade binomial

21

7.2 BoxPlot

> par(mfrow=c(1,2))

> x=rnorm(33)

> boxplot(x)

> boxplot(x,notch=T)

2

1

01

2

1

01

Figura 15: Boxplot

8 Multivariada

Uma funcao que pode ajudar nas analises e manova

> tear gloss opacity Y rate additive fit summary.aov(fit) # Tabela de anova univariadas

Response tear :

Df Sum Sq Mean Sq F value Pr(>F)

rate 1 1.7405 1.74050 15.7868 0.001092 **

additive 1 0.7605 0.76050 6.8980 0.018330 *

rate:additive 1 0.0005 0.00050 0.0045 0.947143

Residuals 16 1.7640 0.11025

---

Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

Response gloss :


rate 1 1.3005 1.30050 7.9178 0.01248 *

additive 1 0.6125 0.61250 3.7291 0.07139 .

rate:additive 1 0.5445 0.54450 3.3151 0.08740 .

Residuals 16 2.6280 0.16425

---

22

Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

Response opacity :


rate 1 0.421 0.4205 0.1036 0.7517

additive 1 4.901 4.9005 1.2077 0.2881

rate:additive 1 3.960 3.9605 0.9760 0.3379

Residuals 16 64.924 4.0578

> summary(fit) # Tabela manova

Df Pillai approx F num Df den Df Pr(>F)

rate 1 0.61814 7.5543 3 14 0.003034 **

additive 1 0.47697 4.2556 3 14 0.024745 *

rate:additive 1 0.22289 1.3385 3 14 0.301782

Residuals 16

---

Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

A extracao de elementos necessarios podem ser feito por fit$Eigenvalues ou fit$SS

9 Modelos lineares generalizados - GLM

Para pesquisar sobre alguma funcao pode ser feito por RSiteSearch(glm) ou ?glm

> counts outcome treatment dados dados

treatment outcome counts

1 1 1 18

2 1 2 17

3 1 3 15

4 2 1 20

5 2 2 10

6 2 3 20

7 3 1 25

8 3 2 13

9 3 3 12

> str(dados)

'data.frame': 9 obs. of 3 variables:$ treatment: Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 3 3 3

$ outcome : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3

$ counts : num 18 17 15 20 10 20 25 13 12

> fit1 anova(fit1,test="Chisq")

Analysis of Deviance Table

Model: poisson, link: log

Response: counts

Terms added sequentially (first to last)

Df Deviance Resid. Df Resid. Dev Pr(>Chi)

NULL 8 10.5814

outcome 2 5.4523 6 5.1291 0.06547 .

treatment 2 0.0000 4 5.1291 1.00000

---

Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

23

> summary(fit1)

Call:

glm(formula = counts ~ outcome + treatment, family = poisson())

Deviance Residuals:

1 2 3 4 5 6 7 8

-0.67125 0.96272 -0.16965 -0.21999 -0.95552 1.04939 0.84715 -0.09167

9

-0.96656

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) 3.045e+00 1.709e-01 17.815 summary(step(fit1,direction="backward"),trace=0)

Start: AIC=56.76

counts ~ outcome + treatment

Df Deviance AIC

- treatment 2 5.1291 52.761

5.1291 56.761

- outcome 2 10.5814 58.214

Step: AIC=52.76

counts ~ outcome

Df Deviance AIC

5.1291 52.761

- outcome 2 10.5814 54.214

Call:

glm(formula = counts ~ outcome, family = poisson())

Deviance Residuals:


-0.9666 -0.6713 -0.1696 0.8471 1.0494

Coefficients:


(Intercept) 3.0445 0.1260 24.165

Null deviance: 10.5814 on 8 degrees of freedom

Residual deviance: 5.1291 on 6 degrees of freedom

AIC: 52.761

Number of Fisher Scoring iterations: 4

> summary(step(fit1,direction="forward"),trace=0)

Start: AIC=56.76

counts ~ outcome + treatment

Call:

glm(formula = counts ~ outcome + treatment, family = poisson())

Deviance Residuals:

1 2 3 4 5 6 7 8

-0.67125 0.96272 -0.16965 -0.21999 -0.95552 1.04939 0.84715 -0.09167

9

-0.96656

Coefficients:


(Intercept) 3.045e+00 1.709e-01 17.815 require(logistf) # Penalizac~ao Firth 1993

> require(brglm) # Penalizac~ao Firth 1993

9.4 Teste estatsticos

Uma tabela de dupla entrada pode ser inserida no R assim.

> matriz colnames(matriz)

> rownames(matriz) matriz

Disease Control

Exposed 156 9421

Unexposed 1531 14797

Teste em tabelas de contingencia Teste de quiquadrado e Fisher.

> x x

B

A 0-6 meses 6-12 meses 12-18 meses

autistas 101 106 77

deficientes 58 32 55

tpicos 113 165 98

> chisq.test(x)

Pearson's Chi-squared test

data: x

X-squared = 22.1991, df = 4, p-value = 0.0001829

> fisher.test(x,workspace=10000000)

Fisher's Exact Test for Count Data

data: x

p-value = 0.0001183

alternative hypothesis: two.sided

Uma opcao muito legal e chisq.test(x)$expected para calcular a esperanca das caselas. Ou usar o pacoteepicalc para montar uma tabela.

> require(epicalc)

> dados dados

Admit Gender Dept Freq

1 Admitted Male A 512

2 Rejected Male A 313

3 Admitted Female A 89

4 Rejected Female A 19

5 Admitted Male B 353

6 Rejected Male B 207

7 Admitted Female B 17

8 Rejected Female B 8

9 Admitted Male C 120

10 Rejected Male C 205

11 Admitted Female C 202

12 Rejected Female C 391

13 Admitted Male D 138

14 Rejected Male D 279

15 Admitted Female D 131

16 Rejected Female D 244

17 Admitted Male E 53

18 Rejected Male E 138

19 Admitted Female E 94

20 Rejected Female E 299

21 Admitted Male F 22

22 Rejected Male F 351

23 Admitted Female F 24

24 Rejected Female F 317

> xtabs(Freq ~ Gender + Admit, dados)

26

Admit

Gender Admitted Rejected

Male 1198 1493

Female 557 1278

9.5 Estatstica Descritivas

> counts outcome treatment dados dados

treatment outcome counts

1 1 1 18

2 1 2 17

3 1 3 15

4 2 1 20

5 2 2 10

6 2 3 20

7 3 1 25

8 3 2 13

9 3 3 12

> str(dados)

'data.frame': 9 obs. of 3 variables:$ treatment: Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 3 3 3

$ outcome : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3

$ counts : num 18 17 15 20 10 20 25 13 12

> tapply(counts,outcome, mean)

1 2 3

21.00000 13.33333 15.66667

> tapply(counts,outcome, sd)

1 2 3

3.605551 3.511885 4.041452

> tapply(counts,outcome, var)

1 2 3

13.00000 12.33333 16.33333

> by(counts,outcome, mean)

outcome: 1

[1] 21

------------------------------------------------------------

outcome: 2

[1] 13.33333

------------------------------------------------------------

outcome: 3

[1] 15.66667

Ou obter as tres estatsticas em uma mesma tabela.

> library(doBy)

> summaryBy(counts~outcome,dados,FUN=c(mean,sd,var))

outcome counts.mean counts.sd counts.var

1 1 21.00000 3.605551 13.00000

2 2 13.33333 3.511885 12.33333

3 3 15.66667 4.041452 16.33333

> require(fdth)

> x fdt(x)

> d=fdt(x)

> d

27

Class limits f rf rf(%) cf cf(%)

[2.99,3.6) 6 0.06 6 6 6

[3.6,4.21) 15 0.15 15 21 21

[4.21,4.81) 18 0.18 18 39 39

[4.81,5.42) 31 0.31 31 70 70

[5.42,6.03) 15 0.15 15 85 85

[6.03,6.64) 11 0.11 11 96 96

[6.64,7.24) 3 0.03 3 99 99

[7.24,7.85) 1 0.01 1 100 100


> require(gdata)

> dad str(dad)

'data.frame': 36 obs. of 5 variables:$ sexo : Factor w/ 2 levels "feminino ","masculino ": 2 2 1 2 1 2 1 2 1 2 ...

$ esc : Factor w/ 3 levels "1 grau","2 grau",..: 1 1 1 2 1 1 1 1 2 2 ...

$ sal : num 4 4.56 5.25 5.73 6.26 6.66 6.86 7.39 7.59 7.84 ...

$ idad : int 26 32 36 20 40 28 41 43 34 23 ...

$ filhos: int 0 1 2 0 1 1 0 0 1 0 ...

> attach(dad)

> require(xtable)

> s print(s,table.placement = "H", caption.placement='top',floating=TRUE)

28

Tabela 2: Tabela de frequencia empresasexo esc sal idad filhos

1 masculino 1 grau 4.00 26 02 masculino 1 grau 4.56 32 13 feminino 1 grau 5.25 36 24 masculino 2 grau 5.73 20 05 feminino 1 grau 6.26 40 16 masculino 1 grau 6.66 28 17 feminino 1 grau 6.86 41 08 masculino 1 grau 7.39 43 09 feminino 2 grau 7.59 34 110 masculino 2 grau 7.84 23 011 feminino 2 grau 8.12 33 212 masculino 1 grau 8.46 27 013 feminino 2 grau 8.74 37 014 masculino 1 grau 8.95 44 315 feminino 2 grau 9.13 30 016 masculino 2 grau 9.35 38 317 feminino 2 grau 9.77 31 118 masculino 1 grau 9.80 39 219 feminino superior 10.53 25 120 masculino 2 grau 10.76 37 021 feminino 2 grau 11.06 30 122 masculino 1 grau 11.59 34 223 feminino 2 grau 12.00 41 424 masculino superior 12.79 26 025 feminino 2 grau 13.23 32 226 masculino 2 grau 13.60 35 227 feminino 1 grau 13.85 46 028 masculino 2 grau 14.69 29 229 masculino 2 grau 14.71 40 530 masculino 2 grau 15.99 35 231 feminino superior 16.22 31 332 masculino 2 grau 16.61 36 133 feminino superior 17.26 43 334 masculino superior 18.75 33 435 feminino 2 grau 19.40 48 236 masculino superior 23.30 42 3

> fit=lm(sal~esc+sexo)

> xtable(anova(fit))

Df Sum Sq Mean Sq F value Pr(>F)esc 2 307.82 153.91 11.85 0.0001sexo 1 10.47 10.47 0.81 0.3759Residuals 32 415.49 12.98

> xtable(summary(fit))

Estimate Std. Error t value Pr(>|t|)(Intercept) 7.0698 1.3221 5.35 0.0000

esc2 grau 3.9540 1.3583 2.91 0.0065escsuperior 8.8557 1.8132 4.88 0.0000

sexomasculino 1.0990 1.2240 0.90 0.3759

> require(fdth)

> tab1=fdt(idad, start=20,h=6,end=50)

> tab1

29


[20,26) 3 0.08 8.33 3 8.33

[26,32) 9 0.25 25.00 12 33.33

[32,38) 12 0.33 33.33 24 66.67

[38,44) 9 0.25 25.00 33 91.67

[44,50) 3 0.08 8.33 36 100.00

> tab2=fdt(idad,breaks="Sturges")

> tab2


[19.8,23.9) 2 0.06 5.56 2 5.56

[23.9,28) 4 0.11 11.11 6 16.67

[28,32.1) 8 0.22 22.22 14 38.89

[32.1,36.2) 8 0.22 22.22 22 61.11

[36.2,40.3) 6 0.17 16.67 28 77.78

[40.3,44.4) 6 0.17 16.67 34 94.44

[44.4,48.5) 2 0.06 5.56 36 100.00

> dados=cbind(rnorm(8000),rnorm(8000))

> tab3=fdt(dados,breaks="Sturges")

> tab3

Column:1


[-3.5,-2.95) 13 0.00 0.16 13 0.16

[-2.95,-2.4) 56 0.01 0.70 69 0.86

[-2.4,-1.85) 174 0.02 2.17 243 3.04

[-1.85,-1.3) 507 0.06 6.34 750 9.38

[-1.3,-0.753) 1094 0.14 13.68 1844 23.05

[-0.753,-0.205) 1545 0.19 19.31 3389 42.36

[-0.205,0.344) 1800 0.22 22.50 5189 64.86

[0.344,0.892) 1388 0.17 17.35 6577 82.21

[0.892,1.44) 875 0.11 10.94 7452 93.15

[1.44,1.99) 381 0.05 4.76 7833 97.91

[1.99,2.54) 122 0.02 1.52 7955 99.44

[2.54,3.09) 37 0.00 0.46 7992 99.90

[3.09,3.64) 5 0.00 0.06 7997 99.96

[3.64,4.18) 3 0.00 0.04 8000 100.00

Column:2


[-4.02,-3.47) 2 0.00 0.02 2 0.02

[-3.47,-2.92) 15 0.00 0.19 17 0.21

[-2.92,-2.37) 66 0.01 0.82 83 1.04

[-2.37,-1.82) 209 0.03 2.61 292 3.65

[-1.82,-1.27) 538 0.07 6.73 830 10.38

[-1.27,-0.723) 1043 0.13 13.04 1873 23.41

[-0.723,-0.174) 1528 0.19 19.10 3401 42.51

[-0.174,0.375) 1732 0.22 21.65 5133 64.16

[0.375,0.924) 1411 0.18 17.64 6544 81.80

[0.924,1.47) 884 0.11 11.05 7428 92.85

[1.47,2.02) 390 0.05 4.88 7818 97.73

[2.02,2.57) 147 0.02 1.84 7965 99.56

[2.57,3.12) 30 0.00 0.38 7995 99.94

[3.12,3.67) 5 0.00 0.06 8000 100.00

> plot(tab1, ylab="Frequencia",xlab="")

30

Freq

uenc

ia

02

46

810

1214

20 26 32 38 44 50

> plot(tab3, ylab="Frequencia",xlab="")

Column:1

Freq

uenc

ia

010

00

3.50 2.40 1.30 0.20 0.89 1.99 3.09 4.18

Column:2

Freq

uenc

ia

010

0020

00

4.02 2.92 1.82 0.72 0.38 1.47 2.57 3.67

> median(sal)

[1] 10.165

> dad2=subset(dad, sexo == "masculino ")

> dad2

sexo esc sal idad filhos

1 masculino 1 grau 4.00 26 0












24 masculino superior 12.79 26 0





31




> summary(dad2)

sexo esc sal idad filhos

feminino : 0 1 grau :8 Min. : 4.000 Min. :20.00 Min. :0.00

masculino :20 2 grau :9 1st Qu.: 7.728 1st Qu.:27.75 1st Qu.:0.00

superior:3 Median :10.280 Median :34.50 Median :1.50

Mean :11.277 Mean :33.35 Mean :1.55

3rd Qu.:14.695 3rd Qu.:38.25 3rd Qu.:2.25

Max. :23.300 Max. :44.00 Max. :5.00

> library(doBy)

> summaryBy(sal~sexo,dad,FUN=c(mean))

sexo sal.mean

1 feminino 10.95438

2 masculino 11.27650

> summaryBy(sal~esc,dad,FUN=c(mean,sd))

esc sal.mean sal.sd

1 1 grau 7.80250 2.905913

2 2 grau 11.57333 3.691888

3 superior 16.47500 4.502438

> fit=aov(sal~sexo)

> popMeans(fit, effect="sexo")

beta0 Estimate Std.Error t.value DF Pr(>|t|) Lower Upper sexo

1 0 10.95437 1.160672 9.437963 34 0 8.595607 13.31314 feminino

2 0 11.27650 1.038136 10.862255 34 0 9.166753 13.38625 masculino

> TukeyHSD(fit)

Tukey multiple comparisons of means

95% family-wise confidence level

Fit: aov(formula = sal ~ sexo)

$sexo

diff lwr upr p adj

masculino -feminino 0.322125 -2.842495 3.486745 0.8373522

32

> require(epicalc)

> xtabs(sal~sexo, dad) # Soma salario por sexo

sexo

feminino masculino

175.27 225.53

> xtabs(sal~esc, dad) # Soma salario por escolaridade

esc

1 grau 2 grau superior

93.63 208.32 98.85

> addmargins(xtabs(sal~esc, dad)) # Adiciona total marginal na tabela

esc

1 grau 2 grau superior Sum

93.63 208.32 98.85 400.80

> summ(sal,box=T,main="Distribuic~ao da variavel resposta",ylab="")

obs. mean median s.d. min. max.

36 11.133 10.165 4.58 4 23.3

> table(sexo)

sexo

feminino masculino

16 20

> prop.table(table(sexo))

sexo

feminino masculino

0.4444444 0.5555556

> addmargins(table(sexo),FUN = mean)

sexo

feminino masculino mean

16 20 18

5 10 15 20

Distribuio da varivel resposta

Figura 16: dist

33

> xtabs(~sexo+idad,dad) # fatores cruzados

idad

sexo 20 23 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44

feminino 0 0 1 0 0 0 0 2 2 1 1 1 0 1 1 0 0 1 2 0 1 0

masculino 1 1 0 2 1 1 1 0 0 1 1 1 2 1 1 1 1 1 0 1 1 1

idad

sexo 46 48

feminino 1 1

masculino 0 0

> aggregate(sal~sexo,data=dad,mean) # media salario por sexo.

sexo sal

1 feminino 10.95438


> aggregate(sal~esc, data=dad,sd) # media salario por ecolaridade.

esc sal

1 1 grau 2.905913

2 2 grau 3.691888

3 superior 4.502438

> table(sexo)

sexo

feminino masculino

16 20

> table(esc)

esc

1 grau 2 grau superior

12 18 6

> require(gmodels)

> CrossTable(esc,sexo, expected = TRUE)

Conteudo das celulas

|-------------------------|

| N |

| N esperado |

| Contribuic~ao para Qui2 |

| N / Total da linha |

| N / Total da coluna |

| N / Total da tabela |

|-------------------------|

==========================================

sexo

esc feminino masculino Total

------------------------------------------

1 grau 4 8 12

5.3 6.7

0.333 0.267

0.333 0.667 0.333

0.250 0.400

0.111 0.222

------------------------------------------

2 grau 9 9 18

8.0 10.0

0.125 0.100

0.500 0.500 0.500

0.562 0.450

0.250 0.250

------------------------------------------

superior 3 3 6

2.7 3.3

0.042 0.033

34

0.500 0.500 0.167

0.188 0.150

0.083 0.083

------------------------------------------

Total 16 20 36

0.444 0.556

==========================================

> require(descr)

> CrossTable(sexo)


|-------------------------|

| N |


|-------------------------|

| feminino | masculino |

|------------|------------|

| 16 | 20 |

| 0.444 | 0.556 |

|------------|------------|

> CrossTable(sal)


|-------------------------|

| N |


|-------------------------|

| 4 | 4.56 | 5.25 | 5.73 | 6.26 | 6.66 | 6.86 | 7.39 | 7.59 |

|-------|-------|-------|-------|-------|-------|-------|-------|-------|

| 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |

| 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 |

|-------|-------|-------|-------|-------|-------|-------|-------|-------|

| 7.84 | 8.12 | 8.46 | 8.74 | 8.95 | 9.13 | 9.35 | 9.77 | 9.8 |

|-------|-------|-------|-------|-------|-------|-------|-------|-------|

| 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |

| 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 |

|-------|-------|-------|-------|-------|-------|-------|-------|-------|

| 10.53 | 10.76 | 11.06 | 11.59 | 12 | 12.79 | 13.23 | 13.6 | 13.85 |

|-------|-------|-------|-------|-------|-------|-------|-------|-------|

| 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |

| 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 |

|-------|-------|-------|-------|-------|-------|-------|-------|-------|

| 14.69 | 14.71 | 15.99 | 16.22 | 16.61 | 17.26 | 18.75 | 19.4 | 23.3 |

|-------|-------|-------|-------|-------|-------|-------|-------|-------|

| 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |

| 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 |

|-------|-------|-------|-------|-------|-------|-------|-------|-------|

> CrossTable(esc)


|-------------------------|

| N |


|-------------------------|

| 1 grau | 2 grau | superior |

|----------|----------|----------|

| 12 | 18 | 6 |

35

| 0.333 | 0.500 | 0.167 |

|----------|----------|----------|

> CrossTable(filhos)


|-------------------------|

| N |


|-------------------------|

| 0 | 1 | 2 | 3 | 4 | 5 |

|--------|--------|--------|--------|--------|--------|

| 11 | 8 | 9 | 5 | 2 | 1 |

| 0.306 | 0.222 | 0.250 | 0.139 | 0.056 | 0.028 |

|--------|--------|--------|--------|--------|--------|

> require(descr)

> descr(sexo)

feminino masculino

16 20

> descr(sal)

Min. 1st Qu. Median Mean 3rd Qu. Max.

4.000 7.778 10.160 11.130 14.060 23.300

> compmeans(sal,sexo)

Valor medio de "sal" segundo "sexo"

Media N Desv. Pd.

feminino 10.95438 16 4.114950

masculino 11.27650 20 5.020285

Total 11.13333 36 4.578760

> crosstab(sexo,esc) # Otimo para cruzar duas variaveis categoricas


|-------------------------|

| Contagem |

|-------------------------|

================================================

esc

sexo 1 grau 2 grau superior Total

------------------------------------------------

feminino 4 9 3 16

------------------------------------------------

masculino 8 9 3 20

------------------------------------------------

Total 12 18 6 36

================================================

> table(sexo,esc) # Otimo para cruzar duas variaveis categoricas

esc

sexo 1 grau 2 grau superior

feminino 4 9 3

masculino 8 9 3

> table(sexo,esc,filhos) # Otimo para cruzar duas variaveis categoricas

, , filhos = 0

esc


feminino 2 2 0

masculino 3 3 1

, , filhos = 1

esc

36


feminino 1 3 1

masculino 2 1 0

, , filhos = 2

esc


feminino 1 3 0

masculino 2 3 0

, , filhos = 3

esc


feminino 0 0 2

masculino 1 1 1

, , filhos = 4

esc


feminino 0 1 0

masculino 0 0 1

, , filhos = 5

esc


feminino 0 0 0

masculino 0 1 0

> require(gregmisc)

> CrossTable(esc, format = "SPSS")


|-------------------------|

| Contagem |

| Percentual por linha |

|-------------------------|


|----------|----------|----------|

| 12 | 18 | 6 |

| 33.333 | 50.000 | 16.667 |

|----------|----------|----------|

> CrossTable(esc, format = "SAS")


|-------------------------|

| N |


|-------------------------|


|----------|----------|----------|

| 12 | 18 | 6 |

| 0.333 | 0.500 | 0.167 |

|----------|----------|----------|

> # Estatsticas descritivas

> library(fBasics)

> basicStats(sal)

37

sal

nobs 36.000000

NAs 0.000000

Minimum 4.000000

Maximum 23.300000

1. Quartile 7.777500

3. Quartile 14.060000

Mean 11.133333

Median 10.165000

Sum 400.800000

SE Mean 0.763127

LCL Mean 9.584104

UCL Mean 12.682563

Variance 20.965046

Stdev 4.578760

Skewness 0.600356

Kurtosis -0.318921

> require(AdequacyModel)

> hist(sal)

> descriptive(sal)

$mean

[1] 11.13333

$median

[1] 10.165

$mode

[1] 9

$variance

[1] 20.96505

$Skewness

[1] 0.62627

$Kurtosis

[1] -0.16353

$minimum

[1] 4

$maximum

[1] 23.3

$n

[1] 36

> ftable(sexo)

sexo feminino masculino

16 20

> require(psych)

> describe(sal)

var n mean sd median trimmed mad min max range skew kurtosis se

1 1 36 11.13 4.58 10.16 10.85 4.72 4 23.3 19.3 0.6 -0.32 0.76

> require(Hmisc)

> describe(sal)

> require(epicalc)

> shapiro.qqnorm(sal)

38

2 1 0 1 2

510

1520

Normal QQ plot of sal

Theoretical Quantiles

Sam

ple

Quan

tiles

ShapiroWilk test P value = 0.2545

> require(fBasics)

> qqnormPlot(sal)

2 1 0 1 2

1

01

2

Normal Quantiles

SS.1

Ord

ered

Dat

a

NORM QQ PLOT

Conf

iden

ce In

terv

als

: 95%

9.6 Biblioteca ExpDes

A biblioteca para analise numa unica rodada pode ser encontrada no site do autor https://sites.google.com/site/ericbferreira/home neste caso baixei a versao ExpDes.pt para linux e compilei via terminal.

R CMD INSTALL ExpDes.pt_1.1.2.tar.gz

ou instale via install.packages(ExpDes.pt, dep=T) nos repositorios estao mantidas as duas versoes em portu-gues e ingles.

> require(ExpDes.pt) # vers~ao portugue^s

> dic(esc, sal, quali = TRUE, mcomp = "tukey", sigT = 0.05, sigF = 0.05)

------------------------------------------------------------------------

Quadro da analise de variancia

------------------------------------------------------------------------

GL SQ QM Fc Pr>Fc

Tratamento 2 307.82 153.909 11.924 0.00012669

Residuo 33 425.96 12.908

Total 35 733.78

------------------------------------------------------------------------

CV = 32.27 %

39

------------------------------------------------------------------------

Teste de normalidade dos residuos (Shapiro-Wilk)

p-valor: 0.4820254

De acordo com o teste de Shapiro-Wilk a 5% de significancia, os residuos podem ser considerados normais.

------------------------------------------------------------------------

Teste de Tukey

------------------------------------------------------------------------

Grupos Tratamentos Medias

a superior 16.475

b 2 grau 11.57333

b 1 grau 7.8025

------------------------------------------------------------------------

> dic(sexo, sal, quali = TRUE, mcomp = "tukey", sigT = 0.05, sigF = 0.05)

------------------------------------------------------------------------

Quadro da analise de variancia

------------------------------------------------------------------------

GL SQ QM Fc Pr>Fc

Tratamento 1 0.92 0.9224 0.042792 0.83735

Residuo 34 732.85 21.5545

Total 35 733.78

------------------------------------------------------------------------

CV = 41.7 %

------------------------------------------------------------------------

Teste de normalidade dos residuos (Shapiro-Wilk)

p-valor: 0.2993022

De acordo com o teste de Shapiro-Wilk a 5% de significancia, os residuos podem ser considerados normais.

------------------------------------------------------------------------

De acordo com o teste F, as medias nao podem ser consideradas diferentes.

Niveis Medias

1 feminino 10.95438


------------------------------------------------------------------------

> require(ExpDes) # Vers~ao ingle^s

> crd(sexo, sal, quali = TRUE, mcomp= "tukey", sigF = 0.05)

------------------------------------------------------------------------

Analysis of Variance Table

------------------------------------------------------------------------

DF SS MS Fc Pr>Fc

Treatament 1 0.92 0.9224 0.042792 0.83735

Residuals 34 732.85 21.5545

Total 35 733.78

------------------------------------------------------------------------

CV = 41.7 %

------------------------------------------------------------------------

Shapiro-Wilk normality test

p-value: 0.2993022

According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.

------------------------------------------------------------------------

According to the F test, the means can not be considered distinct.

Levels Means

1 feminino 10.95438


------------------------------------------------------------------------

40

As funcoes disponveis no pacote ExpDes (nomes em portugues) podem ser aplicadas aos seguintes tipos dedelineamentos experimentais, sempre considerando fatores de efeito fixo:

funcao descricaodic() Delineamento Inteiramente Casualizado balanceado com um so fator.dbc() Delineamento em Blocos Casualizados balanceado com um so fator.dql() Delineamento em Quadrado Latino balanceado com um so fator.fat2.dic() Delineamento Inteiramente Casualizado balanceado em fatorial duplo.fat2.dbc() Delineamento em Blocos Casualizados balanceado em fatorial duplo.fat2.ad.dic() Delineamento Inteiramente Casualizado balanceado em fatorial duplo com um tratamento adicio-nal.fat2.ad.dbc() Delineamento em Blocos Casualizados balanceado em fatorial duplo com um tratamento adicional.fat3.dic() Delineamento Inteiramente Casualizado balanceado em fatorial triplo.fat3.dbc() Delineamento em Blocos Casualizados balanceado em fatorial triplo.fat3.ad.dic() Delineamento Inteiramente Casualizado balanceado em fatorial triplo com um tratamento adicio-nal.fat3.ad.dbc() Delineamento em Blocos Casualizados balanceado em fatorial triplo com um tratamento adicional.psub2.dic() Delineamento Inteiramente Casualizado balanceado em esquema de parcelas subdivididas.psub2.dbc() Delineamento em Blocos Casualizados balanceado em esquema de parcelas subdivididas.

10 Ler funcao externa


> source("fun.R") # Ler func~ao que construi

> f

function (n)

{

y f(10)

> f(100)

> f(1000)

> f(2000)

41

Time

y

2 4 6 8 10

0.0

0.2

0.4

0.6

0.8

1.0

Time

y

0 20 40 60 80 100

0.0

0.2

0.4

0.6

0.8

1.0

Time

y

0 200 400 600 800 1000

0.0

0.2

0.4

0.6

0.8

1.0

Time

y

0 500 1000 1500 2000

0.0

0.2

0.4

0.6

10.1 Dados Multinacional


> require(gdata)

> dados str(dados)

'data.frame': 50 obs. of 9 variables:$ gerente : int 1 2 3 4 5 6 7 8 9 10 ...

$ sexo : Factor w/ 2 levels " Feminino"," Masculino": 2 1 2 2 1 1 2 2 2 1 ...

$ estadocivil : Factor w/ 6 levels " casado"," casado",..: 1 5 2 2 2 2 4 3 3 2 ...

$ idade : int 40 25 28 25 33 42 45 55 44 43 ...

$ filhos : int 2 0 3 2 1 4 2 5 2 2 ...

$ salario : num 15.3 20.4 8.5 16.8 21.7 22.5 29 27.3 19.3 23.4 ...

$ escolaridade: Factor w/ 3 levels " medio"," pos-graduac~ao",..: 3 2 1 3 3 2 2 2 3 3 ...

$ ingle^s : Factor w/ 3 levels " n~ao"," sim",..: 2 3 1 3 3 3 3 3 3 3 ...

$ origem : Factor w/ 4 levels " capital SP",..: 1 2 4 1 2 2 3 4 3 1 ...

> attach(dados)

The following object is masked from dad (position 4):

filhos, sexo

> table(sexo)

sexo

Feminino Masculino

18 32

> pie(table(sexo))

> pie(table(escolaridade))

> pie(table(estadocivil))

> fit=aov(salario~escolaridade+sexo)

> summary(fit)


escolaridade 2 1749 874.5 25.78 0.000000031 ***

sexo 1 56 56.0 1.65 0.205

42

Residuals 46 1561 33.9

---

Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

> library(doBy)

> summaryBy(salario~sexo,dados,FUN=c(mean,sd))

sexo salario.mean salario.sd

1 Feminino 24.13889 8.248789

2 Masculino 24.61250 8.436356

> summaryBy(salario~escolaridade,dados,FUN=c(mean,sd))

escolaridade salario.mean salario.sd

1 medio 11.45000 3.370905

2 pos-graduac~ao 32.60000 7.354998

3 superior 24.08485 5.642535

> library(fBasics)

> basicStats(salario)

salario

nobs 50.000000

NAs 0.000000

Minimum 6.500000

Maximum 41.000000

1. Quartile 19.300000

3. Quartile 29.750000

Mean 24.442000

Median 23.300000

Sum 1222.100000

SE Mean 1.172063

LCL Mean 22.086652

UCL Mean 26.797348

Variance 68.686567

Stdev 8.287736

Skewness 0.205725

Kurtosis -0.492068

> require(descr)

> crosstab(sexo,escolaridade) # otimo cruzar duas variaveis


|-------------------------|

| Contagem |

|-------------------------|

=============================================================

escolaridade

sexo medio pos-graduac~ao superior Total

-------------------------------------------------------------

Feminino 2 6 10 18

-------------------------------------------------------------

Masculino 4 5 23 32

-------------------------------------------------------------

Total 6 11 33 50

=============================================================

> require(gmodels)

> CrossTable(sexo,escolaridade, expected = TRUE)


|-------------------------|

| N |

| N esperado |

| Contribuic~ao para Qui2 |


| N / Total da coluna |

| N / Total da tabela |

43

|-------------------------|

=============================================================

escolaridade

sexo medio pos-graduac~ao superior Total

-------------------------------------------------------------

Feminino 2 6 10 18

2.2 4.0 11.9

0.012 1.051 0.298

0.111 0.333 0.556 0.360

0.333 0.545 0.303

0.040 0.120 0.200

-------------------------------------------------------------

Masculino 4 5 23 32

3.8 7.0 21.1

0.007 0.591 0.167

0.125 0.156 0.719 0.640

0.667 0.455 0.697

0.080 0.100 0.460

-------------------------------------------------------------

Total 6 11 33 50

0.120 0.220 0.660

=============================================================

> boxplot(salario~sexo)

> boxplot(salario~escolaridade)

44


> require(gdata)

> dados str(dados)




$ idade : int 40 25 28 25 33 42 45 55 44 43 ...

$ filhos : int 2 0 3 2 1 4 2 5 2 2 ...

$ salario : num 15.3 20.4 8.5 16.8 21.7 22.5 29 27.3 19.3 23.4 ...




> require(gplots)

> plotmeans(salario~escolaridade,p=0.95,col="red",connect =F,n.label=F,ci.label=F,mean.labels=F, xlab="Escolaridade",ylab="Media")

1015

2025

3035

Escolaridade

Md

ia

mdio psgraduao superior

Figura 17: Intervalo de confianca para media I

45


> require(gdata)

> dados extratificac~ao medias desviopadr~ao n erro require(gplots)

> plotCI(x=medias, uiw=erro)

1.0 1.5 2.0 2.5 3.0

1015

2025

3035

medi

as

Figura 18: Intervalo de confianca para media II

46

> da da$y m0 new new$pred str(new)

'data.frame': 5 obs. of 2 variables:$ trat: Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5

$ pred: num [1:5, 1:3] 0.736 2.176 2.724 4.391 4.903 ...

..- attr(*, "dimnames")=List of 2

.. ..$ : chr "1" "2" "3" "4" ...

.. ..$ : chr "fit" "lwr" "upr"

> ylim bp arrows(bp, new$pred[,2], bp, new$pred[,3], code=3, angle=90)

1 2 3 4 5

01

23

45

Figura 19: Intervalo de confianca para media III

10.2 Grafico Lattice

A latice tem opcoes de histogram(salario) densityplot(salario) algumas outras funcoes.

barchart bar chart x~A or A~x

bwplot boxplot x~A or A~x

cloud 3D scatterplot z~x*y|A

contourplot 3D contour plot z~x*y

densityplot kernal density plot ~x|A*B

dotplot dotplot ~x|A

histogram histogram ~x

levelplot 3D level plot z~y*x

parallel parallel coordinates plot data frame

splom scatterplot matrix data frame

stripplot strip plots A~x or x~A

xyplot scatterplot y~x|A

47

wireframe 3D wireframe graph z~y*x


> require(gdata)

> dados str(dados)




$ idade : int 40 25 28 25 33 42 45 55 44 43 ...

$ filhos : int 2 0 3 2 1 4 2 5 2 2 ...

$ salario : num 15.3 20.4 8.5 16.8 21.7 22.5 29 27.3 19.3 23.4 ...




> require(lattice)

> xyplot(salario~sexo,jitter.x=TRUE)

sexo

sala

rio

10

20

30

40

Feminino Masculino

Figura 20: Graficos Lattice I

48


> require(gdata)

> dados str(dados)




$ idade : int 40 25 28 25 33 42 45 55 44 43 ...

$ filhos : int 2 0 3 2 1 4 2 5 2 2 ...

$ salario : num 15.3 20.4 8.5 16.8 21.7 22.5 29 27.3 19.3 23.4 ...




> require(lattice)

> xyplot(salario~escolaridade,jitter.x=TRUE)

escolaridade

sala

rio

10

20

30

40


Figura 21: Graficos Lattice II


> require(gdata)

> dados str(dados)




$ idade : int 40 25 28 25 33 42 45 55 44 43 ...

$ filhos : int 2 0 3 2 1 4 2 5 2 2 ...

$ salario : num 15.3 20.4 8.5 16.8 21.7 22.5 29 27.3 19.3 23.4 ...




49

> library(doBy)

> require(lattice)

> summaryBy(salario~escolaridade|sexo,dados,FUN=c(mean))

escolaridade sexo salario.mean

1 medio Feminino 12.80000

2 medio Masculino 10.77500

3 pos-graduac~ao Feminino 32.03333

4 pos-graduac~ao Masculino 33.28000

5 superior Feminino 21.67000

6 superior Masculino 25.13478

> aggregate(salario~sexo+escolaridade, data=dados,mean)

sexo escolaridade salario

1 Feminino medio 12.80000

2 Masculino medio 10.77500

3 Feminino pos-graduac~ao 32.03333

4 Masculino pos-graduac~ao 33.28000

5 Feminino superior 21.67000

6 Masculino superior 25.13478

> summaryBy(salario~escolaridade,dados,FUN=c(mean))

escolaridade salario.mean

1 medio 11.45000

2 pos-graduac~ao 32.60000

3 superior 24.08485

> fit=aov(salario~escolaridade)

> popMeans(fit, effect="escolaridade")

beta0 Estimate Std.Error t.value DF Pr(>|t|) Lower Upper

1 0 11.45000 2.394288 4.782215 47 0.0000175 6.633314 16.26669

2 0 32.60000 1.768299 18.435798 47 0.0000000 29.042642 36.15736

3 0 24.08485 1.020928 23.591137 47 0.0000000 22.031007 26.13869

escolaridade

1 medio

2 pos-graduac~ao

3 superior

> fit2=aov(salario~sexo+escolaridade)

> popMeans(fit2, c("sexo", "escolaridade"))

beta0 Estimate Std.Error t.value DF Pr(>|t|) Lower Upper

1 0 9.948168 2.649797 3.754314 46 0.0004869 4.614404 15.28193

2 0 12.200916 2.448708 4.982593 46 0.0000093 7.271922 17.12991

3 0 31.576024 1.928657 16.372030 46 0.0000000 27.693840 35.45821

4 0 33.828771 1.999833 16.915802 46 0.0000000 29.803317 37.85423

5 0 22.514752 1.588136 14.176839 46 0.0000000 19.318000 25.71150

6 0 24.767499 1.144777 21.635212 46 0.0000000 22.463182 27.07182

sexo escolaridade

1 Feminino medio

2 Masculino medio

3 Feminino pos-graduac~ao

4 Masculino pos-graduac~ao

5 Feminino superior

6 Masculino superior

> fit3=lmBy(salario~escolaridade|sexo, data=dados)

> summary(fit3)

Call:

lm(formula = salario ~ escolaridade, data = wd)

Residuals:


-9.835 -4.020 -1.235 4.905 12.865

Coefficients:

50


(Intercept) 10.775 3.021 3.567 0.001278 **

escolaridade pos-graduac~ao 22.505 4.053 5.553 0.00000546 ***

escolaridade superior 14.360 3.273 4.387 0.000139 ***

---

Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1




Call:

lm(formula = salario ~ escolaridade, data = wd)

Residuals:


-11.6333 -2.2950 -0.0017 3.2300 8.5667

Coefficients:


(Intercept) 12.800 3.936 3.252 0.005358 **

escolaridade pos-graduac~ao 19.233 4.545 4.232 0.000724 ***

escolaridade superior 8.870 4.311 2.057 0.057472 .

---

Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1




> require(lsmeans)

> fit4 lsmeans(fit4, pairwise~escolaridade)

$èscolaridade lsmeansèscolaridade lsmean SE df lower.CL upper.CL

medio 11.45000 2.394288 47 6.633314 16.26669

pos-graduac~ao 32.60000 1.768299 47 29.042642 36.15736

superior 24.08485 1.020928 47 22.031007 26.13869

$èscolaridade pairwise differencesèstimate SE df t.ratio p.value

medio - pos-graduac~ao -21.150000 2.976491 47 -7.10568 0.00000

medio - superior -12.634848 2.602865 47 -4.85421 0.00004

pos-graduac~ao - superior 8.515152 2.041856 47 4.17030 0.00038

p values are adjusted using the tukey method for 3 means

51

> barchart(salario~sexo|escolaridade,col="limegreen",layout = c(3, 1))

sala

rio

10

20

30

40

Feminino Masculino

mdio

Feminino Masculino

psgraduao

Feminino Masculino

superior

Figura 22: Graficos Lattice

52

> densityplot(~salario|sexo, plot.points = FALSE,ylab="Densidade")

salario

Den

sida

de

0.00

0.02

0.04

0.06

0 10 20 30 40 50

Feminino

0 10 20 30 40 50

Masculino

Figura 23: Graficos Lattice VI

53

> histogram(~salario|sexo)

salario

Perc

ent o

f Tota

l

0

10

20

30

40

10 20 30 40

Feminino

10 20 30 40

Masculino


54

> histogram(~salario|escolaridade,layout = c(3, 1))

salario

Perc

ent o

f Tota

l

0

10

20

30

40

50

10 20 30 40

mdio

10 20 30 40

psgraduao

10 20 30 40

superior


55

> cloud(salario~filhos*escolaridade)

filhosescolaridade

salario


56


> require(gdata)

> dados str(dados)




$ idade : int 40 25 28 25 33 42 45 55 44 43 ...

$ filhos : int 2 0 3 2 1 4 2 5 2 2 ...

$ salario : num 15.3 20.4 8.5 16.8 21.7 22.5 29 27.3 19.3 23.4 ...




> require(lattice)

> xyplot(salario~sexo|escolaridade,layout = c(3, 1),jitter.x=TRUE)

sexo

sala

rio

10

20

30

40

Feminino Masculino

mdio

Feminino Masculino

psgraduao

Feminino Masculino

superior

Figura 27: Graficos Lattice III

57


> require(gdata)

> dados str(dados)




$ idade : int 40 25 28 25 33 42 45 55 44 43 ...

$ filhos : int 2 0 3 2 1 4 2 5 2 2 ...

$ salario : num 15.3 20.4 8.5 16.8 21.7 22.5 29 27.3 19.3 23.4 ...




> require(lattice)

> bwplot(salario~escolaridade|sexo)

sala

rio

10

20

30

40


Feminino


Masculino

Figura 28: Graficos Lattice IV

58

> qqmath(~salario|escolaridade,layout = c(3, 1))

qnorm

sala

rio

10

20

30

40

2 1 0 1 2

mdio

2 1 0 1 2

psgraduao

2 1 0 1 2

superior

Figura 29: Normalidade por grupo

59


> require(gdata)

> dados str(dados)




$ idade : int 40 25 28 25 33 42 45 55 44 43 ...

$ filhos : int 2 0 3 2 1 4 2 5 2 2 ...

$ salario : num 15.3 20.4 8.5 16.8 21.7 22.5 29 27.3 19.3 23.4 ...




> attach(dados)

The following objects are masked from dados (position 5):

escolaridade, estadocivil, filhos, gerente, idade, ingle^s, origem,

salario, sexo


filhos, sexo

> require(lattice)

> barchart(salario~sexo,col="limegreen")

sala

rio

10

20

30

40

Feminino Masculino

Figura 30: Graficos Lattice V

10.3 ScottKnott, o TukeyC

O pacote facilita bastante a aplicacao do teste de Tukey para os delineamentos simples (DIC, DBC, DQL) eesquemas experimentais simples e usuais (Fatorias com no maximo 3 fatores, parcela sub-dividida e sub-sub-dividida). Para os que ja conhecem o pacote ScottKnott, o TukeyC tem a mesma estrutura basica, porem aplicao teste de Tukey.

60

> require(laercio)

> result=aov(salario~escolaridade)

> LTukey(result,"escolaridade",conf.level=0.95)

TUKEY TEST TO COMPARE MEANS

Confidence level: 0.95

Dependent variable: salario

Variation Coefficient: 23.9947 %

Independent variable: escolaridade

Factors Means

pos-graduac~ao 32.6 a

superior 24.0848484848485 b

medio 11.45 c

> LDuncan(result,"escolaridade",conf.level = 0.95)

DUNCAN TEST TO COMPARE MEANS

Confidence Level: 0.95

Dependent Variable: salario

Variation Coefficient: 23.9947 %

Independent Variable: escolaridade

Factors Means

pos-graduac~ao 32.6 a

superior 24.0848484848485 b

medio 11.45 c

> LScottKnott(result,"escolaridade",conf.level = 0.95)

SCOTT-KNOTT ORIGINAL TEST

Confidence Level: 0.95

Independent variable: escolaridade

FACTORS MEANS

1 pos-graduac~ao 32.60000 A

2 superior 24.08485 A

3 medio 11.45000 A


> require(gdata)

> dad str(dad)

'data.frame': 12 obs. of 3 variables:$ trata: Factor w/ 4 levels "A","B","C","D": 1 1 1 2 2 2 3 3 3 4 ...

$ rep : int 1 2 3 1 2 3 1 2 3 1 ...

$ respo: int 6 14 4 9 10 8 73 22 10 1 ...

> attach(dad)

> dad

trata rep respo

1 A 1 6

2 A 2 14

3 A 3 4

4 B 1 9

5 B 2 10

6 B 3 8

7 C 1 73

8 C 2 22

9 C 3 10

10 D 1 1

61

11 D 2 10

12 D 3 4

> library(doBy)

> summaryBy(respo~trata,dad,FUN=c(mean,sd))

trata respo.mean respo.sd

1 A 8 5.291503

2 B 9 1.000000

3 C 35 33.451457

4 D 5 4.582576

> require(ScottKnott)

> sk1 summary(sk1)

> plot(sk1, title = "Treatamentos", col = rainbow(3))

62


> require(gdata)

> dados str(dados)




$ idade : int 40 25 28 25 33 42 45 55 44 43 ...

$ filhos : int 2 0 3 2 1 4 2 5 2 2 ...

$ salario : num 15.3 20.4 8.5 16.8 21.7 22.5 29 27.3 19.3 23.4 ...




> attach(dados)



salario, sexo



salario, sexo


filhos, sexo

> library(doBy)

> require(lattice)

> barchart(table(sexo,escolaridade))

> summaryBy(salario~sexo|escolaridade,dados,FUN=c(mean))

sexo escolaridade salario.mean







> aggregate(salario~sexo+escolaridade, data=dados,mean)

sexo escolaridade salario







Freq

Feminino

Masculino

0 10 20 30

Figura 31: Plotar cruzamento de variaveis

63

10.4 TukeyC


> library(TukeyC)

> require(gdata)

> dad_b str(dad_b)

> av summary(av)

> tk_b summary(tk_b)

> plot(tk_b)

64

10.5 UsingR


> require(gdata)

> dados str(dados)




$ idade : int 40 25 28 25 33 42 45 55 44 43 ...

$ filhos : int 2 0 3 2 1 4 2 5 2 2 ...

$ salario : num 15.3 20.4 8.5 16.8 21.7 22.5 29 27.3 19.3 23.4 ...




> attach(dados)



salario, sexo



salario, sexo



salario, sexo


filhos, sexo

> require(UsingR)

> par(mfrow=c(2,2))

> boxplot(salario~sexo,col = c("palevioletred1","royalblue2"))

> simple.violinplot(salario~sexo,col="brown")

> simple.densityplot(salario~sexo)

> boxplot(salario~sexo,,notch=T)

Feminino Masculino

1015

2025

3035

40

Feminino Masculino

010

2030

4050

0 10 20 30 40 50

0.00

0.01

0.02

0.03

0.04

0.05

0.06

0.07

Feminino Masculino

Feminino Masculino

1015

2025

3035

40

Figura 32: Funcao UsingR

66

> par(mfrow=c(1,2))

> require(UsingR)

> x=rnorm(200)

> y=x+rnorm(200)

> simple.scatterplot(x,y,col=3)

> simple.hist.and.boxplot(x, main="")

1

2

3

Figura 33: Histograma marginal

67

> par(fig=c(0,0.8,0,0.8), new=TRUE)

> x=rnorm(330)

> y=x+rnorm(330)

> plot(x, y, xlab="Teste",ylab="Teste")

> par(fig=c(0,0.8,0.55,1), new=TRUE)

> boxplot(x, horizontal=TRUE, axes=FALSE)

> par(fig=c(0.65,1,0,0.8),new=TRUE)

> boxplot(y,axes=FALSE)

> mtext("Teste", side=3, outer=TRUE, line=-3)

2 1 0 1 2

4

2

02

4

Teste

Test

e

Teste

Figura 34: Box plot marginal

68

> x=rnorm(3000)

> layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))

> hist(x,main="",ylab="Densidade",prob=T)



x

Den

sida

de

4 2 0 2 4

0.0

0.1

0.2

0.3

0.4

x

Den

sida

de

4 2 0 2 4

0.0

0.1

0.2

0.3

0.4

x

Den

sida

de

4 2 0 2 4

0.0

0.1

0.2

0.3

0.4

Figura 35: Tres graficos em em mesma janela

69

> x=rnorm(3000)

> layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE),widths=c(3,1), heights=c(1,2))




x

Den

sida

de

3 2 1 0 1 2 3

0.0

0.1

0.2

0.3

0.4

x

Den

sida

de

3 2 1 0 1 2 3

0.0

0.1

0.2

0.3

0.4

x

Den

sida

de

3 2 1 0 1 2 3

0.0

0.1

0.2

0.3

0.4

Figura 36: Tres graficos em em mesma janela

70


> require(gdata)

> dados str(dados)




$ idade : int 40 25 28 25 33 42 45 55 44 43 ...

$ filhos : int 2 0 3 2 1 4 2 5 2 2 ...

$ salario : num 15.3 20.4 8.5 16.8 21.7 22.5 29 27.3 19.3 23.4 ...




> attach(dados)



salario, sexo



salario, sexo



salario, sexo



salario, sexo


filhos, sexo

> require(UsingR)

> par(mfrow=c(1,3))

> boxplot(salario,col = "palevioletred1",horizontal=T,xlab="Salario")

> hist(salario,main="",ylab="Freque^ncia",prob=T,col="palevioletred1",

+ xlab="Salario",xlim=c(0,60))

> curve(dnorm(x,mean=mean(salario),sd=sd(salario)),col=4,lty=2,lwd=2,add=TRUE)

> qqnorm(salario,col="palevioletred1")

10 15 20 25 30 35 40

Salrio Salrio

Freq

unc

ia

0 10 20 30 40 50 60

0.00

0.01

0.02

0.03

0.04

0.05

2 1 0 1 2

1015

2025

3035

40

Normal QQ Plot

Theoretical Quantiles

Sam

ple

Quan

tiles

Figura 37: Histograma Box plot71

> par(mfrow=c(2,2))

> means barplot(means,xlab="Escolaridade",ylab="Media por escolaridade",col=c("red","blue","pink"))

> barplot(tapply(salario,list(sexo,escolaridade),mean),beside=T,ylim=c(0,30),

+ col=c("forestgreen", "palegreen"),ylab="Media por escolaridade sexo")

> legend("topleft", ncol=2, legend=c("Feminino", "Masculino"),fill=c("forestgreen", "palegreen"),bty="n",xpd=TRUE)

> barplot(table(sexo,escolaridade),beside=TRUE,col=c("forestgreen", "palegreen"))


> barplot(table(sexo,escolaridade),col=c("forestgreen", "palegreen"))



Escolaridade

Md

ia p

or e

scol

arid

ade

05

1015

2025

30

mdio psgraduao superiorM

dia

por

esc

olar

idad

e se

xo

05

1015

2025

30

Feminino Masculino


05

1015

20

Feminino Masculino


05

1015

2025

30

Feminino Masculino

Figura 38: Media por fator

> barplot(table(sexo,escolaridade),col=c("forestgreen", "palegreen"),legend=T)


Masculino Feminino

05

1015

2025

30

Figura 39: Media por fator

72

> par(mfrow=c(1,3))

> interaction.plot(sexo,escolaridade, salario)

> interaction.plot(escolaridade,sexo, salario)

> interaction.plot(escolaridade,sexo, salario,fun = mean, ylab ="Medias")

1015

2025

30

sexo

mean o

f sa

lario

Feminino Masculino

escolaridade

psgraduao superior mdio

1015

2025

30

escolaridade

mean o

f sa

lario


sexo

Masculino Feminino

1015

2025

30escolaridade

Md

ias


sexo

Masculino Feminino

Figura 40: Interacao

> par(mfrow=c(1,3))

> plot.design(salario~sexo)

> plot.design(salario~escolaridade)

> plot.design(salario~escolaridade+sexo)

24.2

24.3

24.4

24.5

24.6

Factors

mean o

f sal

ario

Feminino

Masculino

sexo

1520

2530

Factors

mean o

f sal

ario

mdio

psgraduao

superior

escolaridade

1520

2530

Factors

mean o

f sal

ario

mdio

psgraduao

superior Feminino

Masculino

escolaridade sexo

Figura 41: Design

73

> coplot(salario~escolaridade|sexo,panel=panel.smooth)

md p sprr

1015

2025

3035

40

md p sprr

escolaridade

sala

rio

Feminino

Masculino

Given : sexo

Figura 42: Coplot

74

10.6 Regressao linear


> require(gdata)

> dados str(dados)

'data.frame': 10 obs. of 2 variables:$ x: num 5 5.4 5.7 5.9 6.3 6.8 7.2 7.3 7.6 7.8

$ y: num 10 10.9 11.4 11.5 12 11.6 11.1 10.5 10.1 9.6

> attach(dados)

The following object is masked _by_ .GlobalEnv:

x, y

> fit=lm(y~poly(x,2,raw = TRUE),data=dados) # Ajuste func~ao segundo grau (modelo quadratico)

> summary(fit)

Call:

lm(formula = y ~ poly(x, 2, raw = TRUE), data = dados)

Residuals:

Documents

apostila_R2_André Oliveira Souza