# Estatística, lista 01

[Índice: Estatística IBI5086](/8z457TaXRmKsDwO5iHG6-w)

Feito pra disciplina IBI-5086, em Agosto/2020 (lista 01).

Alunos:
- Maria Fernanda Zaneli Campanari (nº USP: 12099257)
- Rafael Correia da Silva (nº USP: 9786869)
- Marília Viana Albuquerque de Almeida (nº USP: 12084806)
- Leonardo Yoshida (nº USP: 10399769)
- Gabriella Berwig Möller (nº USP: 11387062)
- Isabella Blagiem de Campos (nº USP: 12099240)

---

[toc]

## Preparo do R.

Antes de rodar as análises, setamos um seed:

set.seed(123)

## Problema 01.

### Suposições

Para os cálculos abaixo, supõe-se:
- Homoscedasticidade (homogeneidade de variâncias);
- Distribuição normal da população;
- Indepêndencia das observações.

## Lista 1

### Valores brutos
values_a <- c(89.7, 81.4, 84.5,
              84.8, 87.3, 79.7,
              85.1, 81.7, 83.7,
              84.5)

values_b <- c(84.7, 86.1, 83.2,
              91.9, 86.3, 79.3,
              82.6, 89.1, 83.7,
              88.5)

  
### Médias
mean_a <- mean(values_a)
mean_b <- mean(values_b)

### Aparecer
print(c(mean_a, mean_b))

### Desvios-padrões
sd_a <- sd(values_a)
sd_b <- sd(values_b)

### Intervalos
interval_a = c((mean_a - 2*sd_a), (mean_a + 2*sd_a))
interval_b = c((mean_b - 2*sd_b), (mean_b + 2*sd_b))

### Exibe
print(interval_a)
print(interval_b)

### Veja o intervalo
par(mfrow=c(1, 2))
hist(values_a)
hist(values_b)


### Resposta

$$
\bar{Y}_a = 84.24 \\
\bar{Y}_b = 85.54 \\
s_a = 2.9018 \\
s_b = 3.6503 \\
$$

Para $j = a$:

$$
(Y_a-2s_a;Y_a+2s_a) = \\
(84.24-2*2.90; 84.24+2*2.90) = \\
(78.44; 90.04)
$$

Para $j = b$:

$$
(Y_b-2s_b;Y_b+2s_b) = \\
(85.54-2*3.65; 85.54+2*3.65) = \\
(78.24; 92.84)
$$

### Discussão

- Interprete os dados.
- Qual é sua utilidade?
- Faça as suposições que achar necessárias.

:::warning
- O intervalo para $j = a$ encompassa todos os valores, sobrando ainda espaço pela margem esquerda;
- O mesmo é válido para é $j=b$;
- Neste conjunto, é util para enxergar até onde os dados vão tanto à esquerda quanto à direita.
- Em conjuntos com uma distribuição mais extrema, os valores irão ultrapassar estes limites, então esse número diz o quão "comprimidos" os dados estarão.
- Esse intervalo $j=a$ e $j=b$ equivale ao intervalo de variação dos valores do Grupo A e Grupo B. A sua utilidade é saber o comportamento dos valores quando adicionado os desvios padrão da amostra, para que possamos identificar se os valores encontrados na coleta de dados está dentro da margem esperada.
:::


## Problema 02.

### Suposições

Para os cálculos abaixo, supõe-se:
- Homoscedasticidade (homogeneidade de variâncias);
- Distribuição normal da população;
- Indepêndencia das observações.

#### Cálculo do primeiro intervalo

### Primeiro intervalo
#### Para A:
interval_2_a_p = mean_a - (2*sd_a / sqrt(length(values_a)))
interval_2_a_n = mean_a + (2*sd_a / sqrt(length(values_a)))
print(c(interval_2_a_p, interval_2_a_n))

