---
# System prepended metadata

title: 'Estatística, lista 01'

---

# 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:

```r=
set.seed(123)
```
## Problema 01.
![](https://i.imgur.com/TIgVGw5.png)

### 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.


```r=
## 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)
```

![](https://i.imgur.com/6rpLhYK.png)


### 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.
:::

![](https://i.imgur.com/gvzi0NC.png)


## Problema 02.
![](https://i.imgur.com/xpLsdlw.png)

### 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

```r=
### 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.

```r=
## 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.
![](https://i.imgur.com/TOxEloA.png)

### 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.

```r=
## 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?

```r=
## 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.

```r
## 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.
![](https://i.imgur.com/W7jTnEX.png)

### Algoritmo explícito

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

![](https://i.imgur.com/fCtirxZ.png)
![](https://i.imgur.com/nHxd8ig.png)


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

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

**Comprimento:**
(72.85579;72.88334)

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

### Código:

```r=
### 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.
![](https://i.imgur.com/rbhmpJ8.png)

### 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.
:::

![](https://i.imgur.com/wWBf8k6.png)

```r=
## 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.
![](https://i.imgur.com/MuR23o4.png)

### 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:

![](https://i.imgur.com/Sig4ThA.png)

#### Plot de M4:

![](https://i.imgur.com/guSJP5Z.png)


### 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

![](https://i.imgur.com/ON82sWF.png)

### Resíduos:

![](https://i.imgur.com/exVa5rb.png)


### Código em R

```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?
