# Methods 02 - Portfolio Assignment 1
## Study group 7
#### Amélie Kretschmer, Aya Sofia Nikolajsen, Carl Emil Grum-Nymann, Laura Sørine Voldgaard and Sofie Bøjgaard Thomsen

knitr::opts_chunk$set(message=FALSE, error=FALSE, warning=FALSE, comment=NA)

## 1. Continuous probability simulation
#### Exercise 5.2 from ROS

* log weights for men= M=5.13, sd=0.17
* log weights for women= M=4.96, sd=0.20
* 10 adults selected at random step on an elevator with a capacity of 1750 pound

*What is the probability of the weight of 10 elevator passengers exceeding the weight limit?*

exceeding <- 0
all_totals <- c()

for(i in 1:1000){
female <- sum(rbinom(10,1,0.5))
male <- 10 - female

totfemale <- sum(rlnorm(female, meanlog = 4.96, sdlog = 0.2))
totmale <- sum(rlnorm(male, meanlog = 5.13, sdlog = 0.17))

totalweight <- totfemale + totmale

all_totals[i] <- totalweight

if(totalweight > 1750){
    exceeding <- exceeding + 1
}

}

print(exceeding/1000)

hist_break <- hist(all_totals)$breaks

color_list <- rep('darkolivegreen1', length(hist_break))
color_list[hist_break >= 1750] <- 'deeppink'

hist(all_totals, col=color_list, xlab="Weight of 10 passengers (in pounds)", ylab="Frequence of occuring weight", main="Distribution of weightcomposition of elevator passengers")
abline(v=1750, col="blue", lwd=3, lty=2)

For exercise 1, we created a for loop to simulate 1000 possible outcomes of the two genders and the total weight when 10 people are in the elevator. For each time the weight exceeded 1750 pounds, the simulation noted that - and therefore the (simulated) probability of the weight exceeding the limit has in one of our simulations been estimated to ... - Future simulations should reveal a similar result.

We created a histogram to visualize the distribution of simulations and the number of simulations that exceed the limit of 1750 pounds.

## 2. Propagation of uncertainty
#### Exercise 5.6 from ROS

* cost savings at \$5 per unit
with a standard error of \$4

* size of the market 40000
with a standard error of 10 000

*Estimate of the total amount of money saved by the new product:*

all_savings <- c()
for(i in 1:1000){
unit_savings <- rnorm(1, mean = 5, sd = 4)
market_size <- rnorm(1, mean = 40000, sd = 10000)

total_savings <- market_size * unit_savings
all_savings[i] <- total_savings
}

hist(all_savings, col="pink", xlab="Amount of money saved by new product (in $)", ylab="Frequency of saved amount", main="Distribution of total savings")
mean(all_savings)
sd(all_savings)

In order to estimate the total amount of money saved by the new product, we created a for loop to simulate 1000 possible scenarios of savings per unit times number of sold products. In one of our simulations we got a mean value of \$ of savings in total with a standard deviation of \$. Future simulations should reveal a similar result.

We again visualised the distribution of simulations in a histogram.

## 3. Inference for a ratio of parameters
#### Exercise 5.10 from ROS

* difference in costs between treatments A and B is estimated at \$600 per patient:
* with a standard error of \$400
* based on a regression with 50 degrees of freedom

* difference in effectiveness is estimated at 3.0 (on some relevant measure)
* with a standard error of 1.0
* based on a regression with 100 degrees of freedom

*a. Create 1000 simulation draws of the cost difference and the effectiveness difference, and make a scatterplot of these draws:*

sd_cost <- sqrt(51)*400
sd_effect <- sqrt(101)*1

cost_dif <- c()
effect_dif <- c()

for(i in 1:1000){
cost_difference <- rnorm(1, mean = 600, sd = sd_cost)
effect_difference <- rnorm(1, mean = 3.0, sd = sd_effect)

cost_dif[i] <- cost_difference
effect_dif[i] <- effect_difference

}

plot (cost_dif, effect_dif, col="orange", xlab="Cost difference", ylab="Effect difference", main="Effectiveness difference compared to cost difference")

We calculated the sd for cost and effectiveness difference between the treatments A and B using the equation:
>SD = sqrt(n+1)*SE

Then we created a simulation again, to simulate the cost and effectiveness difference and visualized this in a scatterplot.

*b. Use simulation to come up with an estimate, 50% interval, and 95% interval for the incremental cost-effectiveness ratio:*

effect_of_cost <- cost_dif/effect_dif

print(interval50 <- quantile(effect_of_cost, c(.25, .75)))
print(interval95 <- quantile(effect_of_cost, c(.025, .975)))
	

Then we calculated the necessary increase in cost to increase effectiveness of the treatment by one unit.
Using the quantile(), we found the estimate for a 50% interval, and a 95% interval for the incremental cost-effectiveness ratio.

*c. Repeat, changing the standard error on the difference in effectiveness to 2.0:*

sd_cost_c <- sqrt(51)*400
sd_effect_c <- sqrt(101)*2

cost_dif_c <- c()
effect_dif_c <- c()

for(i in 1:1000){
cost_difference_c <- rnorm(1, mean = 600, sd = sd_cost_c)
effect_difference_c <- rnorm(1, mean = 3.0, sd = sd_effect_c)

cost_dif_c[i] <- cost_difference_c
effect_dif_c[i] <- effect_difference_c

}

plot (cost_dif_c, effect_dif_c, col="orange", main="Effectiveness difference compared to cost difference")

effect_of_cost_c <- cost_dif_c/effect_dif_c

print(interval50_c <- quantile(effect_of_cost_c, c(.25, .75)))
print(interval95_c <- quantile(effect_of_cost_c, c(.025, .975)))

Here we repeated the process as described above just by changing the calculated variables.