#### Para B:
interval_2_b_p = mean_b - (2*sd_b / sqrt(length(values_b)))
interval_2_b_n = mean_b + (2*sd_b / sqrt(length(values_b)))
print(c(interval_2_b_p, interval_2_b_n))

#### Segundo intervalo, somente com os valores de A

Algoritmo:
- Calcular as variâncias;
- Tomar das duas variância o $s_\bar{d}$;

### Intervalos

- Intervalo 01. Seja $j = A$;

$$
(\hat{Y}_A -2s_A / \sqrt{n_A}; \hat{Y}_A +2s_A / \sqrt{n_A}) = \\
(84.24 - 5.80/3.16; 84.24 + 5.80/3.16) = \\
(82.40;86.08)
$$

- Intervalo 02. Seja $j = B$;

$$
(\hat{Y}_B -2s_B / \sqrt{n_B}; \hat{Y}_B +2s_B / \sqrt{n_B}) = \\
(85.54 - 7.30/3.16; 85.54 + 7.30/3.16) = \\
(83.23;87.84)
$$

A seguir, foi necessário determinar $\bar{d}$ e $s_\bar{d}$:

$$
\bar{d} = d_b - d_a \\
d = 85.54 - 84.24 \\
d = 1.3
$$

$$
s^2 = \frac{s_1^2*(n_1-1)+s_2^2*(n_2-1)}{n_1+n_2-2} = 83.13
$$

$$
s^2_\bar{d} = s^2(1/n_1+1/n_2) = 16.62
$$

$$
\sqrt{s^2_\bar{d}} = 4.07
$$

Para calcular, por fim, os intervalos (Intervalo 03):

$$
(1.3-2*4.07;1.3+2*4.07) = \\
(-6.855; 9.455)
$$

### Interpretação dos resultados

:::success
- Qual é a utilidade de cada um destes intervalos? Interprete-os. Com base no intervalo, decida se a média de B é maior que a de A.
- **Intervalos 1 e 2**: sempre que forem tomadas amostras aleatórias desta população, é provável que sua média esteja entre esses intervalos.
- **Intervalo 3**: é útil para observar o quanto variam ao máximo as diferenças de médias, com a observação de que, no nosso caso, a média de B é maior que a de A ($d_b > d_a$). Como o valor à direita é maior que o módulo do valor a esquerda, observamos também por este intervalo que a média de B parece ser maior que a de A.
:::

### Código completo do problema 2.

## Problema 2.

### Primeiro intervalo
#### Para A:
interval_2_a_p = mean_a - (2*sd_a / sqrt(length(values_a)))
interval_2_a_n = mean_a + (2*sd_a / sqrt(length(values_a)))
print(c(interval_2_a_p, interval_2_a_n))

#### Para B:
interval_2_b_p = mean_b - (2*sd_b / sqrt(length(values_b)))
interval_2_b_n = mean_b + (2*sd_b / sqrt(length(values_b)))
print(c(interval_2_b_p, interval_2_b_n))

### Segundo intervalo correto
#### Variância de A
#### Variância de B
var_a = sum((values_a - mean_a)*(values_a - mean_a))/(length(values_a)-1)
var_b = sum((values_b - mean_b)*(values_b - mean_b))/(length(values_b)-1)

### Re-calculo!
s_squared = (var_a*(rol_a-1))+(var_b*rol_b-1)/(rol_a+rol_a-2)
s_squared_d = s_squared * (1/rol_a + 1/rol_b)
s_d = sqrt(s_squared_d)

### E então, os intervalos
interval_2 = c((mean_b - mean_a)-(2*s_d), (mean_b - mean_a)+(2*s_d))
print(interval_2)

## Problema 03.

### Respostas

:::success
- $\bar{d} = -1.3$
- $\bar{d}_b = -3.66$
- $s_{\bar{d}b} = 0.54$
- $t_b = -4.38$
- Intervalo:
$$
(\bar{d}+t_{b2.5\%}s_\bar{d};\bar{d}-t_{b97.5\%}s_\bar{d} = X) = \\
(-1.565; -1.218)
$$

