tutorialzinho básico

##Load##

library (devtools)

library (dplyr)

library (ggplot2)

library (shiny)

library(statsr)

library(huxtable)

data("nlschools")

---

library(GGally)

esse é um tutorialzinho - básico- começa com regressão e depois segue com outros - se tiver ideias ou dúvidas me mande gregoriosilva1986@gmail.com

FONTE e Codebook

FONTE Snijders, T. A. B. and Bosker, R. J. (1999) Multilevel Analysis. An Introduction to Basic and Advanced Multilevel Modelling. London: Sage.

References Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth edition. Springer.

CODEBOOOK: lang language test score.

IQ verbal IQ.

class class ID.

GS class size: number of eighth-grade pupils recorded in the class (there may be others: see COMB, and some may have been omitted with missing values).

SES social-economic status of pupil’s family.

COMB were the pupils taught in a multi-grade class (0/1)? Classes which contained pupils from grades 7 and 8 are coded 1, but only eighth-graders were tested.

r str(nlschools)

'data.frame': 2287 obs. of 6 variables: $ lang : int 46 45 33 46 20 30 30 57 36 36 ... $ IQ : num 15 14.5 9.5 11 8 9.5 9.5 13 9.5 11 ... $ class: Factor w/ 133 levels "180","280","1082",..: 1 1 1 1 1 1 1 1 1 1 ... $ GS : int 29 29 29 29 29 29 29 29 29 29 ... $ SES : int 23 10 15 23 10 10 23 10 13 15 ... $ COMB : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...

O que será feito (ou o indíce)

A variável a ser explicada é a nota no teste de linguaguem(lang). Para cumprir tal propósito essa atividade terá as seguintes atribuições:

1- Selecionar em uma novo “quadro de dados” chamado(q1) contendo as variáveis: lang,IQ, e SES.

2- Primeiro serão testes de dispersão e frequência das 3 variáveis

3- Depois associar-se-à a variável depedente com as outras.

4- Será rodada regressão linear simples **lang ~ IQ" e testes de diagnóstico

5- Será rodada regressão linear simples **lang ~ SES" e testes de diagnóstico

6- Será rodada regressão linear múltipla completa

7- Decidirá qual das hipóteses é a correta

#H1 - Há associação entre lang e IQ

There were 22 warnings (use warnings() to see them)

#H2 - Há associação entre lang e SES

#H3 - Há associção entre lang e ambas

#h0 - não há associação.
glance(q1)

Data frame tidiers are deprecated and will be removed in an upcoming release of broom.`data_frame()` is deprecated as of tibble 1.1.0.
Please use `tibble()` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

daí em diante outros tutoriais para mais info ver depois do sete consta novo indíce

8, 9, 10, ….. veja lá. outros testes. ‘sampling’ ‘intervalo de confiança’ ‘níveis de confiança’ ‘inferência de variável numérica sem ser regressão’ daí em diante outros tutoriais para mais info ver depois do sete consta novo indíce

na sessão 11 - encontra testes de associação entre numérica e caractere bem legal.

---

1- Criando o quadro de dados

q1 <- subset(nlschools, select = c(lang, IQ, SES)) %>% na.omit()
dim(nlschools)

[1] 2287    6

dim(q1)

[1] 2287    3

Ok não houve perda de casos - Q1 criado.

---

2- testes de dispersão e frequência das 3 variáveis

summary(q1)

      lang             IQ             SES       
 Min.   : 9.00   Min.   : 4.00   Min.   :10.00  
 1st Qu.:35.00   1st Qu.:10.50   1st Qu.:20.00  
 Median :42.00   Median :12.00   Median :27.00  
 Mean   :40.93   Mean   :11.83   Mean   :27.81  
 3rd Qu.:48.00   3rd Qu.:13.00   3rd Qu.:35.00  
 Max.   :58.00   Max.   :18.00   Max.   :50.00  

hist(q1$lang)

hist(q1$IQ)

