## Este script contem o exercicio sobre pressão sanguinea 
 # e também sobre analise de variancia e lm com diagnosticos
  # aqui serao utilizados os bancos pressao sanguinea, cafe, milhas, e altura2022

## pacotes para baixar ou elicitar 
library(psych)
library (arsenal)
library (knitr)
## Exercício 1 - Bioestatística 


# Questões 1 a Questão 3
  # importar dados, 



setwd ("E:/2021/Curso-R/Pressao-sanguinea") # set working directory = setwd

          ## outros bancos que serão utilizados 
          pressaotxt <- read.table (file = 'PRESSAO_SANGUINEA.txt', header = TRUE)
          cafe <- read.table (file = 'cafe-anova-oneway.txt', header = TRUE)
          alt <- read.table (file = 'altura2022.txt', header = TRUE )
          
          
                        ######### RECODIFICACAO DE VARIAVEIS #### 
                        #cada vez que eu recodificar algo para um exercicio vou colocar aqui
                        # para facilitar a vida e até exportar o banco se quiser
                        pressaotxt$pesokg <- pressaotxt$PESO*0.454
                        
                        pressaotxt$alturam <- pressaotxt$ALTURA*0.305
                        
                        pressaotxt$sexocat <- as.factor(pressaotxt$SEXO)
                        
                        #########################################
                      
                        # se quiser salvar definitivamente
                        write.csv (pressaotxt,"E:/2021/Curso-R/Pressao-sanguinea/pressao18-11-22", row.names = FALSE )
                        write.table(pressaotxt,"E:/2021/Curso-R/Pressao-sanguinea/pressao18-11-22txt", row.names = FALSE )
                        write.table(pressaotxt,"E:/2021/Curso-R/Pressao-sanguinea/pressao18-11-22txt", na='.')
                          # este ultimo formato de exportação pode ser util porque exporta os NA como dot ou o que quiser. 
                          # isso é util para o SAs que nao reconhece NA e sim dot. 
                        


View(pressaotxt) # pegadinha V tem letra maiuscula
dim (pressaotxt)
nrow (pressaotxt)
ncol (pressaotxt)

# Questão 4 - monte o livro de codigos
  #preciso saber os tipos de variaveis e seus nomes e seus valores

str(pressaotxt)


# Questao 6 - classifique as variáveis se quantitativas e tipo:
str (pressaotxt) # 

# Questão 7  - transformação das variaveis peso e altura
  # transformação de variaveis peso e altura

 pressaotxt$pesokg <- pressaotxt$PESO*0.454

 pressaotxt$alturam <- pressaotxt$ALTURA*0.305
 
 
# Questao 8 e 9 - Analise descritiva de todas as variaveis (medidas de tendência central + dispersão) e questao 9 graficos
   # preciso ver se há necessidade de limpar as variaveis e decidir qual melhor medida 
    # especificamente para cada variavel
 