---
- Decida se a média do tratamento B é maior que a de A: I. Considerando que o valor original da média está incluso no intervalo; II. considerando que o zero não está no intervalo, por I. e II. concluímos que há diferença entre médias, e que a média de B ($\bar{Y}_B$) maior que a de A ($\bar{Y}_A$).
:::

### Subproblema 1. Cálculo dos valores.

## Problema 3.
### Cálculo dos valores
#### Diferença das médias
mean_diff = mean_a-mean_b
print(mean_diff)

#### Desvio padrão das diferenças
#### soma? de sd(a) e sd(b)
sd_diff = (sd(values_a)/10) + (sd(values_b)/10)
print(sd_diff)

### Subproblema 2. Bootstrap ou reamostragem?

## Re-amostragem e boostrap
revalues_a <- sample(values_a, size = 10, replace = TRUE)
revalues_b <- sample(values_b, size = 10, replace = TRUE)

### Médias e SD
mean_re_a = mean(revalues_a)
mean_re_b = mean(revalues_b)
sd_re_a = sd(revalues_a)
sd_re_b = sd(revalues_b)
mean_revalues = mean_re_a-mean_re_b
sd_revalues = sd_re_a/10 + sd_re_b/10

### Cálculo do T
t_bootstrap = (mean_revalues-mean_diff)/sd_revalues

### Subproblema 3. Bootstrap enorme.

## Re-amostrar mil vezes
resample_bootstrap = vector()

for (sample in 1:1000) {
  ## 1. Re-amostra A e B
  reval_a = sample(values_a, size = 10, replace = T)
  reval_b = sample(values_b, size = 10, replace = T)
  ## 2. Calcula as médias de A e B
  mean_reval_a = mean(reval_a)
  mean_reval_b = mean(reval_b)
  ## 3. Calcula os d.p. de A e B
  sd_reval_a = sd(reval_a)
  sd_reval_b = sd(reval_b)
  ## 4. Calcula o desvio padrão da diferença
  sd_diffs_a_b = sd_reval_a/10 + sd_reval_b/10
  ## 5. Calcula a diferença entre:
  ### - A nova diferença das médias, tirando a antiga
  diffs_mean = (mean_reval_a - mean_reval_b) - mean_diff
  # OBS: loop_bootstrap é a "estatística t"
  loop_bootstrap = diffs_mean/sd_diffs_a_b
  # Guarda os novos valores
  resample_bootstrap = c(loop_bootstrap, loop_bootstrap_a)
}

## Calcula os quantis
quantiles = quantile(x = resample_bootstrap,
                       probs = seq(0, 1, 0.025))

quantile_025 = quantiles["2.5%"]
quantile_975 = quantiles["97.5%"]

## Finalmente, os intervalos
ex3_interval_left = mean_diff - (quantile_025 * sd_diff)
ex3_interval_right = mean_diff + (quantile_975 * sd_diff)
print(c(ex3_interval_left, ex3_interval_right))

## Problema 04.

### Algoritmo explícito

[Referência](https://stattrek.com/estimation/confidence-interval.aspx)


### Subproblema 01. Intervalo de confiança.

:::success
**Intervalo:**
72.86 +- 0.0137

**Comprimento:**
(72.85579;72.88334)

**Margem de erro:**
0.01377
:::

### Código:

### Passos
### 0. Importar os dados (duh)
prob4_set = read.delim(file = "../learnerboi/06_USP_stat/data/L01/pulse.csv",
                       header = T,
                       sep = ";")

head(prob4_set)

### 1. Decidir a estatística (média)
data_vect_repouso <- prob4_set$P1
mean_rep = mean(data_vect_repouso)
data_rol = length(data_vect_repouso)

### 2. Nível de confiança
data_input_conf = 0.95

### 3. Margem de erro.
data_standard_error = sd(data_vect_repouso)/sqrt(data_rol))

