Econometria 2, Proyecto

5.1

Given the fitted response function

\[\begin{aligned} \hat{y} = 72.0 + 3.6x_1 - 2.5x_2 \end{aligned}\]

which is found to fit well in a localized region of \((x_1, x_2)\).

1. Plot contours of constant response, \(y\), in the \((x_1, x_2)\) plane.

library(plotly)
x1 <- seq(0, 10, by=0.1)
x2 <- seq(0, 10, by=0.1)
y <- 72 +3.6*x1-2.5*x2
z <- cbind(x1, x2, y)
p <- plot_ly(z = z) %>% add_surface()
p

2. Calculate and plot the path of the steepest ascend generated from this response function. Use \(x_1\) to determine the path, with \(\Delta x_1 = 1\).

library(qualityTools)
set <- facDesign(k=2)

implicit list embedding of S4 objects is deprecatedimplicit list embedding of S4 objects is deprecated

dset <- as.data.frame(set)
y <- 72 +3.6*dset$A-2.5*dset$B
response(set) <- y
set

sa <- steepAscent(factors = c("A", "B"), response="y", data=set)

Steepest Ascent for set 

sa_set <- as.data.frame(sa)
y_hat <- 72 +3.6*sa_set$A.real-2.5*sa_set$B.real
response(sa) <- y_hat
sa

plot(sa)

5.2

Suppose a fitted first-order model is given by:

\[\begin{aligned} y=25+4x_1+3x_2-2.5x_3 \end{aligned}\]

Suppose we wish to maximize the response. Calculate several points along the path of the steepest ascent

set <- facDesign(k=3)

implicit list embedding of S4 objects is deprecatedimplicit list embedding of S4 objects is deprecatedimplicit list embedding of S4 objects is deprecated

dset <- as.data.frame(set)
y <- 25+4*dset$A+3*dset$B-2.5*dset$C
response(set) <- y
set

sa <- steepAscent(factors = c("A", "B", "C"), response="y", data=set)

Steepest Ascent for set 

d_sa <- as.data.frame(sa)
y_hat <- 25+4*d_sa$A.real+3*d_sa$B.real-2.5*d_sa$C.real
response(sa) <- y_hat
sa

plot(sa)

5.3

Reconsider the model in Exercise 5.1. Construct a path of steepest ascent, using \(x_2\) to define the path, with \(\Delta x_2=1\)

set <- facDesign(k=2)

implicit list embedding of S4 objects is deprecatedimplicit list embedding of S4 objects is deprecated

dset <- as.data.frame(set)
y <- 72 +3.6*dset$A-2.5*dset$B
response(set) <- y
sa_1 <- steepAscent(factors = c("A", "B"), response="y", data=set)

Steepest Ascent for set 

sa_1_set <- as.data.frame(sa_1)
y_hat <- 72 +3.6*sa_1_set$A.real-2.5*sa_1_set$B.real
response(sa_1) <- y_hat
sa_2 <- steepAscent(factors = c("B", "A"), response="y", data=set)

Steepest Ascent for set 

sa_2_set <- as.data.frame(sa_2)
y_hat <- 72 +3.6*sa_2_set$A.real-2.5*sa_2_set$B.real
response(sa_2) <- y_hat
d_sa2 <- as.data.frame(sa_2)
plot(sa_1, type="b")
lines(d_sa2$Delta,d_sa2$y_hat, col="blue", type="b", lty=2)

legend(0, 76,legend = c("Delta X1", "Delta X2"), col=c("black", "blue"), lty=1:2)

2. How different are the two paths that you have constructed?

\(\Delta x_2\) has bigger values faster than \(\Delta x_1\)

5.4

Consider the injection molding experiment in Example 5.2

1. Use the coded design units to determine how far the points Base \(+4\Delta\) is from the design center.

set <- facDesign(k=3)

implicit list embedding of S4 objects is deprecatedimplicit list embedding of S4 objects is deprecatedimplicit list embedding of S4 objects is deprecated

dset <- as.data.frame(set)
y <- 25+4*dset$A+3*dset$B-2.5*dset$C
response(set) <- y
sa <- steepAscent(factors = c("A", "B", "C"), size=4, response="y", data=set)

Steepest Ascent for set 

d_sa <- as.data.frame(sa)
y_hat <- 25+4*d_sa$A.real+3*d_sa$B.real-2.5*d_sa$C.real
response(sa) <- y_hat
sa

sa_original <- steepAscent(factors = c("A", "B", "C"), response="y", data=set)

Steepest Ascent for set 

d_sa <- as.data.frame(sa_original)
y_hat <- 25+4*d_sa$A.real+3*d_sa$B.real-2.5*d_sa$C.real
response(sa_original) <- y_hat
sa_original

