The restriction on randomization mentioned in the splitplot designs can be extended to more than one factor. For the case where the restriction is on two factors the resulting design is called a splitsplitplot design. These designs usually have three different sizes or types of experimental units.
Example 14.4 of the text book (Design and Analysis of Experiments, Douglas C. Montgomery, 7th and 8th Edition) discusses an experiment in which a researcher is interested in studying the effect of technicians, dosage strength and wall thickness of the capsule on absorption time of a particular type of antibiotic. There are three technicians, three dosage strengths and four capsule wall thicknesses resulting in 36 observations per replicate and the experimenter wants to perform four replicates on different days. To do so, first, technicians are randomly assigned to units of antibiotics which are the whole plots. Next, the three dosage strengths are randomly assigned to splitplots. Finally, for each dosage strength, the capsules are created with different wall thicknesses, which is the splitsplit factor and then tested in random order.
First notice the restrictions that exist on randomization. Here, we can not simply randomize the 36 runs in a single block (or replicate) because we have our first hard to change factor, named Technician. Furthermore, even after selecting a level for this hard to change factor (say technician 2) we can not randomize the 12 runs under this technician because we have another hard to change factor, named dosage strength. After we select a random level for this second factor, say dosage strength of level 3, we can then randomize the four runs under these two combinations of two factors and randomly run the experiments for different wall thicknesses as our third factor.
The linear statistical model for the splitsplitplot design would be:
\(y_{ijkh}=\mu+\tau_i+\beta_j+(\tau\beta)_{ij}+\gamma_k+(\tau\gamma)_{ik}+(\beta\gamma)_{jk}+(\tau\beta\gamma)_{ijk} +\delta_h\)
\(+(\tau\delta)_{ih}+(\beta\delta)_{jh}+(\tau\beta\delta)_{ijh}+(\gamma\delta)_{kh} +(\tau\gamma\delta)_{ikh}+(\beta\gamma\delta)_{jkh}+(\tau\beta\gamma\delta)_{ijkh}+\varepsilon_{ijkh} \left\{\begin{array}{c} i=1,2,\ldots,r \\ j=1,2,\ldots,a \\ k=1,2,\ldots,b \\ h=1,2,\ldots,c \end{array}\right. \)
Using the Expected Mean Square approach mentioned earlier for splitplot designs, we can proceed and analyze the splitsplitplot designs, as well. Based on Expected Mean Squares given in Table 14.25 to build test statistics (assuming the block factor to be random and the other factors to be fixed), , and are whole plot, splitplot and splitsplitplot errors, respectively. Minitab handles this model exactly in this way by GLM. (This was Table 14.22 in the 7th edition. The 8th edition has only the factors and EMS without the list of subscripts.)
Factor  rRi  aFj  bFk  cFh  lRl  Expected Mean Square  

Whole plot  \(\tau_i\)  1  a  b  c  1  \(\sigma^2 + abcd \sigma_{\tau}^2\) 
\(\beta_j\)  r  0  b  c  1  \(\sigma^2 + \sigma_{\tau \beta}^2 + \dfrac{rbc \sum \beta_{j}^2}{(a  1)}\)  
\((\tau \beta)_{ij}\)  1  0  b  c  1  \(\sigma^2 + bc \sigma_{\tau \beta}^2\)  
Subplot  \(\gamma_k\)  r  a  0  c  1  \(\sigma^{2}+a c \sigma_{\tau \gamma}^{2}+\dfrac{r a c \sum \gamma_{k}^{2}}{(b1)}\) 
\((\tau \gamma)_{ik}\)  1  a  0  c  1  \(\sigma^2 + ac\sigma_{\tau \gamma}^2\)  
\((\beta \gamma)_{jk}\)  r  0  0  c  1 
\(\sigma^2 + c\sigma_{\tau \beta \gamma}^2 + \dfrac{rc \sum \sum(\beta \gamma)_{jh}^2}{(a  1)(b  1)}\) 

\(\tau \beta \gamma)_{ijk}\)  1  0  0  c  1 
\(\sigma^2 + c\sigma_{\tau \beta \gamma}^2\) 

Subsubplot  \(\delta_h\)  r  a  b  0  1  \(\sigma^{2}+a b \sigma_{\tau\delta}^{2}+\dfrac{r a b \sum \delta_{k}^{2}}{(c1)}\) 
\((\tau \delta)_{ih}\)  1  a  b  0  1  \(\sigma^{2}+a b \sigma_{\tau s}^{2}\)  
\((\beta \delta)_{jh}\)  r  0  b  0  1  \(\sigma^{2}+b \sigma_{\tau \beta s}^{2}+\dfrac{r b \sum \sum(\beta \delta)_{j k}^{2}}{(a1)(c1)}\)  
\((\tau \beta\delta)_{ijh}\)  1  0  b  0  1  \(\sigma^{2}+b \sigma_{\tau \beta \delta}^{2}\)  
\((\gamma \delta)_{kh}\) 
r  a  0  0  1  \(\sigma^{2}+a \sigma_{\tau \gamma \beta s}^{2}+\dfrac{r a \sum \sum(\gamma\delta)_{kh}^{2}}{(b1)(c1)}\)  
\((\tau \gamma \delta)_{ikh}\)  1  a  0  0  1  \(\sigma^2 + a\sigma_{\tau \gamma \delta}^2\)  
\((\beta \gamma \delta)_{jkh}\) 
r  0  0  0  1 
\(\sigma^2 + \sigma_{\tau \beta \gamma \delta}^2 + \dfrac{r \sum \sum (\gamma \delta)_{ijk}^2}{(a  1)(b  1)(c  1)}\) 

\((\tau \beta \gamma \delta)_{ijkh}\)  1  0  0  0  1  \(\sigma^2 + \sigma_{\tau \beta \gamma \delta}^2\)  
\(\epsilon_{l(ijkh)}\)  1  1  1  1  1  \(\sigma^2\) (not estimable)  
Table 14.25 (Design and Analysis of Experiments, Douglas C. Montgomery, 8th Edition) 
However, we can use the traditional splitplot approach and extend it to the case of splitsplitplot designs as well. Keep in mind, as mentioned earlier, we should pool all the interaction terms with the block factor into the error term used to test for significance of the effects, in each section of the design, separately.