# questao 10 e 11 -descrever as distribuicoes das variaveis continuas vamos fazer junto + kurtose e skweness - cv
  ## CV = desvio padrao divido pela media, o desvio é grade em relacao a media? 
 
 # elicitar arsenal e knitr que eu gosto..
 
       library (arsenal)
       library (knitr)
 
 
 # IDADE
 
       summary (freqlist(~ IDADE , data = pressaotxt ))
       
       hist(pressaotxt$IDADE)
       boxplot (pressaotxt$IDADE)
       qqnorm (pressaotxt$IDADE, pch = 1)
       qqline (pressaotxt$IDADE)
       shapiro.test(pressaotxt$IDADE)
       # distribuição não é normal pelo jeitao mas pelo teste da normal borderline (logo mediana e q1, q3)
       
       summary (pressaotxt$IDADE) # se quiser apenas reportar mediana e q1,q3
       describe(pressaotxt$IDADE) # em termos de normalidade nao ta das piores
       
       ## Coeficiente de variação nao tem no basico do R assim podemos resolver escrevendo uma função
       cv <- function(variable) { sd(variable) / mean(variable) }
       cv (pressaotxt$IDADE)     
  
  # PESOKG
       
       summary (freqlist(~ pressaotxt$pesokg ))
       
       hist(pressaotxt$pesokg)
       boxplot (pressaotxt$pesokg)
       qqnorm (pressaotxt$pesokg, pch = 1)
       qqline (pressaotxt$pesokg)
       
       str(pressaotxt)
       
       ## testes para normalidade
        (pressaotxt$pesokg)
       shapiro.test(pressaotxt$pesokg) # para amostras pequenas
       
       
       summary (pressaotxt$pesokg)
       describe (pressaotxt$pesokg)       
       cv(pressaotxt$pesokg)
       
       
  # Alturam
       
       summary (freqlist(~ pressaotxt$alturam ))
       
       hist(pressaotxt$alturam)
       boxplot (pressaotxt$alturam)
       # distribuição não é normal  
       
       summary (pressaotxt$alturam)
       describe (pressaotxt$alturam)       
       shapiro.test(pressaotxt$alturam)
       
       
   # Sistolica
       
       summary (freqlist(~ pressaotxt$SISTOLICA))
       
       hist(pressaotxt$SISTOLICA)
       boxplot (pressaotxt$SISTOLICA)
       # distribuição não é normal  
       
       summary (pressaotxt$SISTOLICA)
       describe (pressaotxt$SISTOLICA)       
       shapiro.test(pressaotxt$SISTOLICA)
       
    # Diastolica
       
       summary (freqlist(~ pressaotxt$DIASTOLICA))
       
       hist(pressaotxt$DIASTOLICA)
       boxplot (pressaotxt$DIASTOLICA)
       # distribuição não é normal  
       
       summary (pressaotxt$DIASTOLICA)
       describe (pressaotxt$DIASTOLICA)  
       shapiro.test(pressaotxt$DIASTOLICA)
       
       
      
  ##Altura é diferente entre os sexos? 
          # grafivos
       
       str (pressaotxt)
       # sexo esta como variável numerica tem que criar uma categorica
       
       pressaotxt$sexocat <- as.factor(pressaotxt$SEXO)
       summary (freqlist(~ pressaotxt$sexocat))
       
       boxplot(pressaotxt$alturam ~ pressaotxt$sexocat) 
      
        str(pressaotxt$sexocat) 
       
       
       describe(pressaotxt$alturam)
       # variancias sao iguais? 
       var.test(alturam ~ sexocat, alternative='two.sided', conf.level=.95, 
                data= pressaotxt) # p = 0.4037 
       
       t.test(pressaotxt$alturam ~ pressaotxt$sexocat)#default Welch´s test; bicaudal
       t.test(pressaotxt$alturam ~ pressaotxt$sexocat, var.equal = TRUE)
       
        # nao parametrico
       wilcox.test(pressaotxt$alturam ~ pressaotxt$sexocat, exact = FALSE)

         library (ggplot2)
         ggplot(data = pressaotxt,
          aes(x= sexocat, y = alturam, fill= sexocat )) + 
           geom_boxplot() +
           xlab ("Sexo") +
           ylab ("Altura") 
         summary(tableby(sexocat ~ alturam ,  data = pressaotxt))
         
  ##Peso é diferente entre os sexos? 
         # graficos
         
         
         hist (pressaotxt$pesokg) 
         boxplot(pressaotxt$pesokg ~ pressaotxt$sexocat) 
        
         
         library (ggplot2)
         ggplot(data = pressaotxt,
                aes(x= sexocat, y = pesokg, fill= sexocat )) + 
           geom_boxplot() +
           xlab ("Sexo") +
           ylab ("Peso(Kg)") 
         
         var.test(pesokg ~ sexocat, alternative='two.sided', conf.level=.95, 
                  data= pressaotxt) # p = 0.3528 
         
       
       t.test(pressaotxt$pesokg ~ pressaotxt$sexocat, var.equal = TRUE)
       
       wilcox.test(pressaotxt$pesokg ~ pressaotxt$sexocat, exact = FALSE)
       
       
       summary(tableby(sexocat ~ pesokg ,  data = pressaotxt), text= TRUE)
       
         
         
