---
title: "Power"
author: "Prof. Eric A. Suess"
format: html
---

## CI Simulation

```{r}
library(BSDA)
```


```{r}
CIsim(100, 30, 100, 10)
    # Simulates 100 samples of size 30 from 
    # a normal distribution with mean 100
    # and standard deviation 10.  From the
    # 100 simulated samples, 95% confidence
    # intervals for the Mean are constructed 
    # and depicted in the graph. 

CIsim(100, 30, 100, 10, type="Var")
    # Simulates 100 samples of size 30 from 
    # a normal distribution with mean 100
    # and standard deviation 10.  From the
    # 100 simulated samples, 95% confidence
    # intervals for the variance are constructed 
    # and depicted in the graph.
    
CIsim(100, 50, .5, type="Pi", conf.level=.90)     
    # Simulates 100 samples of size 50 from 
    # a binomial distribution where the population
    # proportion of successes is 0.5.  From the
    # 100 simulated samples, 90% confidence
    # intervals for Pi are constructed 
    # and depicted in the graph.  
```

## Power

```{r}
library(stats)

power.t.test(n = 20, delta = 1)
power.t.test(power = .90, delta = 1)
power.t.test(power = .90, delta = 1, alternative = "one.sided")
```

```{r}
power.prop.test(n = 50, p1 = .50, p2 = .75)      ## => power = 0.740
power.prop.test(p1 = .50, p2 = .75, power = .90) ## =>     n = 76.7
power.prop.test(n = 50, p1 = .5, power = .90)    ## =>    p2 = 0.8026
power.prop.test(n = 50, p1 = .5, p2 = 0.9, power = .90, sig.level=NULL)
                                                 ## => sig.l = 0.00131
power.prop.test(p1 = .5, p2 = 0.501, sig.level=.001, power=0.90)
                                                 ## => n = 10451937
try(
 power.prop.test(n=30, p1=0.90, p2=NULL, power=0.8)
) # a warning  (which may become an error)
## Reason:
power.prop.test(      p1=0.90, p2= 1.0, power=0.8) ##-> n = 73.37
```


```{r}
library(PASWR2)
```

Simulates 100 samples of size 30 from  a normal distribution with mean 100
and a standard deviation of 10.  From the 100 simulated samples, 90% confidence
intervals for the Mean are constructed and depicted in the graph. 


```{r}
cisim(samples = 100, n = 30, parameter = 100, sigma = 10, conf.level = 0.90)
```

Simulates 100 sample of size 30 from a normal distribution with mean 100
and a standard deviation of 10.  From the 100 simulated samples, 95% confidence
intervals for the variance are constructed and depicted in the graph.

```{r}
cisim(100, 30, 100, 10, type = "Var")

```

Simulates 100 samples of size 50 from a binomial distribution where the 
population proportion of successes is 0.5.  From the 100 simulated samples,
92% confidence intervals for Pi are constructed and depicted in the graph.

```{r}
cisim(100, 50, 0.5, type = "Pi", conf.level = 0.92)

```


## Power


```{r}
library(pwr)

help(pwr)
```


```{r}
## Power at mu=105 for H0:mu=100 vs. H1:mu>100 (sigma=15) 20 obs. (alpha=0.05) 
sigma <- 15
c <- 100
mu <- 105
d <- (mu-c)/sigma
pwr.norm.test(d=d, n=20, sig.level=0.05, alternative="greater")
```


```{r}
## Sample size of the test for power=0.80
pwr.norm.test(d=d, power=0.8, sig.level=0.05, alternative="greater")
```

```{r}
## Power function of the same test
mu <- seq(75,125,l=100)
d <- (mu-c)/sigma
plot(d, pwr.norm.test(d=d, n=20, sig.level=0.05, alternative="greater")$power,
    type = "l", ylim = c(0,1))
abline(h = 0.05)
abline(h = 0.80)
```

```{r}
## Power function for the two-sided alternative
plot(d,pwr.norm.test(d=d, n=20, sig.level=0.05, alternative="two.sided")$power,
    type="l", ylim=c(0,1))
abline(h = 0.05)
abline(h = 0.80)
```

```{r}
library(pwr2)
```

```{r}
## Example 3
n <- seq(2, 30, by=4)
f <- seq(0.1, 1.0, length.out=10)
pwr.plot(n=n, k=5, f=f, alpha=0.05)
```

```{r}
## Example 1
pwr.2way(a=3, b=3, alpha=0.05, size.A=4, size.B=5, f.A=0.8, f.B=0.4)
pwr.2way(a=3, b=3, alpha=0.05, size.A=4, size.B=5, delta.A=4, delta.B=2, sigma.A=2, sigma.B=2)

## Example 2
ss.2way(a=3, b=3, alpha=0.05, beta=0.1, delta.A=1, delta.B=2, sigma.A=2, sigma.B=2, B=100)
```

