För undersöka ett samband mellan mätresultat och några försöksvariabler används en försöksplan (DoE) och nedan presenteras en dylik. Med något datorprogram beräknas ett matematiskt samband, en ‘modell’. Denna kan i grafisk form ofta anges som ett ‘fisknät’ i 3D och man söker maximum (eller minimum) som då visar på bästa inställning av försöksvariablerna.
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 + 1.3X_{1} + 2.7X_{2} -3.2X_{1}^{2} -2.5X_{2}^{2} + 1.9X_{1}X_{2} + \epsilon \end{align} \]


Skapa försöksplanen. CCD-kommandot skapar en försöksplan där ‘2’ anger antal X-variabler. ‘n0’ anger här att planen skall innehålla 0 punkter i x1 = x2 = 0 för ‘box’-planen men 3 i x1 = x2 = 0 för ‘stjärn’-planen (se diagram nedan). ‘randomize = FALSE’ anger att utskriften sker i ‘run order’. (I verkligheten skall mätningar göras i slumpmässig ordning.) Längst till höger anges ‘block’ som anger att mätning i block 1 kan med fördel göras och analyseras först. (Om inget intressant visas kanske resten av mätningarna inte utförs.): 


library(rsm)                                                            # Funktioner för RSM

minPlan <- ccd (2, n0 = c(0,3), alpha = "rotatable", randomize = FALSE) # 'Central Composite Design' CCD

x1 <- minPlan$x1                                                        # 'Lyfter ut' x1 o x2-variablerna
x2 <- minPlan$x2                                                        # ur försöksplanen.

minPlan                                                                 # Skriver ut försöksplanen.
##    run.order std.order  x1.as.is  x2.as.is Block
## 1          1         1 -1.000000 -1.000000     1
## 2          2         2  1.000000 -1.000000     1
## 3          3         3 -1.000000  1.000000     1
## 4          4         4  1.000000  1.000000     1
## 5          1         1 -1.414214  0.000000     2
## 6          2         2  1.414214  0.000000     2
## 7          3         3  0.000000 -1.414214     2
## 8          4         4  0.000000  1.414214     2
## 9          5         5  0.000000  0.000000     2
## 10         6         6  0.000000  0.000000     2
## 11         7         7  0.000000  0.000000     2
## 
## Data are stored in coded form using these coding formulas ...
## x1 ~ x1.as.is
## x2 ~ x2.as.is

Planen som graf. Följande kommando skapar en graf över försöksplanen. De röda punkterna utgör block 1 och de röda är ‘stjärnformationen’. Det finns tre punkter i mitten av planen: 


library(ggplot2)                           # Rutiner som behövs för att rita graferna.

Fplan <- data.frame(x1, x2)                # Används för att rita x1x2-diagrammet.

farg <- c(rep("red", each = 4),rep("black", each = 7))

diaFplan <- ggplot(Fplan, aes(x = x1, y = x2)) + geom_point(color = farg, size = 3) + 
            xlab("Variabel x1") + ylab("Variabel x2") + scale_x_continuous(breaks=c(-1, 1),
            limits=c(-1.7, 1.7), labels=c("-1", "+1")) + scale_y_continuous(breaks=c(-1, 1),
            limits=c(-1.7, 1.7), labels=c("-1", "+1")) +
            theme(axis.text = element_text(color = "black", size = 13)) +
            theme(axis.title.x = element_text(size = 11, color = "blue", face = "bold"),
            axis.title.y = element_text(size = 11, color = "blue", face = "bold")) + coord_fixed()
diaFplan


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’).
En synnerligen intressant uppgift i tabellen är ‘Stationary point of response surface’ som anger optimum på responsytan (här ett maximum) uttryckt i ‘x1’- och ‘x2’-variablerna. (Detta kan också ses på nivådiagrammet nedan): : 


antrader <- length(x1)
epsilon <- rnorm(antrader, 0, 1)                                         # Slumpterm.
y <- 51.5 + 1.3*x1 + 2.7*x2 - 3.2*x1^2 - 2.5*x2^2 + 1.9*x1*x2 + epsilon  # Skapar mätvärdena för analys.

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

library(rsm)                                                             # Anropar funktioner för 'rsm'.
analys <- rsm(y ~ SO(x1, x2), data = minadata)                           # SO = 'Standard Order'.
summary(analys)                                                          # En sammanfattande tabell.
## 
## Call:
## rsm(formula = y ~ SO(x1, x2), data = minadata)
## 
##             Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 51.82594    0.50173 103.2945 1.612e-09 ***
## x1           1.64144    0.30725   5.3424 0.0030832 ** 
## x2           2.31676    0.30725   7.5404 0.0006498 ***
## x1:x2        1.30744    0.43451   3.0090 0.0297903 *  
## x1^2        -3.88879    0.36570 -10.6340 0.0001272 ***
## x2^2        -2.43524    0.36570  -6.6592 0.0011524 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Multiple R-squared:  0.9779, Adjusted R-squared:  0.9558 
## F-statistic: 44.24 on 5 and 5 DF,  p-value: 0.0003854
## 
## Analysis of Variance Table
## 
## Response: y
##             Df Sum Sq Mean Sq F value    Pr(>F)
## FO(x1, x2)   2 64.493  32.247 42.6997 0.0007195
## TWI(x1, x2)  1  6.838   6.838  9.0540 0.0297903
## PQ(x1, x2)   2 95.710  47.855 63.3671 0.0002807
## Residuals    5  3.776   0.755                  
## Lack of fit  3  1.618   0.539  0.4999 0.7194866
## Pure error   2  2.158   1.079                  
## 
## Stationary point of response surface:
##        x1        x2 
## 0.3047618 0.5574840 
## 
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] -2.184495 -4.139538
## 
## $vectors
##          [,1]       [,2]
## x1 -0.3581290 -0.9336721
## x2 -0.9336721  0.3581290

Ytterligare kommentar om utskriften. Om en term är ‘icke-signifikant’ bör analysen göras om med en enklare modell. Man börjar med att eliminera icke-signifikanta andragrads- och samspelstermer (x1:x2). Om en ‘x’-term är icke-signifikant medan motsvarande kvadratterm är signifikant, skall ‘x’-termen behållas i modellen.



Nivådiagram. Diagrammet avbildar responsytan som en karta. Det framgår att högsta nivån är ungefär 50 för de värden på ‘x1’- och ‘x2’-variablerna som angetts under ‘Stationary…’ ovan: 

contour(analys, ~ x1 + x2)                        # Skapar ett nivådiagram.

Responsyta. Diagrammet avbildar responsytan som ett 3D-diagram. Generellt gäller att ett 3D-diagram är svårare att tolka: 


persp(analys, x1 ~ x2, zlab="Y", contours=list(z = "bottom"))     # Skapar 3D-diagram inkl nivådiagram.

Avslutningsvis
Ovanstående text avser inte att ge en uttömmande diskussion om analys av data från en CCD-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.)