# IDADe é diferente entre os sexos? 
       
       boxplot(pressaotxt$IDADE ~ pressaotxt$sexocat) 
       hist (pressaotxt$IDADE ~ pressaotxt$sexocat) 
       
       var.test(IDADE ~ sexocat, alternative='two.sided', conf.level=.95, 
                data= pressaotxt)
       
       t.test (pressaotxt$IDADE ~pressaotxt$sexocat, var.equal = TRUE)
       wilcox.test(pressaotxt$IDADE ~ pressaotxt$sexocat, exact = FALSE)
       
       
       summary(tableby(sexocat ~ IDADE ,  data = pressaotxt))
       
       
       
       tabelaexemplo <- tableby(sexocat~ alturam  + IDADE + pesokg, 
          numeric.stats=c("median","q1q3", "mean", "sd"), data=pressaotxt)
       
       summary (tabelaexemplo, text = TRUE)
       
       
       
      #  Histograma com dois grupos sobrepostos. 
       
       library (ggplot2)
       
       ggplot (data = pressaotxt,  aes (x = SISTOLICA, fill = sexocat)) + geom_histogram(aes (color= sexocat))
      
       
       
       #histograma 
         ggplot (pressaotxt, aes(x= alturam)) +
         geom_histogram(position = 'identity', bins =30, alpha = 0.7) + facet_wrap(~ sexocat)
       
       
       
       ## histograma com superposicao para homens e mulheres
       ggplot (pressaotxt, aes(x= alturam , fill= sexocat)) +
       geom_histogram (position  = 'identity', bins = 30, alpha = 0.7) 
       
       
       wilcox.test(pressaotxt$alturam ~ pressaotxt$sexocat, exact = FALSE) # visualizar resultados no console
       # resultado nao parametrico p = 0.01354 - ver tb as diferencas de medianas 
       
       
       boxplot(pressaotxt$IDADE ~ pressaotxt$sexocat) 
       
       library (ggplot2)
       ggplot(data = pressaotxt,
              aes(x= sexocat, y = IDADE, fill= sexocat )) + 
         geom_boxplot() +
         xlab ("Sexo") +
         ylab ("Idade(anos)") 
       
       var.test(IDADE ~ sexocat, alternative='two.sided', conf.level=.95, 
                data= pressaotxt) # p = 0.3528 
       
       t.test(pressaotxt$IDADE ~ pressaotxt$sexocat, var.equal = TRUE)
       ttestidade <- t.test(pressaotxt$IDADE ~ pressaotxt$sexocat, var.equal = TRUE)
         
         
#Pressao sistolica é diferente entre os sexos? não está no exercicio mas vamos verificar
       
       
  ## Altura e peso - correlacao
       
       
       ggplot(pressaotxt, aes(alturam, pesokg)) +
         geom_point() + geom_smooth(method = "lm", 
                                    se = FALSE)
       
       
      #grafico de dispersao por sexo - apenas de curiosidade
       
        ggplot(data = subset(pressaotxt, !is.na (sexocat)), 
              aes(x = alturam, y = pesokg)) +
         geom_point(aes(shape = sexocat)) +
         geom_point(aes(color = sexocat)) +
         geom_smooth(method = "lm", 
                     se = FALSE, 
                     aes(color = sexocat)) + 
         labs(title = "Gráfico de dispersao altura vs peso por sexo",
              x = "Altura ",
              y = "Peso (Kg)")
       
                    #apenas de curiosidade separando subset 
                   ggplot(subset (pressaotxt, sexocat == "1"),
                          aes(x = alturam, y = pesokg)) +
                     geom_point(aes(shape = sexocat)) +
                     geom_point(aes(color = sexocat)) +
                     geom_smooth(method = "lm", 
                                 se = FALSE, 
                                 aes(color = sexocat)) + 
                     labs(title = "Gráfico de dispersao altura vs peso por sexo",
                          x = "Altura ",
                          y = "Peso (Kg)")
                   
         
                   
        # ANALISE DE REGRESSÃO LINEAR    
                    
                    
         regaltpeso <- lm( pressaotxt$alturam ~ pressaotxt$pesokg) # regressao com duas variaveis
         summary (regaltpeso) 
        
         residuoaltpeso = resid(regaltpeso)
          
         plot (pressaotxt$pesokg, residuoaltpeso) 
         abline (0,0)

                  
         par(mfrow=c(2,2))
         
         
         plot (regaltpeso, which = 1:6)
         
        
         summary (influence.measures(regaltpeso))