### 4. Valor crítico
data_alpha = 1 - (data_input_conf/100)

### 5. Achar a probabilidade cumulativa/crítica
data_p_star = 1 - (data_alpha/2)

### 6. Achar o grau de liberdade.
data_df = data_rol - 1

### 7. Encontrar o valor crítico da estatística t.
#### OBS: olhado na distribuição T:
#### https://stattrek.com/online-calculator/t-distribution.aspx
#### Nosso cálculo: https://drive.google.com/uc?id=1kyZNVroezzA1YhuHAtV-9Vhen5pugrLY
data_critical_value = 0.012

### 8. Cálculo da margem de erro.
data_ME = data_critical_value * data_standard_error

### 9. Exibe o intervalo.
### média +- data_me

## Problema 05.

### Resposta:

::: success
- Em 06 simulações não foi observada a média verdadeira ($\mu$)
- Em 94 simulações a média verdadeira $\mu$ foi incluída.
:::


## Problema 5.
library(RcmdrPlugin.TeachingDemos)

### Código da professora
ci.examp(mean.sim = 72, sd = 11,
         n = 92, reps = 100,
         conf.level = 0.95, method = "t")

## Problema 06.

### Respostas

:::success
- Comente os resultados:
- **Valores dos resíduos**: a partir dos resultados, observamos o menor resíduo para o modelo M3, isto é, o ajuste dos dados foi o mais próximo da linha de regressão;
- **Significância dos coeficientes**: Todos rejeitam a hipótese nula ($H_0$ rejeitada); na ordem: $M1 < M2 < M3 < M4$;
- **Coeficiente de determinação** $R^2$: Os melhores modelos foram M3 e M4, ou seja, aqueles cujos dados ajustados foram mais proximamente ajustado ao modelo linear
- **Estatística F**: Considerando que o `F-value` é a razão entre o ajuste do modelo e o erro, e quanto maior mais ajustado é o modelo, ele favoreceu desproporcionalmente o modelo **M3**, que teve o maior valor de F. A ordem foi: $M3 > M4 > M1 > M2$
- Qual modelo você escolheria?
- R: Ficaríamos entre os modelos M3 e M4, tendendo ao modelo M3. Ambos tiveram valor similares de: (1) resíduos, (2) valores de significância e (3) $R^2$. Entretanto, o modelo M3 teve um maior valor de estatística F, o que indica que a variância entre as médias foi maior que a variância dentro dos tratamentos (erro), o que aponta para um ajuste melhor dos dados ($93.31 > 62.21$).
- Houve dúvida, porém, se o modelo de M4 não seria um ajuste melhor justamente por considerar uma interação entre: (1) a diferença entre a média da pessoa em repouso $P1c$ interagindo com o status de corrida $Ran$. Não houve consenso no grupo se biologicamente esse é um fator que explica ou não o comportamento dos batimentos.
:::

#### Plot de M3:


#### Plot de M4:


### Cálculo dos valores de interesse

- OBS: Ran é se correu ou não

| | M1 | M2 | M3 | M4 |
|--------------|:--------:|:--------:|:--------:|:--------:|
| $R^2$ | 0.379 | 0.3327 | **0.6771** | **0.6795** |
| residuals.se | 13.54 | 14.04 | **9.822** | **9.84** |
| p-val | 6.21E-11 | 1.76E-09 | **2.20E-16** | **2.20E-16** |
| F-value | 55.09 | 44.88 | **93.31** | 62.21 |
| P1c | 0.9567 | | 0.9124 | 0.8489 |
| Ran | NA | -20.19 | 19.1227 | 19.082 |
| P1cRan | NA | NA | NA | 0.157 |
| Intercept | 80 | 112.7 | 73 | 72.69 |

### Sobre os dados pedidos (observação nossa)