d_sa_original <- as.data.frame(sa_original)
plot(sa, type="b")
lines(d_sa_original$Delta,d_sa_original$y_hat, col="blue", type="b", lty=2)

legend(0, 76,legend = c("+4 Delta", "+1 Delta"), col=c("black", "blue"), lty=1:2)

2. Suppose that we had selected \(x_2\) to express the points on the path of steepest descent with \(\Delta x_2=1\). Recompute the points in Table E5.4. Use the point Base \(+4\Delta\), and compute the distance of this point from the design center.

set <- facDesign(k=3)

implicit list embedding of S4 objects is deprecatedimplicit list embedding of S4 objects is deprecatedimplicit list embedding of S4 objects is deprecated

dset <- as.data.frame(set)
y <- 25+4*dset$A+3*dset$B-2.5*dset$C
response(set) <- y
sa <- steepAscent(factors = c("B", "A", "C"), size=4, response="y", data=set)

Steepest Ascent for set 

d_sa <- as.data.frame(sa)
y_hat <- 25+4*d_sa$A.real+3*d_sa$B.real-2.5*d_sa$C.real
response(sa) <- y_hat
sa

d_sa <- as.data.frame(sa)
y_hat <- 25+4*d_sa$A.real+3*d_sa$B.real-2.5*d_sa$C.real
response(sa) <- y_hat
sa

sa_original <- steepAscent(factors = c("B", "A", "C"), response="y", data=set)

Steepest Ascent for set 

d_sa <- as.data.frame(sa_original)
y_hat <- 25+4*d_sa$A.real+3*d_sa$B.real-2.5*d_sa$C.real
response(sa_original) <- y_hat
sa_original

d_sa_original <- as.data.frame(sa_original)
plot(sa, type="b")
lines(d_sa_original$Delta,d_sa_original$y_hat, col="blue", type="b", lty=2)

legend(0, 76,legend = c("+4 Delta", "+1 Delta"), col=c("black", "blue"), lty=1:2)

3. How different are the two points you have computed?

The points have the same center in 25, but then the dispersion is really big being \(+4\Delta\) having a value of 225 units, meanwhile the default \(\Delta\) only reaches until 33 untis.

5.5

In a metallurgy experiment it is desired to test the effect of four factors and their interactions on the concentration (percent by weight) of a particular phosphorus compound in costing material. The variables are: A, percent phosphorus in the refinement; B, percent remelted material; C, fluxing time, and D, holding time. The four factors are varied in a 24 factorial experiment with two castings taken at each factor combination. The 32 castings were made in random order, and the data are shown in this table:

Treatment	Replication I	Replication II	Total
(1)	30.3	28.6	58.9
a	28.5	31.4	59.9
b	24.5	25.6	50.1
ab	25.9	27.2	53.1
c	24.8	23.4	48.2
ac	26.9	23.8	50.7
bc	24.8	27.8	52.6
abc	22.2	24.9	47.1
d	31.7	33.5	65.2
ad	24.6	26.2	50.8
bd	27.6	30.6	58.2
abd	26.3	27.8	54.1
cd	29.9	27.7	57.6
acd	26.8	24.2	51.0
bcd	26.4	24.9	51.3
abcd	26.9	29.3	56.2
Total	428.1	436.9	865.0

results <- c (58.9, 59.9, 50.1, 53.1, 48.2, 50.7, 52.6, 47.1, 65.2, 50.8, 58.2, 54.1, 57.6, 51, 51.3, 56.2, 865)
set <- facDesign(k=4)

implicit list embedding of S4 objects is deprecatedimplicit list embedding of S4 objects is deprecatedimplicit list embedding of S4 objects is deprecatedimplicit list embedding of S4 objects is deprecated

total <- results[unlist(runOrd(set))]
response(set) <- total
set

With the table we need to create the function using lm

df <- as.data.frame(data)
func <- lm(total ~ A + B + C + D, data=df)
func

Call:
lm(formula = total ~ A + B + C + D, data = df)

Coefficients:
(Intercept)            A            B            C            D  
    54.0625      -1.4375      -0.3750       0.3875       0.3500  

With this we get the final equation which is:

\[\begin{aligned} y = 54.0625 -1.4375x_1 - 0.3750x_2 + 0.3875x_3 + 0.3500x_4 \end{aligned}\]

2. Construct a table of the path of steepest ascent in the coded design variables.

sa <- steepAscent(factors = c("A", "B", "C","D"), response="total", data=data)

Steepest Ascent for data 

d_sa <- as.data.frame(sa)
y_hat <- 54.0625 - 1.4375*d_sa$A.coded - 0.3750*d_sa$B.coded + 0.3875*d_sa$C.coded + 0.35*d_sa$D.coded
response(sa) <- y_hat
sa

