Vid ett första experiment kan man inte förvänta att genast hitta optimala inställningar på försöksvariablerna. ‘Steepest Ascent’ är en metod att utnytta den första modellen och se vart den pekar för ett ännu bättre resultat.
Man använder helt enkelt de från analysen skattade β-parametrarna i modellen (som här bara är ett plan dvs två βX-termer) och med lite matematik erhålles i vilken riktning försöksplanen bör flyttas.
Skärmdumpen nedan visar ett exempel från https://ovn.ing-stat.se/ [41b] och används som illustration i detta dokument.


Koder, grafer och övrig programmering kan utföras på flera olika sätt men har här hållits till ett minimum för tydlighetens skull.



 En responsyta från https://ovn.ing-stat.se/ [41b]

Sliderna till vänster kan användas för att illustrera hur ytan förändras med olika inställningar. De representerar koefficienterna i den linjära modellen (Z är mätresultatet, X och Y är försöksvariablerna). Om t.ex. alla sliderna (utom den översta) sätts till 0 bildas ett plan i en riktning (som egentligen kan ritas som en vanlig 2D-graf). Det finns ytterligare exempel och kommentarer med knapparna nederst till vänster.

I verkligheten känner man naturligtvis inte till koefficienterna, dessa behöver skattas från data i försöksplanen. Nedan adderas en simulerad slumpkomponent till modellen för en realistisk situation.
Varför en linjär modell bestående av flera andragradstermer? Ett komplicerat matematiskt uttryck kan motsvaras av ett litet antal termer som lätt kan analyseras. Krökta samband kan oftast modelleras med andragradstermer (och om sambandet skulle vara ännu mer komplicerat skulle det behövas långt större datamängder för att hitta kanske små irrlevanta förbättringar av en kanske redan tillfredsställande modell).



























Följande modell används nedan och innehåller samma koefficienter som skärmdumpen ovan. Sista termen ε är slumpkomponenten som innhåller summan av påverkan av alla andra förklaringsvariabler (kända och okända) som inte finns i modellen:


\[ \begin{align} \large Y = 51.5 + 2.1X_{1} + 2.6X_{2} + \epsilon \end{align} \]

Försöksvariablerna. Nedan anges försöksvariablernas faktiska värden samt deras kodning i analysen:


library(kableExtra) 

Variabel <- c("X1", "X2")
Låg   <- c(50, 120)
Hög   <- c(70, 150)
Kodning <- c( " -1, 1", " -1, 1")

data <- data.frame(Variabel, Låg, Hög, Kodning)               # Skapar en sk. 'dataframe'.

data <- kable(head(data), "html")
kable_styling(data, "striped", position = "left", font_size = 17, full_width = FALSE) %>%
  column_spec(1, width = "2cm") %>%
  column_spec(2, width = "3cm") %>%
  column_spec(3, width = "3cm") %>% 
  column_spec(4, width = "5cm")
Variabel Låg Hög Kodning
X1 50 70 -1, 1
X2 120 150 -1, 1

Kodning av variablerna sker med följande två enkla formler:


\[ \begin{align} \large \frac{50-60}{10}=-1 & \phantom{litet tomt uterymme} \large \frac{70-60}{10}=1\phantom{invisible}\text{(Låg/Hög nivå på variabel X1)} \\ \\ \large \frac{120-135}{15}=-1 & \phantom{litet tomt uterymme} \large \frac{150-135}{15}=1\phantom{invisib}\text{(Låg/Hög nivå på variabel X2)} \end{align} \]


Skapa försöksplanen. FrF2-kommandot skapar en försöksplan med nruns = 4 rader och med variabelnamn X1 och X2 i icke-slumpmässig ordning:


library(FrF2)                                                        # Funktioner för att skapa försöksplaner.
minPlan <- FrF2(nruns = 4, factor.names = c("X1", "X2"), randomize = FALSE) 
minPlan                                                              # Skriver ut försöksplanen.
##   X1 X2
## 1 -1 -1
## 2  1 -1
## 3 -1  1
## 4  1  1
## class=design, type= full factorial