hist(q1$SES)

---

3- Depois associar-se-à a variável depedente com as outras.

3.1 lang e IQ

plot(q1$lang ~ q1$IQ)

q1 %>%
  summarise(cor(lang, IQ))

obs - tanto visual, quanto no valor de cor, parece haver associação.

3.2 lang e SES

plot(q1$lang ~ q1$SES)

q1 %>%
  summarise(cor(lang, SES))

Não parece haver associação (ou bem fraca)

---

4- Será rodada regressão linear simples com lang ~ IQ" e testes de diagnóstico

m1 <- lm(lang ~ IQ, data = q1)
summary(m1)

Call:
lm(formula = lang ~ IQ, data = q1)

Residuals:
     Min       1Q   Median       3Q      Max 
-28.7022  -4.3944   0.6056   5.2595  26.2212 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.52848    0.86682   10.99   <2e-16 ***
IQ           2.65390    0.07215   36.78   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.137 on 2285 degrees of freedom
Multiple R-squared:  0.3719,    Adjusted R-squared:  0.3716 
F-statistic:  1353 on 1 and 2285 DF,  p-value: < 2.2e-16

Bom resultado, - cada indicar de Qi parece fazer subir 2.6 em lang.

vejamos a análise

ggplot(data = q1, aes(x = IQ, y = lang)) +
  geom_point() +
  stat_smooth(method = "lm", se = FALSE)

Diagnóstico

