Organizing Code with Structured Programming
1) Use the code that you worked on in Homework #8 (creating fake data sets), and re-organize it following the principles of structured programming. Do all the work in a single chunk in your R markdown file, just as if you were writing a single R script. Start with all of your annotated functions, preliminary calls, and global variables. The program body should be only a few lines of code that call the appropriate functions and run them in the correct order. Make sure that the output from one function serves as the input to the next. You can either daisy-chain the functions or write separate lines of code to hold elements in temporary variables and pass them along.
For more information on purpose of Code - Homework 8
# Modeled Chl A response to Canopy Disturbances
# 23 March 2022
# LSC
# Initialize ---------------------------------
library(ggplot2)
##################################################
# FUNCTION: make_df
# purpose: read in a .csv file
# input: x and y vectors of numeric. Must be same length
# output: data frame
#-------------------------------------------------
get_data <- function(beta = 0.2, n=20, shape = 1.235, scale =0.55, mean=0.2, sd=0.1){
beta0 <- beta
biomass <- rgamma(n = n, shape = shape, scale = scale)
interc <- rnorm(n = n , mean = 0.2, sd = 0.1)
gap_frac <- biomass*beta0 + interc
d_frame <- data.frame(biomass, gap_frac)
return(d_frame)
}
new_df <- get_data()
print(new_df)## biomass gap_frac
## 1 0.12442684 0.3631867
## 2 0.40669832 0.1079599
## 3 0.33937379 0.3711872
## 4 0.23248396 0.3386595
## 5 0.69400541 0.4019568
## 6 0.28666097 0.2865990
## 7 0.26222509 0.4186471
## 8 0.93170428 0.3465431
## 9 0.11924829 0.1184332
## 10 0.13353699 0.1799112
## 11 0.06849154 0.2827753
## 12 1.15303004 0.4763122
## 13 0.53339434 0.3374594
## 14 0.86938635 0.4280215
## 15 0.59471493 0.2637286
## 16 0.41301595 0.1278743
## 17 1.44679014 0.3915449
## 18 1.56632751 0.1958016
## 19 1.71834607 0.5513315
## 20 1.50614393 0.5138105##################################################
# FUNCTION: calculate_model;
# purpose: Fits an ordinary least squares regression model
# input: x and y vectors of numeric. Must be same length
# output: entire model summary from lm
#-------------------------------------------------
calculate_model <- function(d_frame=new_df) {
reg_model <-lm(gap_frac~biomass, data = d_frame)
return(reg_model)
}
new_model <- calculate_model()
print(new_model)##
## Call:
## lm(formula = gap_frac ~ biomass, data = d_frame)
##
## Coefficients:
## (Intercept) biomass
## 0.2416 0.1246##################################################
# FUNCTION: summarize_output
# purpose: pull elements from model summary list
# input: list from summary call of lm
# output: vector of regression residuals
#-------------------------------------------------
summarize_output <- function(reg_model = new_model) {
z <- summary(reg_model)
return(z$residuals)
}
summarize_output()## 1 2 3 4 5 6
## 0.106090395 -0.184313844 0.087303625 0.068096804 0.073878038 0.009284682
## 7 8 9 10 11 12
## 0.144378018 -0.011158409 -0.138017678 -0.078320413 0.032649868 0.091028524
## 13 14 15 16 17 18
## 0.029396425 0.078086298 -0.051976295 -0.165186747 -0.030348025 -0.240988465
## 19 20
## 0.095596490 0.084520711##################################################
# FUNCTION: graph_results
# purpose: graph data and fitted OLS line
# input: x and y vectors of numeric. Must be same length
# output: creates graph
#-------------------------------------------------
graph_results <- function(d_frame = new_df) {
reg_plot <- ggplot(data=d_frame) +
aes(x=biomass,y=gap_frac) +
geom_point() +
stat_smooth(method=lm,se=0.95) # default se=0.95
print(reg_plot)
}
graph_results()## `geom_smooth()` using formula 'y ~ x'
2) Once your code is up and working, modify your program to do something else: record a new summary variable, code a new statistical analysis, or create a different set of random variables or output graph. Do not rewrite any of your existing functions. Instead, copy them, rename them, and then modify them to do new things. Once your new functions are written, add some more lines of program code, calling a mixture of your previous functions and your new functions to get the job done.
### Changing biomass from gamma to poison
##################################################
# FUNCTION: get_data_pois
# purpose: read in a .csv file
# input: x and y vectors of numeric. Must be same length
# output: data frame
#-------------------------------------------------
get_data_pois <- function(beta = 0.2, n=20, lambda= 0.2, mean=0.2, sd=0.1){
beta0 <- beta
biomass <- rpois(n = n, lambda = lambda)
interc <- rnorm(n = n , mean = 0.2, sd = 0.1)
gap_frac <- biomass*beta0 + interc
d_frame <- data.frame(biomass, gap_frac)
return(d_frame)
}
new_df <- get_data_pois()
print(new_df)## biomass gap_frac
## 1 0 0.23857903
## 2 1 0.54943844
## 3 0 0.07575111
## 4 0 0.23748987
## 5 0 0.17332533
## 6 1 0.32548689
## 7 1 0.23230815
## 8 0 0.31443315
## 9 0 0.23163424
## 10 0 0.10721016
## 11 0 0.45593412
## 12 0 0.27074302
## 13 0 0.26781020
## 14 0 0.15223631
## 15 1 0.53275318
## 16 0 0.12681595
## 17 0 0.27569565
## 18 0 0.24510644
## 19 0 0.23975776
## 20 0 0.11116169##################################################
# FUNCTION: calculate_model_pois;
# purpose: Fits an ordinary least squares regression model
# input: x and y vectors of numeric. Must be same length
# output: entire model summary from lm
#-------------------------------------------------
calculate_model_pois <- function(d_frame=new_df) {
reg_model <-lm(gap_frac~biomass, data = d_frame)
return(reg_model)
}
new_model <- calculate_model_pois()
print(new_model)##
## Call:
## lm(formula = gap_frac ~ biomass, data = d_frame)
##
## Coefficients:
## (Intercept) biomass
## 0.2202 0.1898##################################################
# FUNCTION: summarize_output_pois
# purpose: pull elements from model summary list
# input: list from summary call of lm
# output: vector of regression residuals
#-------------------------------------------------
summarize_output_pois <- function(reg_model = new_model) {
z <- summary(reg_model)
return(z$residuals)
}
summarize_output_pois()## 1 2 3 4 5 6
## 0.01834878 0.13944178 -0.14447914 0.01725962 -0.04690492 -0.08450977
## 7 8 9 10 11 12
## -0.17768852 0.09420290 0.01140399 -0.11302009 0.23570387 0.05051277
## 13 14 15 16 17 18
## 0.04757995 -0.06799394 0.12275651 -0.09341431 0.05546540 0.02487619
## 19 20
## 0.01952751 -0.10906856##################################################
# FUNCTION: graph_results_pois
# purpose: graph data and fitted OLS line
# input: x and y vectors of numeric. Must be same length
# output: creates graph
#-------------------------------------------------
graph_results_pois <- function(d_frame = new_df) {
reg_plot <- ggplot(data=d_frame) +
aes(x=biomass,y=gap_frac) +
geom_point() +
stat_smooth(method=lm,se=0.95) # default se=0.95
print(reg_plot)
}
graph_results()## `geom_smooth()` using formula 'y ~ x'
#################################################
# change the sample size from 20 to 30
get_data(n=30, shape = 1.235, scale =0.55)## biomass gap_frac
## 1 0.342315845 0.4121720
## 2 1.951744634 0.5757330
## 3 0.508799643 0.1813206
## 4 0.094880404 0.1792969
## 5 0.009994031 0.2568756
## 6 1.110212899 0.2499884
## 7 0.120596242 0.3271959
## 8 0.875289996 0.4293181
## 9 0.839077712 0.4313533
## 10 0.914237449 0.4834356
## 11 0.742438944 0.2481419
## 12 2.819261424 0.7602160
## 13 0.198279745 0.1811358
## 14 0.397963519 0.3774891
## 15 0.317504493 0.4150632
## 16 0.562604165 0.4957129
## 17 0.640084950 0.2725481
## 18 0.136483620 0.2163718
## 19 1.251532698 0.3309921
## 20 0.209541541 0.1296989
## 21 1.734296945 0.5684372
## 22 0.259296613 0.4132105
## 23 0.718890126 0.4566350
## 24 0.245655597 0.3406776
## 25 0.614788870 0.3224699
## 26 0.231736602 0.3580435
## 27 1.041496474 0.3174002
## 28 0.563036793 0.2950244
## 29 0.693093530 0.2721701
## 30 0.806243521 0.2223082new_df_large <- get_data()
calculate_model(d_frame = new_df_large)##
## Call:
## lm(formula = gap_frac ~ biomass, data = d_frame)
##
## Coefficients:
## (Intercept) biomass
## 0.1662 0.2403new_model_large <- calculate_model()
summarize_output(reg_model = new_model_large)## 1 2 3 4 5 6
## 0.01834878 0.13944178 -0.14447914 0.01725962 -0.04690492 -0.08450977
## 7 8 9 10 11 12
## -0.17768852 0.09420290 0.01140399 -0.11302009 0.23570387 0.05051277
## 13 14 15 16 17 18
## 0.04757995 -0.06799394 0.12275651 -0.09341431 0.05546540 0.02487619
## 19 20
## 0.01952751 -0.10906856graph_results(d_frame = new_df_large)## `geom_smooth()` using formula 'y ~ x'