#### Estatística F


### Resíduos:


### Código em R

### Código da professora:
#### Meus ajustes
P1 <- prob4_set$P1
P2 <- prob4_set$P2
dados <- prob4_set

#Ajuste do modelo M1
# y = b0 + b1P1c + e
### Calcula o desvio da média
P1c <- P1 - mean(P1)
dados$P1c <- P1c

### Desvio da média em relação à P2
plot(P2 ~ P1c)

### Faz o modelo linear lm()
m1 <- lm(P2~P1c, dados)

### Descritivos do modelo 1
summary(m1)
names(m1)
m1$coefficients

### Interpretação do modelo 1
plot(m1$fitted.values,m1$residuals)
lines(lowess(m1$fitted.values,m1$residuals),col="red")
plot(P2~P1c)
abline(m1)

#Ajuste do modelo M2
#y=b0+b2Ran+e
Ran <- dados$Ran
plot(P2~Ran)
plot(P2~factor(Ran))
m2<-lm(P2~Ran)
summary(m2)
m2$coefficients #igual às médias dos grupos
mean(P2[Ran=="0"]) #intercept
mean(P2[Ran=="1"]) #intercept + Ran
plot(m2$fitted.values,m2$residuals)
lines(lowess(m2$fitted.values,m2$residuals),col="red")
plot(P2~Ran)
abline(m2)

#Ajuste do modelo M3
#y=b0+b1P1c+b2Ran+e
Ran
library(car)
Ran<-recode(Ran,"1='1';2='0'")
Ran
plot(P2~P1c, pch=23, bg=c('red', 'blue')[factor(Ran)])
m3<-lm(P2~P1c+Ran)
summary(m3)
m3$coefficients
plot(m3$fitted.values,m3$residuals)
lines(lowess(m3$fitted.values,m3$residuals),col="red")

b0<-m3$coefficients[1]
b0
names(b0)<-NULL
b0
b1<-m3$coefficients[2]
names(b1)<-NULL
b2<-m3$coefficients[3]
names(b2)<-NULL
plot(P2~P1c, pch=23, bg=c('red', 'blue')[factor(Ran)])
abline(b0,b1, col="red")
abline((b0+b2),b1,col="blue")

#Ajuste do modelo M4
#y=b0+b1P1c+b2Ran+b3(P1c*Ran)+e
P1cRan<-P1c*Ran
plot(P2~P1c, pch=23, bg=c('red', 'blue')[factor(Ran)])
m4<-lm(P2~P1c+Ran+P1cRan)
summary(m4)
m4$coefficients
plot(m4$fitted.values,m4$residuals)
lines(lowess(m4$fitted.values,m4$residuals),col="red")

b0<-m4$coefficients[1]
b0
names(b0)<-NULL
b0
b1<-m4$coefficients[2]
names(b1)<-NULL
b2<-m4$coefficients[3]
names(b2)<-NULL
b3<-m4$coefficients[4]
names(b3)<-NULL
plot(P2~P1c, pch=23, bg=c('red', 'blue')[factor(Ran)])
abline(b0,b1, col="red")
abline((b0+b2),(b1+b3),col="blue")

## Extraindo os valores
m1$coefficients
m2$coefficients
m3$coefficients
m4$coefficients
summary(m1)
summary(m2)
summary(m3)
summary(m4)


summary(m1)$r.squared
summary(m2)$r.squared
summary(m3)$r.squared
summary(m4)$r.squared

## Perguntas.

- Como se interpreta um intervalo? Fizemos direito?
- Qual é a fórmula do desvio padrão da diferença?
- No problema 3, é chamado $sd(a) + sd(b)$
- Na internet, é outra coisa completamente diferente
- Eu não acho que "bootstrap" do problema 3 seja de fato bootstrap.
- Interpretar valores de resíduos e coeficientes
- O intervalo 2 inclui o 0 ou é discreto?