#Linearidade: você já verificou se a relação entre as corridas e as batidas é linear usando um gráfico de dispersão. Devemos também verificar esta condição com um gráfico dos valores residuais vs. ajustados (previstos).
ggplot(data = m1, aes(x = .fitted, y = .resid)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  xlab("Fitted values") +
  ylab("Residuals")

#Resíduos quase normais: para verificar esta condição, podemos olhar um histograma
ggplot(data = m1, aes(x = .resid)) +
  geom_histogram(binwidth = 25) +
  xlab("Residuals")

ggplot(data = m1, aes(sample = .resid)) +
  stat_qq()

Em suma, parece bom. O R_adju_m1 é 0.3716

r_adju_m1 <- 0.3716

---

5- Será rodada regressão linear simples com lang ~ SES" e testes de diagnóstico

m2 <- lm(lang ~ SES, data = q1)
summary(m2)

Call:
lm(formula = lang ~ SES, data = q1)

Residuals:
    Min      1Q  Median      3Q     Max 
-32.577  -5.483   0.557   6.363  21.410 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 32.77758    0.48219   67.98   <2e-16 ***
SES          0.29330    0.01614   18.17   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.418 on 2285 degrees of freedom
Multiple R-squared:  0.1263,    Adjusted R-squared:  0.1259 
F-statistic: 330.2 on 1 and 2285 DF,  p-value: < 2.2e-16

A variável SES parece ter menos poder explicativo, embora apresente algum. o valor de r adjusted é menor que na sessão anterior, mas não é desprezível.

veja a linha:

ggplot(data = q1, aes(x = SES, y = lang)) +
  geom_point() +
  stat_smooth(method = "lm", se = FALSE)

testes de diagnóstico dá medo

#Linearidade: você já verificou se a relação entre as corridas e as batidas é linear usando um gráfico de dispersão. Devemos também verificar esta condição com um gráfico dos valores residuais vs. ajustados (previstos).
ggplot(data = m2, aes(x = .fitted, y = .resid)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  xlab("Fitted values") +
  ylab("Residuals")

#Resíduos quase normais: para verificar esta condição, podemos olhar um histograma
ggplot(data = m2, aes(x = .resid)) +
  geom_histogram(binwidth = 25) +
  xlab("Residuals")

ggplot(data = m1, aes(sample = .resid)) +
  stat_qq()

Em suma:

testes meia boca

e

valor de r adjusted para m2 é 0.125

r_adju_m2 <- 0.125

---

6- Será rodada regressão linear múltipla completa*

A primeira coisa é avaliar a função gg_pairs para ver se há colinearidade

ggpairs(q1)

Não há correlação entre as variáveis indepedentes a ponto do risco de colinearidade

AGORA A REGRESSÃO CHAMADA DE m3

m3 <- lm(lang ~ IQ + SES, data = q1)
summary(m3)

Call:
lm(formula = lang ~ IQ + SES, data = q1)

Residuals:
     Min       1Q   Median       3Q      Max 
-27.8208  -4.5375   0.4499   4.9173  25.6117 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  8.33128    0.85416   9.754   <2e-16 ***
IQ           2.40542    0.07430  32.376   <2e-16 ***
SES          0.14877    0.01409  10.557   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.971 on 2284 degrees of freedom
Multiple R-squared:  0.4011,    Adjusted R-squared:  0.4006 
F-statistic: 764.8 on 2 and 2284 DF,  p-value: < 2.2e-16

Adjusted de 0,40 melhor que o de m1- portanto o melhor modelo: Logo: parece ser h3, mas antes vamos aos testes de diagnóstico

dignóstico de m3(Reg.Linear Múltipla)

#passo 1 - LINEAR RELATIONSHIPS BETWEEN X AND Y - obs tem que ser random para dar boa
cog_final = lm(lang ~ IQ + SES, data = q1)
plot(cog_final$residuals ~ q1$IQ)

plot(cog_final$residuals ~ q1$SES)

#passo 2 - nearly normal residuals
hist(cog_final$residuals)

qqnorm(cog_final$residuals)
qqline(cog_final$residuals)

#passo 3 - constant variability of residualsa
plot(cog_final$residuals ~ cog_final$fitted)

plot(abs(cog_final$residuals) ~ cog_final$fitted)

#PASSO 4 - iNDEPENDENCE OF RESIDUALS
plot(cog_final$residuals)

Vou olhar em C:R _DUKE33 (procurar no meu Dropbox se estiver em outro PC) imagem a28 em diante para ver se está correto

---

7- qual das hipóteses é a correta?

H3 - Há associção entre lang e ambas

FIm??????

se quer ver como faz regressão só com variáveis numéricas okm

o Novo indíce

8 - brincando de sampling lang

9- intervalo de confiança sampling com lang

10- Confidence levels em lang

11- ‘inferência de variável numérica sem ser regressão’

---

8- brincando de sampling lang

#os dados "reais"
q1 %>%
  summarise(mu = mean(lang), pop_med = median(lang), 
            sigma = sd(lang), pop_iqr = IQR(lang),
            pop_min = min(lang), pop_max = max(lang),
            pop_q1 = quantile(lang, 0.25),  # first quartile, 25th percentile
            pop_q3 = quantile(lang, 0.75))  # third quartile, 75th percentile

##criando amostrinha de 50
samp1 <- q1 %>%
  sample_n(size = 50)
samp1 %>%
  summarise(mean_samp_1 = mean(lang), samp_1_med = median(lang), 
            s_samp_1 = sd(lang), samp_1_iqr = IQR(lang),
            samp_1_min = min(lang), samp1_max = max(lang),
            samp_1_q1 = quantile(lang, 0.25),  # first quartile, 25th percentile
            samp_1_q3 = quantile(lang, 0.75))  # third quartile, 75th percentile

esse amostrinha(nome real = samp1) aleatória de 50 deu parecida com o dado real, a maior diferença é o minimo de 9 na “população” e 21 na amostrinha

#Se estivermos interessados em estimar a lang média em Ames usando a amostrinha, nossa melhor estimativa é a média da amostra.acima 41,4. vou chamála de x_bar - vide abaixo

There were 21 warnings (use warnings() to see them)

samp1 %>%
  summarise(x_bar = mean(lang))

NA

#Aqui iremos gerar 15.000 amostras e calcular a média da amostra de cada um. Observe que estamos amostrando com substituição, `substituir = TRUE`, uma vez que as distribuições de amostragem são construídas com amostragem com substituição.
sample_means50 <- q1 %>%
                    rep_sample_n(size = 50, reps = 15000, replace = TRUE) %>%
                    summarise(x_bar = mean(lang))

`summarise()` ungrouping output (override with `.groups` argument)

ggplot(data = sample_means50, aes(x = x_bar)) +
  geom_histogram(binwidth = 1)

acima fiz módulo 2 week 1———daí módulo módulo 2 weeek 2 Foundations for inference - Confidence intervals…

---

9 - intervalos de confiança com lang repito que

acima fiz módulo 2 week 1

daí módulo módulo 2 weeek 2 Foundations for inference - Confidence intervals

então vemmmmmmmmmmmm

#começa criando uma de 60
n <- 60
samp <- sample_n(q1, n)
#aí se mostra as medidas de dispersão dessa amostrinha aleatória de 60(como fiz com o de 50 na sessão anterior)
#obs - para comparar com os dados "populacionais" e a amostrinha de 50 (chamada de samp1) basta ir a sessão anterior(8).
samp %>%
summarise (mean_samp60 = mean(lang), med_s60 = median(lang),
se_s60 = sd(lang))

confidence intervals de 95%

primeiro aplico a fórmulo

xbar +ou- z* s/raiz quadrada de n

a famosa fórmula para esse cáclulo

para saber o valor de z* ´é só olhar a famosa tabela

#Podemos encontrar o valor crítico para um intervalo de confiança de 95% usando
z_star_95 <- 1.959964
z_star_95

[1] 1.959964

aplica-se a fórmula:

samp %>%
  summarise(lower = mean(lang) - z_star_95 * (sd(lang) / sqrt(n)),
            upper = mean(lang) + z_star_95 * (sd(lang) / sqrt(n)))

---

10- níveis de confiança em lang

Confidence levels…………

#params é o nome dado para fixar o valor real da média populacional
params <- q1 %>%
  summarise(mu = mean(lang))
params

vamo continua trabaiando com a nossa amostrinha de 60 (samp)

antes só uma observação - o intervalo gerado acima não capta a média “real”(mu)

por isso ñ é 100%

né

vamos então criar vários intervalos de confiança…………………

ci <- q1 %>%
        rep_sample_n(size = n, reps = 50, replace = TRUE) %>%
        summarise(lower = mean(lang) - z_star_95 * (sd(lang) / sqrt(n)),
                  upper = mean(lang) + z_star_95 * (sd(lang) / sqrt(n)))

`summarise()` ungrouping output (override with `.groups` argument)

#eu fiz cinquenta, vamos ver os cinco primeiros
ci %>%
  slice(1:5)

todos os cinco acima captam o valor real mu - conforme se vê abaixo

ci <- ci %>%
  mutate(capture_mu = ifelse(lower < params$mu & upper > params$mu, "yes", "no"))
ci %>%
  slice(1:5)

vamos ver da linha 21 a 32

ci %>%
  slice(21:32)

como se vê acima a linha 25 (vigésima quinta simulação não contemplou o valor real de mu, as outras sim)

A seguir gerarei um ótimo gráfico com todas as 50

ci_data <- data.frame(ci_id = c(1:50, 1:50),
                      ci_bounds = c(ci$lower, ci$upper),
                      capture_mu = c(ci$capture_mu, ci$capture_mu))
ggplot(data = ci_data, aes(x = ci_bounds, y = ci_id, 
                           group = ci_id, color = capture_mu)) +
  geom_point(size = 2) +  # add points at the ends, size = 2
  geom_line() +           # connect with lines
  geom_vline(xintercept = params$mu, color = "darkgray") # draw vertical line

somente duas não capturaram o valor real. nossa amostra de 60 seria um bom número para 95% de confiança

---

11- ‘inferência de variável numérica sem ser regressão’

cserá necessária outra base - uma que tenha - uma variável numérica sendo explicada por uma categorial

conforme exemplo da módulo 2 week 3 - (peso do bebê explicado pelo hábito de fumar da gestante[sim ou não dicotômica])

para esse exercício procurarei outra relaçaõ simular (com tipos de variáveis parecidas) criarei q11 com as duas variáveis da brincadeira

a idade da mãe gestante (mage) pelo fato de ela ser casada ou não (marital - sim ou não dicotômica) no dia do nascimento da criança

acredito que a idade da mãe seja explicada por esse estado marital

omiti os missings

veja tudo abaixo

data(nc)
q11 <- subset(nc, select = c(mage,marital)) %>% na.omit()
dim(nc)

[1] 1000   13

dim(q11)

[1] 999   2

perdi um caso que era missing, de boas , ainda tem 999 na amostra.

summary(q11$mage)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  13.00   22.00   27.00   26.99   32.00   50.00 

summary(q11$marital)

    married not married 
        386         613 

hist(q11$mage)

plot(q11$marital)

primeiro teste de associação é o boxplot (visual)

boxplot(q11$mage ~ q11$marital, col = "orange", main="idade da mãe e casamento", ylab="idade da mãe", xlab="estado civil")

as não casadas parecem mães mais velhas que as casadas. nossa!

Os boxplots mostram como as medianas das duas distribuições se comparam, mas também podemos comparar as médias das distribuições usando o seguinte para primeiro agrupar os dados pela variável de marital e, em seguida, calcular a idade média nesses grupos usando a função média.

q11 %>%
  group_by(marital) %>%
  summarise(mean_mage = mean(mage))

`summarise()` ungrouping output (override with `.groups` argument)

as casadas tem o filho com 23 anos em média, ao passado que as não casados são mães com cerca de 29 anos.

Inference

inference(y = mage, x = marital, data = q11, statistic = "mean", type = "ht", null = 0, 
          alternative = "twosided", method = "theoretical")

Response variable: numerical
Explanatory variable: categorical (2 levels) 
n_married = 386, y_bar_married = 23.5622, s_married = 5.658
n_not married = 613, y_bar_not married = 29.1419, s_not married = 5.524
H0: mu_married =  mu_not married
HA: mu_married != mu_not married
t = -15.3164, df = 385
p_value = < 0.0001

Let’s pause for a moment to go through the arguments of this custom function. The first argument is y, which is the response variable that we are interested in: weight. The second argument is the explanatory variable, x, which is the variable that splits the data into two groups, smokers and non-smokers: habit. The third argument, data, is the data frame these variables are stored in. Next is statistic, which is the sample statistic we’re using, or similarly, the population parameter we’re estimating. In future labs we can also work with “median” and “proportion”. Next we decide on the type of inference we want: a hypothesis test (“ht”) or a confidence interval (“ci”). When performing a hypothesis test, we also need to supply the null value, which in this case is 0, since the null hypothesis sets the two population means equal to each other. The alternative hypothesis can be “less”, “greater”, or “twosided”. Lastly, the method of inference can be “theoretical” or “simulation” based.

For more information on the inference function see the help file with ?inference.

Exercise: What is the conclusion of the hypothesis test?

REJECT the NULL HIPOTHESIS p value bom

---

constituindo um intervalo de confiança (seguindo com dados e interação de cima)

inference(y = mage, x = marital, data = q11, statistic = "mean", type = "ci", method = "theoretical", order = c("not married","married"))

Response variable: numerical, Explanatory variable: categorical (2 levels)
n_not married = 613, y_bar_not married = 29.1419, s_not married = 5.524
n_married = 386, y_bar_married = 23.5622, s_married = 5.658
95% CI (not married - married): (4.8635 , 6.296)

CONCLUSÃO : Estamos 95% confiantes de que mulheres não casadas que tem seus filhos são em média 4,8 a 6,2 anos mais velhas que as casadas.

xxx