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:
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):
Sedan är det lätt att med lite matematik få en rät linje, här en linjen med lutningen -0.807:
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.)