#Podemos tambem  pedir alem dos graficos as medidas influenciadoras 
         
         
         
         
         
  #ALTURA E PRESSAO SISTOLICA
  
         summary (freqlist(~ pressaotxt$sexocat))
  
          
    ggplot(pressaotxt, aes(alturam, SISTOLICA)) +
            geom_point() + geom_smooth(method = "lm", 
                                       se = FALSE)
          
          
          #grafico de dispersao por sexo - apenas de curiosidade
          
          ggplot(data = subset(pressaotxt, !is.na (sexocat)), 
                 aes(x = alturam, y = SISTOLICA)) +
            geom_point(aes(shape = sexocat)) +
            geom_point(aes(color = sexocat)) +
            geom_smooth(method = "lm", 
                        se = FALSE, 
                        aes(color = sexocat)) + 
            labs(title = "Gráfico de dispersao altura vs pressao sistólica por sexo",
                 x = "Altura ",
                 y = "Pressão Sistolica (mmHg)")
          
          
              # Analise de regressoa linear       
  
              
              
              regaltsist <- lm( pressaotxt$SISTOLICA ~ pressaotxt$alturam) # regressao com duas variaveis
              summary (regaltsist) 
              
                    residuoaltsist = resid(regaltsist)
                    
                    plot (pressaotxt$alturam, residuoaltsist) 
                    abline (0,0)
                    
                    
                    par(mfrow=c(2,2))
                    
                    plot (regaltsist, which = 1:6)
                    
                    summary (influence.measures(regaltsist))
                              #Podemos tambem  pedir alem dos graficos as medidas influenciadoras. 
                               # valores 6 e 42
                        
                                  ## removendo valores. 
                                pressaotxt$SISTOLICA [pressaotxt$OBS == 6] <- NA
                                pressaotxt$SISTOLICA [pressaotxt$OBS == 42] <- NA
                                
                                ggplot(pressaotxt, aes(alturam, SISTOLICA)) +
                                  geom_point() + geom_smooth(method = "lm", 
                                                             se = FALSE)
                          
         
          ## adicionando sexo como interacao
                      
                      regaltsist2 <- lm( pressaotxt$SISTOLICA ~ pressaotxt$alturam + pressaotxt$sexocat +
                                        pressaotxt$alturam*pressaotxt$sexocat) # regressao com 3 variaveis
                      summary (regaltsist2) 
                      
                      
                    sistaltsexo1 <-  lm (pressaotxt$SISTOLICA ~ pressaotxt$alturam, subset = (pressaotxt$sexocat == '1'))
                    summary (sistaltsexo1)
                    anova (sistaltsexo1 ) 
                      
                    
                    sistaltsexo2 <-  lm (pressaotxt$SISTOLICA ~ pressaotxt$alturam, subset = (pressaotxt$sexocat == '0'))
                    summary (sistaltsexo2)
                    
                      
                      
            
         