```{r}
library(pwrAB)
```

```{r}
# Search for power given other parameters
AB_t2n(N = 3000, percent_B = .3, mean_diff = .15, sd_A = 1,
sd_B = 2, sig_level = .05, alternative = 'two_sided')

# Search for sample size required to satisfy other parameters
AB_t2n(percent_B = .3, mean_diff = .15, sd_A = 1,
sd_B = 2, sig_level = .05, power = .8, alternative = 'two_sided')
```

[pwrss](https://pwrss.shinyapps.io/index/)    

[A Practical Guide to Statistical Power and Sample Size Calculations in R](https://cran.r-project.org/web/packages/pwrss/vignettes/examples.html)
    
```{r}
library(pwrss)
```
    
```{r}
#| warnings = FALSE
suppressWarnings(power.t.test(ncp = 3, plot = TRUE,
                              df = 99, alpha = 0.05, alternative = "greater"))
```
    
```{r}
power.t.test(ncp = c(0.50, 1.00, 1.50, 2.00, 2.50), plot = FALSE,
             df = 99, alpha = 0.05, alternative = "not equal")
```
    
    


```{r}
pwrss.t.2means(mu1 = 30, mu2 = 28, sd1 = 12, sd2 = 8, kappa = 1, 
               n2 = 50, alpha = 0.05,
               alternative = "not equal")
```

```{r}
pwrss.t.2means(mu1 = 30, mu2 = 28, sd1 = 12, sd2 = 8, kappa = 1, 
               power = .80, alpha = 0.05,
               alternative = "not equal")
```

```{r}
# normal approximation
pwrss.z.prop(p = 0.45, p0 = 0.50,
             alpha = 0.05, power = 0.80,
             alternative = "less",
             arcsin.trans = FALSE)
```
  
```{r}
pwrss.z.2props(p1 = 0.45, p2 = 0.50,
               alpha = 0.05, power = 0.80,
               alternative = "not equal")
```
  
 
## Effect Size

```{r}
#| warnings = FALSE

# generate two vectors with numbers using the means and sd from above 
# means = 101, 99, sd = 10
num1 <- rnorm(1000000, mean = 101, sd = 10)
num2 <- rnorm(1000000, mean = 99, sd = 10)
# perform t-test
tt <- t.test(num1, num2, alternative = "two.sided")
# extract effect size
effectsize::effectsize(tt)
```

    
    
    
[FAQ: HOW ARE THE LIKELIHOOD RATIO, WALD, AND LAGRANGE MULTIPLIER (SCORE) TESTS DIFFERENT AND/OR SIMILAR?](https://stats.oarc.ucla.edu/other/mult-pkg/faq/general/faqhow-are-the-likelihood-ratio-wald-and-lagrange-multiplier-score-tests-different-andor-similar/)
    
```{r}
library(lmtest)

#fit full model
model_full <- lm(mpg ~ disp + carb + hp + cyl, data = mtcars)

#fit reduced model
model_reduced <- lm(mpg ~ disp + carb, data = mtcars)

#perform likelihood ratio test for differences in models
lrtest(model_full, model_reduced)
```

## Report

```{r}
library(tidyverse)
library(report)
```

```{r}
report(iris)
```

```{r}
iris %>%
  select(-starts_with("Sepal")) %>%
  group_by(Species) %>%
  report() %>%
  summary()
```

```{r}
report(t.test(mtcars$mpg ~ mtcars$am))
```

```{r}
aov(Sepal.Length ~ Species, data = iris) %>%
  report()

```

```{r}
model <- glm(vs ~ mpg * drat, data = mtcars, family = "binomial")

report(model)
```

```{r}
model <- lm(Sepal.Length ~ Species, data = iris)

report_model(model)
# linear model (estimated using OLS) to predict Sepal.Length with Species (formula: Sepal.Length ~ Species)

report_performance(model)
# The model explains a statistically significant and substantial proportion of
# variance (R2 = 0.62, F(2, 147) = 119.26, p < .001, adj. R2 = 0.61)

report_statistics(model)
# beta = 5.01, 95% CI [4.86, 5.15], t(147) = 68.76, p < .001; Std. beta = -1.01, 95% CI [-1.18, -0.84]
# beta = 0.93, 95% CI [0.73, 1.13], t(147) = 9.03, p < .001; Std. beta = 1.12, 95% CI [0.88, 1.37]
# beta = 1.58, 95% CI [1.38, 1.79], t(147) = 15.37, p < .001; Std. beta = 1.91, 95% CI [1.66, 2.16]
```