Analys av data. Med hjälp av modellen ovan skapas mätvärdena (Y) som lagras tillsammans med X-variablerna i ‘minadata’-tabellen. Analysen av data skapar en mer eller mindre standardutskrift, likt den som finns vid s.k. regressionsanlys.
Varje rad ‘Intercept’, ‘x1’, osv är en term i den linjära modellen. Om P-värdet (‘Pr(>|t[)’) längst till höger är lågt, säg, < 0.05, anses termen platsa i den modell som bäst förklarar sambandet (termen är ‘signifikant’):


x1 <- as.numeric(as.character(minPlan$X1))
x2 <- as.numeric(as.character(minPlan$X2))

antrader <- length(x1)
epsilon <- rnorm(antrader, 0, 0.1)                                       # Slumpterm.
Y <- 51.5 + 2.1*x1 + 2.7*x2 + epsilon                                    # Skapar mätvärdena för analys.

minadata <- cbind(minPlan, Y)                                            # Lägger 'Y' i 'X'-tabellen.                                     

analys <- lm(Y ~ x1 + x2, data = minadata)                    
summary(analys)                                                          # En sammanfattande tabell.
## 
## Call:
## lm.default(formula = Y ~ x1 + x2, data = minadata)
## 
## Residuals:
##        1        2        3        4 
##  0.02727 -0.02727 -0.02727  0.02727 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 51.54054    0.02727 1889.81 0.000337 ***
## x1           2.11847    0.02727   77.68 0.008195 ** 
## x2           2.70927    0.02727   99.34 0.006408 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05455 on 1 degrees of freedom
## Multiple R-squared:  0.9999, Adjusted R-squared:  0.9998 
## F-statistic:  7951 on 2 and 1 DF,  p-value: 0.00793


Nivådiagram. Diagrammet nedan avbildar responsytan som en karta. Det framgår att högsta nivån är drygt 50:


xkontur <- seq(-1, 1, length.out = 100)
ykontur <- seq(-1, 1, length.out = 100)

mat_coeff <- summary(analys)$coefficients              # All estimated coefficients.

beta0 <- analys$coefficients[1]                        # (Intercept, beta0)
beta1 <- analys$coefficients[2]                        # (beta1)
beta2 <- analys$coefficients[3]                        # (beta2)

zkontur <- outer(xkontur, ykontur, function(xkontur, ykontur) { beta0 + beta1*xkontur + beta2*ykontur })
contour(xkontur, ykontur, zkontur, main = "Konturdiagram responsyta, X1 och X2", xlab = "X1-variabel", 
        ylab = "X2-variabel")

Kommentar till nivådiagrammet. Nivådiagrammet är en ‘höjdkarta’ (här ett lutande plan) och t.ex. en linje märkt ‘50’ binder ihop punkter med Y = 50. Det framgår att diagrammet lutar uppåt i både X1- och X2-variabeln (precis som dokumentets översta graf ovan).
Antag att man strävar efter Y-värden, säg, > 56, då skall man förflytta planen i en riktning vinkelrätt mot nivålinjerna (röd pil, beräkningar sker nedan), och ta Y-mätningar för de nya X1- och X2-värderna (och kanske senare ett behov av en större modell som kan plocka upp eventuell krökning i responsytan (dvs en modell med andragradstemer)):





Beräkning av riktning. Antag att analysen givit uttrycket till vänster nedan. En linje med Y = 50 skapas genom att sätta in 50 i uttrycket (till höger nedan):


\[ \begin{align} \large Y = 51.5 + 2.1X_{1} + 2.6X_{2} & \phantom{litet tomt ut}\large 50 = 51.5 + 2.1X_{1} + 2.6X_{2} \end{align} \]

Sedan är det lätt att med lite matematik få en rät linje, här en linjen med lutningen -0.807:


\[ \begin{align} \large X_{2} = -0.538 - 0.807X_{1} \end{align} \]

Denna beräkning har utförts tre gånger och dessa tre linjer har ritats in i diagrammet nedan. För att snabbast öka utfallet bör alltså nya mätningar tas i den riktning (uppåt höger) som den röda linjen visar.
Det är lätt att i litteraturen eller på nätet hitta hur man beräknar lutningen på en linje som är vinkelrät mot en given linje. Antag att de två linjernas lutningar anges med k-värden. Då gäller följande:

\[ \begin{align} \large k_{1}\cdot k_{2} = -1 & \phantom{litet tomt ut}\large -0.807\cdot k_{2} = -1\phantom{litet tomt} \phantom{invie}\text{som ger lutningen ->} \phantom{litet}\large k_{2} = 1.239 \end{align} \]

library(ggplot2)

X1    <- c(+1, +1, +1, +1,   -1, -1, -1, -1)
X2    <- c(-1, -1, +1, +1,   -1, -1, +1, +1)

kontur <- ggplot() + geom_point(aes(x = X1, y = X2), color = "red", size = 3) + 
                     xlab("Variabel X1") + ylab("Variabel X2") + scale_x_continuous(breaks=c(-1, 1), limits=c(-1.2, 1.2),
                                              labels=c("-1", "+1")) + theme(axis.text = element_text(color = "black", size = 16)) +
                                              scale_y_continuous(breaks=c(-1, 1), limits=c(-1.2, 1.2),
                                              labels=c("-1", "+1")) + theme(axis.text = element_text(color = "black", size = 16)) +
                                              coord_fixed()
kontur <- kontur + geom_abline(intercept = -0.538, slope = -0.807, linetype = 1, colour = "black") +
                   geom_abline(intercept = -0.154, slope = -0.807, linetype = 1, colour = "black") +
                   geom_abline(intercept =  0.231, slope = -0.807, linetype = 1, colour = "black") +
                   geom_abline(intercept =      0, slope =  1.239, linetype = 2, linewidth = 2, colour = "red")

kontur <- kontur +  theme(axis.title.x = element_text(size = 12, color = "blue", face = "bold"),
                            axis.title.y = element_text(size = 12, color = "blue", face = "bold"))
kontur


Tolkning av riktning. Om de övre gränserna på försöksvariablerna är vad som är tekniskt möljigt blir maximum utfall det som erhålles då X1 = 1 och X2 = 1 sätts in i ekvationen (alltså övre högra hörnet). Men kanske försöksledaren funderar på att flytta försöksplanen enligt lutningen på den röda linjen.
Alternativt väljer försöksledaren i stället att mäta utfallet på en eller flera punkter utefter den röda linjen för att se eventuell förändring och möjligtvis behov av att modellera eventuell krökning i responsytan. Man bör också vara medveten om att riktningen är beräknad med en viss statistik osäkerhet och bör naturligtvis inte dras ut alltför långt.

Med bara två försöksvariabler är det möjligt att praktiskt illustrera med diagram som ovan. Om man i stället har flera försöksvariabler som blir signifikanta i analysen, krävs det en matematisk behandling för att hitta ‘steepest ascent’. Det finns en rik litteratur i bokform eller på nätet (t.ex. Emperical Model-Building and Response Surfaces av George Box och Norman Draper).



Avslutningsvis
Ovanstående text avser inte att ge en uttömmande diskussion om analys av data från en FrF2-plan. Det finns många kommando som kan ge mer information och flera sätt att automatisera en analys. Även graferna kan göras mer utförliga med t.ex. textinformation, färger m.m. Men här har framställningen poängterat enkelhet.

(Analyserna har gjorts med datorprogrammet ‘R’ och med det grafiska gränssnittet ‘R-studio’ och bägge är gratis tillgängliga på nätet. Se https://www.indstat.se och knappen [Statistikprogram - R] för installation.)

(Se https://www.indstat.se för många andra simuleringsövningar.)