## Exemplo da importancia dos resíduos. 
         
         
         milhas <- read.table (file = 'milhas-desgaste.txt', header = TRUE)
         
         ggplot(milhas, aes(milhas, desgaste)) +
           geom_point() + geom_smooth(method = "lm", 
                                      se = FALSE)
         
         milha <- lm( milhas$milhas ~ milhas$desgaste) # regressao com duas variaveis
         summary (milha)
         plot (milha)
         
         residmilhas = resid(milha)
         
         plot (milhas$milhas, residmilhas) 
         abline (0,0)
         


         
         
         
  # Pressão sistolica tem diferenca entre cor/etinia? 
         # ANALISE DE VARIANCIA 
         
        #1. verificar normalidade do Y 
             hist(pressaotxt$SISTOLICA)
             qqnorm(pressaotxt$SISTOLICA)
             qqline (pressaotxt$SISTOLICA)
             shapiro.test(pressaotxt$SISTOLICA)
              #admitimos normalidade 
             
         # 2.Boxplot para ter noçao 
                  boxplot(pressaotxt$alturam ~ pressaotxt$ETNIA)
                  
         # 3. Obter media e desvio para cada categoria 
         
        
         
         # 4. As variancias são iguais entre os grupos? 
                  
                 # pacote car baixar por causa do teste de levene
         library (car)
                                  pressaotxt$etniacat [pressaotxt$ETNIA == "A"] <- 1
                                  pressaotxt$etniacat [pressaotxt$ETNIA == "B"] <- 2
                                  pressaotxt$etniacat [pressaotxt$ETNIA == "C"] <- 3
                                  pressaotxt$etniacat <- as.factor(pressaotxt$etniacat)  
                                  
                                str (pressaotxt)   
                                
                  bartlett.test (data = pressaotxt, SISTOLICA ~ ETNIA) # muito sensivel a desvio de normalidade
                  
                  leveneTest(data = pressaotxt, SISTOLICA ~ ETNIA) # levene by defaut testa varibilidade ao redor da mediana
                  leveneTest(data = pressaotxt, SISTOLICA ~ etniacat) # levene by defaut testa varibilidade ao redor da mediana
                  
                  leveneTest(data = pressaotxt, SISTOLICA ~ etniacat, center = mean) # menos sensivel que o Barlett
                  
                  fligner.test (data = pressaotxt, SISTOLICA ~ etniacat) # teste não parametrico muito robusto
                  
                  
         
         
         #5. Anova (lm ou aov)
        
        ## lm 
        
                lm (SISTOLICA ~ etniacat, data = pressaotxt) # fornece apenas intercepto
                
                # as informações outras de um modelo que incluem a analise de variancia são guardadas no objeto 
                #modelo.anova que ao pedir summary retorna as informacoes gerais
                
                anova.sistolica <- lm (SISTOLICA ~ etniacat, data = pressaotxt)
                summary (anova.sistolica )      
                anova (anova.sistolica ) # este sim me fornece sum of sq etc. 
                # faz anova mas nao consegue aplicar tukey ##;
                
                
                ## aov
                anova.sistolica2 <- aov (SISTOLICA ~ etniacat, data = pressaotxt)
                summary (anova.sistolica2)
                
                # precisamos das médias de cada grupo e o erro padrão e depois fazer o pos teste, neste caso 
                # nao precisaria porque as medias não foram diferentes
                
                shapiro.test(resid(anova.sistolica2)) #
                # nessa situação os residuos normais. 
                # aqui usei o shapiro emobra nao seja apropriado pq os residuos nao sao independentes
                 
                ## pos testes
                
                TukeyHSD(x =anova.sistolica2,  conf.level=0.95)
                plot (TukeyHSD(anova.sistolica2,  conf.level=0.95))
               
                install.packages("DescTools")
                library(DescTools)
                
                # Scheffe´s test
                
                ScheffeTest(anova.sistolica2)
                
                pairwise.t.test(data = pressaotxt,  SISTOLICA, etniacat, p.adjust = none, pool.sd = T)
                
          # exemplo cafe de manha atenção 
                
                cafe <- read.table (file = 'cafe-anova-oneway.txt', header = TRUE)
                
                hist(cafe$atencao)
                qqnorm (cafe$atencao)
                qqline (cafe$atencao)
                
                anova.cafe <- aov (atencao ~ jejum, data = cafe)
                summary (anova.cafe)
                
               
                # deu errado... 
                str(cafe) # cafe inteiro
                cafe$jejum <- as.factor(cafe$jejum)
                
                anova.cafe <- aov (atencao ~ jejum, data = cafe)
                summary (anova.cafe)
                
                library (car)
              # igualdade das variancias 
                leveneTest(data = cafe, atencao ~ jejum)
                
              # Pos test  
               TukeyHSD(x =anova.cafe,  conf.level=0.95)
               plot (TukeyHSD(anova.cafe,  conf.level=0.95))
               
               pairwise.t.test(cafe$atencao, cafe$jejum, p.adjust="holm") 
                
               # Dunnet - comprar cada media apenas com um controle
               
               dunnet_comparison <- glht(anova.cafe, linfct = mcp(group = "Dunnett"))
               
                library(DescTools)
               DunnettTest(x=cafe$atencao, g=cafe$jejum)
               
               # kruskall-wallis nao parametrico
               
               kruskal.test(data = cafe, atencao ~ jejum)
               
               
                
         # 6. Anova twoway etnia e sexo
               fone <-  read.table (file = 'telefone.txt', header = TRUE)
               hist (fone$custo)
                
         
         
         
                
                
                
                
  ## Outro exemplo
         
         
         
         alt <- read.table (file = 'altura2022.txt', header = TRUE)
         
    #1. Verificar normalidade
         
         hist (alt$aluno)
         boxplot (alt$aluno)
         qqnorm(alt$aluno)
         qqline (alt$aluno)
         shapiro.test(alt$aluno)
         
         describe (alt$aluno)
    
    #2. Gafico de dispersão
         
         
         ggplot(alt, aes(pai , aluno)) +
           geom_point() + geom_smooth(method = "lm", 
                                      se = FALSE)
    #3. Regressao  
         altura <- lm( alt$aluno ~ alt$pai) # regressao com duas variaveis
         summary (altura)
        
    #4. Resíduos graficos
    
         
         plot (altura, which = 1:6)
         
         
         residaltura = resid(altura)
         
         plot (alt$mae, residaltura)
         abline (0,0)
         
         ## residuo identificando meninos e meninas
         ggplot(alt, aes(pai , residaltura , color = sexocat)) +
           geom_point() 
         
             
             
    
    #5. Observações influentes
         
         summary (influence.measures(altura))
         
         
         # vamos remover o 68 e 0 72.
         
         alt$pai [alt$indi == 72] <- NA
         
         ggplot(alt, aes(pai , aluno)) +
           geom_point() + geom_smooth(method = "lm", 
                                      se = FALSE)
         
         
         alt$pai [alt$indi == 68] <- NA
         
         ggplot(alt, aes(pai , aluno)) +
           geom_point() + geom_smooth(method = "lm", 
                                      se = FALSE)
         
         
         altura <- lm( alt$aluno ~ alt$pai) # 
         summary (altura)
         
         
      # Altura do aluno e pai por sexo
         
         str(alt)
         alt$sexocat <- as.factor(alt$sexo)
         
         ggplot(data = subset(alt, !is.na (sexocat)), 
                aes(x = aluno, y = pai)) +
           geom_point(aes(shape = sexocat)) +
           geom_point(aes(color = sexocat)) +
           geom_smooth(method = "lm", 
                       se = FALSE, 
                       aes(color = sexocat)) + 
           labs(title = "Gráfico de dispersao altura do aluno vs altura do pai por sexo",
                x = "Altura (PAi) ",
                y = "Altura(aluno")
         
         
         
         altsexo <- lm( alt$aluno ~ alt$pai + alt$sexocat + alt$sexocat*alt$pai)
         summary (altsexo)
         
         
         
         altsexob <- lm( alt$aluno ~ alt$pai + alt$sexocat )
         summary (altsexob)
         
         plot (altsexo, which = 1:6)

         
         residaltsexob = resid(altsexob)
         
         plot (alt$pai, residaltsexob)
         abline (0,0)
 
         ggplot(alt, aes(pai , residaltsexob , color = sexocat)) +
           geom_point() 
         
         
         
 
         
         altsexoc <- lm( alt$aluno ~ alt$pai + alt$sexocat + alt$mae)
         summary (altsexoc)
         
         plot (altsexo, which = 1:6)
         
         
         residaltsexob = resid(altsexob)
         
         plot (alt$pai, residaltsexob)
         abline (0,0)
         
         
         
         
         p <- ggplot(pressaotxt, aes(x=ETNIA, y= SISTOLICA)) + 
           geom_dotplot(binaxis='y', stackdir  = "center")
        
         p
         
         p <- ggplot(pressaotxt, aes(x=etniacat, y= SISTOLICA, fill = etniacat)) + 
           geom_dotplot(binaxis='y', stackdir='center',
                        stackratio = 1, dotsize = 0.9, )
         p
         
         # Rotate the dot plot
         p + coord_flip()
         
 
         
         p + stat_summary(fun.y=mean, geom="point", shape=18,
                          size=3, color="red")
         # dot plot with median points
         p + stat_summary(fun.y=median, geom="point", shape=18,
                          size=3, color="orange")

         
         p + scale_x_discrete(limits=c("2", "3", "1"))