5.6

Consider the \(2^2\) factorial design featured in the illustration in Figure 3.1 in Chapter 3. Factor A is the concentration of a reactant, factor B is the feed rate, and the response is the viscosity of the output material. As one can tell by the computed effects and the analysis of variance shown, the main effects dominate in the analysis. Using the following table:

Combination	I	II	III	IV	Total
A low, B low	145	148	147	140	580
A high, B low	158	152	155	152	617
A low, B high	135	138	141	139	553
A high, B high	150	152	146	149	597

1. Fit a first-order model

set <- facDesign(k=2)

implicit list embedding of S4 objects is deprecatedimplicit list embedding of S4 objects is deprecated

total <- c(553, 617, 597, 580)
response(set) <- total
set

Now we create the function using a lineal regresion

lm_source <- as.data.frame(set)
lm_model <- lm(total~A+B, data=lm_source)
lm_model

Call:
lm(formula = total ~ A + B, data = lm_source)

Coefficients:
(Intercept)            A            B  
     586.75        20.25       -11.75  

So the function is:

\[\begin{aligned} y = 586.75 + 20.25x_1 - 11.75x_2 \end{aligned}\]

2. Compute and plot the path of steepest ascent.

sa <- steepAscent(factors = c("A", "B"), response="total", data=set)

Steepest Ascent for set 

sa_source <- as.data.frame(sa)
y_hat <- 586.75 + 20.25*sa_source$A.real - 11.75*sa_source$B.real
response(sa) <- y_hat
sa

c)Suppose that in the metric of the coded variables, the concentration cannot exceed 3.5 and the feed rate cannot go below -2.7. Show the practical path of steepest ascent—that is, the path that takes account of these constraints.

set <- facDesign(k=2)

implicit list embedding of S4 objects is deprecatedimplicit list embedding of S4 objects is deprecated

total <- c(553, 617, 597, 580)
response(set) <- total
lows(set) <- c(-2.7, -2.7)
highs(set) <- c(3.5, 3.5)
set

sa <- steepAscent(factors = c("A", "B"), size=0.1, response="total", data=set)

Steepest Ascent for set 

sa_source <- as.data.frame(sa)
y_hat <- 586.75 + 20.25*sa_source$A.real - 11.75*sa_source$B.real
response(sa) <- y_hat
sa

5.7

It is stated in the text that the development of the path of steepest ascent makes use of the assumption that the model is truly first-order in nature. However, even if there is a modest amount of curvature or interaction in the system, the use of steepest ascent can be extremely useful in determining a future experimental region. Suppose that in a system involving \(x_1\) and \(x_2\) the actual model is given by:

\[\begin{aligned} E(y) = 14 + 5x_1 - 10x_2 + 3x_1x_2 \end{aligned}\]

Assume that \(x_1\) and \(x_2\) are in coded form:

1. Show a plot of the path of steepest ascent (based on actual parameters) if the interaction is ignored.

set <- facDesign(k=2)

implicit list embedding of S4 objects is deprecatedimplicit list embedding of S4 objects is deprecated

set_source <- as.data.frame(set)
e_y <- 14 + 5*set_source$A - 10*set_source$B 
response(set) <- e_y
sa_linear <- steepAscent(factors = c("A", "B"), response="e_y", data=set)

Steepest Ascent for set 

sa_source <- as.data.frame(sa_linear)
y_hat <- 14 + 5*sa_source$A.coded - 10*sa_source$B.coded 
response(sa_linear) <- y_hat
plot(sa_linear, type="b")

2. Show a plot of the path of steepest ascent for the model with interaction. Note that this path is not linear

set <- facDesign(k=2)

implicit list embedding of S4 objects is deprecatedimplicit list embedding of S4 objects is deprecated

set_source <- as.data.frame(set)
e_y <- 14 + 5*set_source$A - 10*set_source$B + 3*set_source$A*set_source$B
response(set) <- e_y
sa_interaction <- steepAscent(factors = c("A", "B"), response="e_y", data=set)

Steepest Ascent for set 

sa_source <- as.data.frame(sa_interaction)
y_hat <- 14 + 5*sa_source$A.coded - 10*sa_source$B.coded + 3*sa_source$A.coded*sa_source$B.coded
response(sa_interaction) <- y_hat
plot(sa_interaction, type="b")

3. Comment on the difference in the two paths.

This Graph explains better the difference

sa_interaction_source <- as.data.frame(sa_interaction)
plot(sa_linear, type="b")
lines(sa_interaction_source$Delta,sa_interaction_source$y_hat, col="blue", type="b", lty=2)

legend(0, 35,legend = c("No Interaction", "Interaction"), col=c("black", "blue"), lty=1:2)