In this lab, we will focus on linear and non-linear programming.
Linear programming, as discussed in the previous lab, works with simple and multiple linear regression techniques; sometimes the variables have completely direct or completely non-direct relationships and these techniques can model them.
Sometimes, however, the variables do not predict each other in a linear way. For example, looking at the stock market vs. time, we know that generally the market was booming before the crash, then the market crashed and the great depression hit, and slowly the market started to rise again.
This pattern is not linear, and in fact a non-linear programming technique can be used to model it and predict the value of the market based on the year.
In this lab, we will explore topics like optimization, solve a marketing model, and perform linear and non-linear regression on the cost of servers.
We are going to use tidyverse a collection of R packages designed for data science.
## Loading required package: lpSolveAPI
library("tidyverse")
## -- Attaching packages ---------------------------------------------------------------------------------------------------------------------- tidyverse 1.2.1 --
## v ggplot2 2.2.1 v purrr 0.2.4
## v tibble 1.4.2 v dplyr 0.7.4
## v tidyr 0.7.2 v stringr 1.2.0
## v readr 1.1.1 v forcats 0.2.0
## -- Conflicts ------------------------------------------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
lprec <- make.lp(0, 2)
lp.control(lprec, sense="max")
## $anti.degen
## [1] "fixedvars" "stalling"
##
## $basis.crash
## [1] "none"
##
## $bb.depthlimit
## [1] -50
##
## $bb.floorfirst
## [1] "automatic"
##
## $bb.rule
## [1] "pseudononint" "greedy" "dynamic" "rcostfixing"
##
## $break.at.first
## [1] FALSE
##
## $break.at.value
## [1] 1e+30
##
## $epsilon
## epsb epsd epsel epsint epsperturb epspivot
## 1e-10 1e-09 1e-12 1e-07 1e-05 2e-07
##
## $improve
## [1] "dualfeas" "thetagap"
##
## $infinite
## [1] 1e+30
##
## $maxpivot
## [1] 250
##
## $mip.gap
## absolute relative
## 1e-11 1e-11
##
## $negrange
## [1] -1e+06
##
## $obj.in.basis
## [1] TRUE
##
## $pivoting
## [1] "devex" "adaptive"
##
## $presolve
## [1] "none"
##
## $scalelimit
## [1] 5
##
## $scaling
## [1] "geometric" "equilibrate" "integers"
##
## $sense
## [1] "maximize"
##
## $simplextype
## [1] "dual" "primal"
##
## $timeout
## [1] 0
##
## $verbose
## [1] "neutral"
set.objfn(lprec, c(275.691, 48.341))
add.constraint(lprec, c(1, 1), "<=", 350000)
add.constraint(lprec, c(1, 0), ">=", 15000)
add.constraint(lprec, c(0, 1), ">=", 75000)
add.constraint(lprec, c(2, -1), "=", 0)
lprec
## Model name:
## C1 C2
## Maximize 275.691 48.341
## R1 1 1 <= 350000
## R2 1 0 >= 15000
## R3 0 1 >= 75000
## R4 2 -1 = 0
## Kind Std Std
## Type Real Real
## Upper Inf Inf
## Lower 0 0
# solve
solve(lprec)
## [1] 0
get.objective(lprec)
## [1] 43443517
get.variables(lprec)
## [1] 116666.7 233333.3
Name your dataset ‘mydata’ so it easy to work with.
Commands: read_csv() head()
mydata <- read.csv(file = "data/ServersCost.csv")
head(mydata)
## servers cost
## 1 1 27654
## 2 2 24789
## 3 3 21890
## 4 4 21633
## 5 5 15843
## 6 6 12567
servers <- mydata$servers
head(servers)
## [1] 1 2 3 4 5 6
cost <- mydata$cost
head(cost)
## [1] 27654 24789 21890 21633 15843 12567
cor(mydata)
## servers cost
## servers 1.00000000 0.03356606
## cost 0.03356606 1.00000000
The correlation between the two variables, cost and servers, is very weak since it is around .034. ### 2C) Create a plot for the dependent (y) and independent (x) variables. Note any patterns or relation between the two variables describe the trend line.
Commands: p <- qplot( x = INDEPENDENT, y = DEPENDENT, data = mydata) + geom_point()
p <- qplot( x = servers, y = cost, data = mydata) + geom_point()
p
The relationship between cost and servers decreases at a steep rate when the number of servers is 0 to 9. Then it levels out from 9 to around 13 and then increases at a steep rate. Commmand: p + geom_smooth(method = “lm”)
p + geom_smooth(method = "lm")
linear_model <- lm(cost ~ servers)
summary(linear_model)
##
## Call:
## lm(formula = cost ~ servers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10646.2 -8646.2 -544.7 7066.0 12858.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14747.2 4035.5 3.654 0.00181 **
## servers 48.0 336.9 0.142 0.88828
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8687 on 18 degrees of freedom
## Multiple R-squared: 0.001127, Adjusted R-squared: -0.05437
## F-statistic: 0.0203 on 1 and 18 DF, p-value: 0.8883
We use a transformation and use a nonlinear quadratic model to see how the model fits to the data.
Quadratic Model: y = x + x^2
sqservers <- servers^2
quad_model <- lm(cost ~ servers + sqservers)
quad_model
##
## Call:
## lm(formula = cost ~ servers + sqservers)
##
## Coefficients:
## (Intercept) servers sqservers
## 35417.8 -5589.4 268.4
summary(quad_model)
##
## Call:
## lm(formula = cost ~ servers + sqservers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2897.8 -1553.4 -513.2 1152.4 4752.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 35417.77 1742.64 20.32 2.30e-13 ***
## servers -5589.43 382.19 -14.62 4.62e-11 ***
## sqservers 268.45 17.68 15.19 2.55e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2342 on 17 degrees of freedom
## Multiple R-squared: 0.9314, Adjusted R-squared: 0.9233
## F-statistic: 115.4 on 2 and 17 DF, p-value: 1.282e-10
This is a great fit for the data because the r-squared is .9314 and the adjusted r-squared is .9233. These numbers are a lot higher than that of the linear model meaning this model is a much better fit.
Commands: predicted_2 <- predict( quad_model, data = mydata )
sqservers = servers^2
quad_model = lm(cost ~ servers + sqservers )
predicted2 = predict(quad_model,data=mydata)
predicted2
## 1 2 3 4 5 6 7
## 30096.790 25312.706 21065.520 17355.233 14181.844 11545.354 9445.762
## 8 9 10 11 12 13 14
## 7883.068 6857.273 6368.376 6416.377 7001.277 8123.076 9781.772
## 15 16 17 18 19 20
## 11977.367 14709.861 17979.252 21785.543 26128.731 31008.818
Commands: qplot( x = DEPENDENT, y = INDEPENDENT/PREDICTED, colour = “red” )
qplot( x = servers, y = predicted2, colour = "red" )
This graph looks very similar to the real data and it has the same negative slope to around 10 then increases afterwards.
cubeservers <- servers^3
cubic_model <- lm(cost ~ servers + sqservers + cubeservers)
cubic_model
##
## Call:
## lm(formula = cost ~ servers + sqservers + cubeservers)
##
## Coefficients:
## (Intercept) servers sqservers cubeservers
## 36133.696 -5954.738 310.895 -1.347
summary(cubic_model)
##
## Call:
## lm(formula = cost ~ servers + sqservers + cubeservers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2871.0 -1435.1 -473.6 1271.8 4600.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36133.696 2625.976 13.760 2.77e-10 ***
## servers -5954.738 1056.596 -5.636 3.72e-05 ***
## sqservers 310.895 115.431 2.693 0.016 *
## cubeservers -1.347 3.619 -0.372 0.715
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2404 on 16 degrees of freedom
## Multiple R-squared: 0.932, Adjusted R-squared: 0.9193
## F-statistic: 73.11 on 3 and 16 DF, p-value: 1.478e-09
This is a great model for the data as well because r-squared is .932 and adjusted r-squared is .9193 which are both very high. The r-square is higher than the second model but the adjusted r-squared is a bit lower.
Commands: predicted3 <- predict( cubic_model, data = mydata )
predicted3 <- predict( cubic_model, data = mydata )
Commands: qplot( x = DEPENDENT, y = INDEPENDENT/PREDICTED, colour = “red” )
qplot( x = servers, y = predicted3, colour = "red" )
This data looks very similar to the second model because they both had a decreasing slope to around 10 then increased from then on.
variables: LINEAR_MODEL , PREDICTED_QUADRATIC, PREDICTED_CUBIC
# Black = Actual Data
plot(servers, cost, pch = 16)
# Blue = Linear Line based on Linear Regression Model
abline(linear_model, col = "blue", lwd = 2)
# Red = Quadratic Model based on Quadratric Regression found above
# Needed to overlay new points without the labels and annotations
par(new = TRUE, xaxt = "n", yaxt = "n", ann = FALSE)
plot(predicted2, col = "red", pch = 16)
# Green = Cubic Model based on Cubic Regression found above
# Overlay new points without the labels and annotations
par(new = TRUE, xaxt = "n", yaxt = "n", ann = FALSE)
plot(predicted3, col = "green", pch = 16)
The green dots which are the cubic model seem to be slightly closer to the black dots which are the real data points than the red dots of the quadratic model. Therefore, the cubic model is the best fit for this data set. The linear model does not predict any of the data properly.