So, we’re in the process of replacing our kitchen with an IKEA Euro-trash kitchen. This means that I have a lot of IKEA to build. It occurred to me as I built SEKTION cabinet after cabinet, each faster than the last, that it might be interesting to log my progression in IKEA building: an IKEA learning curve in the most literal sense.

It turns out that the kitchen we created was fraught with drawers. I’m not sure why we preferred drawers to doors so much, but I wasn’t about to let the opportunity go to waste. So, armed with a stack MAXIMERA drawers of various sizes, I sat down for an IKEA building and data collecting (and let’s face it, beer drinking) extravaganza.

My objectives were two-fold. First, I wanted to see just how much improvement I would see as I built drawer after drawer, memorizing the instructions and mastering the nuanced steps of hardware searching and proper pre-drilled hole selection. Second, if the time for each drawer assembly indeed decreased as I expected, I wanted to estimate ~~the theoretical human~~ my theoretical limit of IKEA drawer assembly performance. Noble cause.

I’m going to warn you now that there will come a point in this post that I will suggest you open Microsoft Excel and plot your data. I know. Be calm. In the words of Old Luke Skywalker: “Just breathe”. I’ll give you a minute here in case you need to find something to bite down on when this occurs.

*Non-linear model fitting in R has flummoxed me before, so if you have suggestions on how to do this better, faster, or more reliably, I would love to hear about it.*

### Methods

~~I had 21 drawers to build.~~ After our basement flooded mid-renovation (seriously), I had 19 drawers to build and 2 drawers to claim on insurance. I selected drawers in a random order (actually, the order that they went into the cabinets, which is pretty random) and timed drawer construction on a (borderline antique) iPhone 4S. Timing started once all drawer parts were taken out of the box and ended when the drawer was mounted on runners in the cabinet. Drawer leveling and centering were not included in the timing as the time required for this process is highly stochastic (sometimes drawers line up really well on the first shot and sometimes you need to fiddle forever to get them straight).

The basic drawer building steps were:

- Assemble drawer
- Attach hardware to drawer front
- Attach drawer front to drawer
- Determine mounting location in cabinet for runners
- Mount runners in cabinet
- Mount drawer on runners

I didn’t anticipate that the building process is slightly different for the two types of MAXIMERA drawers: low and medium/high (hereafter “high”). The high drawers have a few extra parts and an extra step or two to mount the side rails. Accordingly, I recorded the drawer type in the data and preserved it in the analysis as a categorical variable.

Times were plotted and curves were fitted in R…with help from Excel…damn.

### Data

Here’s what the final data look like. Note that, while I recorded the time (in mm:ss), I also include a column for decimal minutes, which are way easier to work with in R than actual times (Excel plays mean tricks on time variables too…because Excel sucks…even though I use it here…damn).

Order Type Time Min Sec DecMin 1 Low 20:17 20 17 20.28 2 High 21:29 21 29 21.48 3 High 23:54 23 54 23.90 4 High 15:32 15 32 15.53 5 High 20:48 20 48 20.80 6 Low 8:33 8 33 8.55 7 Low 7:30 7 30 7.50 8 High 12:51 12 51 12.85 9 Low 7:24 7 24 7.40 10 High 11:14 11 14 11.23 11 High 9:05 9 5 9.08 12 Low 7:18 7 18 7.30 13 High 10:48 10 48 10.80 14 Low 5:45 5 45 5.75 15 Low 4:25 4 25 4.42 16 High 8:40 8 40 8.67 17 High 8:32 8 32 8.53 18 Low 5:15 5 15 5.25 19 Low 4:35 4 35 4.58

### Analysis & Results

At first glance, you might be forgiven (well, probably not by a decent reviewer) for thinking that all the data together actually appear to fit a fairly reasonable (albeit a bit noisy) curve:

Of course, I know this not to be the case: the high drawers have extra building steps that must be accounted for. Adding the drawer type to the visualisation obviously demonstrates this difference:

Okay great, but what about those fitted curves? Well, they’re automatically generated by the stat_smooth command in ggplot2. I’ve guessed that they’re a second-order polynomial, so have given provided the command with the formula *y~poly(x,2*). But herein lies the first problem: **that curve formula doesn’t make any conceptual sense**.

Non-linear model fitting in R can be a bit tricky. I find it to be largely a bit of playing around and visualising various data and curves. There are certainly heaps of resources (publications, tutorials, etc.) for non-linear model fitting, using sophisticated statistical tools. However, for a simple analysis like this, I can use a simpler and quicker approach.

**Start by picking the right conceptual formula.** The model you choose has to make conceptual sense. In this case, I know that there will be some kind of physical limit on just how fast anyone can build and IKEA drawer. Even a robot specifically designed for MAXIMERA building will encounter limitations on just how fast pieces can be picked up and put together.

*Related Tangent: the same kind of thinking has generated some interesting philosophy around solving a Rubik’s Cube. While any cube can be solved in 26 quarter turns or less (the so-called “God’s Algorithm”), and while a suitably powerful computer may be able to solve this algorithm virtually instantly, it still takes time to turn the bloody cube. Hence, the world record speed (even for a computer and cube-turning robot) will always be >0 seconds. How much >0 is a matter of fascinating debate. There are similarly fascinating discussions in sport, of course.*

The same holds true for IKEA drawer building (though to my knowledge there are no MAXIMERA Tournaments). Knowing that our curve should always asymptote at some time >0 substantially narrows the options of curve formulae.

The most obvious choice here, given this criterion, is an exponential decay model. Such a model decreases exponentially and eventually asymptotes. The basic model formula is *y~a*e^(x*b)+c* (or *y~a*exp(x*b)+c* in R), where *y* is the response (time), *x* is the predictor (build order, a proxy for “experience”), *a* is the initial quantity of *y*, *b* is the decay constant (controls the steepness of the curve), and *c* is the asymptote (the value of *y* as *x* reaches infinity).

### Fitting the models using *nls* in R

To fit the models, I used the *nls* command from the base stats package in R. This command requires you to provide a formula for the model (see above) as well as “starting values” for all the coefficients (*a*, *b*, and *c*). The command uses these starting values to attempt to solve the curve for the data and, put simply, if the curve isn’t solved, it makes adjustments to these values and tries again. The command will, by default, only try a given number of times to solve the model before it gives up and also has a default threshold for what constitutes “solved”. You can control all these options by defining an *nls.control* object before running the *nls* command.

Some *nls* tips:

- If your model is fairly simple and your data not massive, crank up the number of iterations.
- Tell nls to only provide a warning when the model fit fails, rather than returning an error.
- If you get warnings related to the minimum factor, decrease it. This is usually because the model steps (changes to the starting values) are too large to accommodate fitting one of your coefficients.

# Make the nls.control object nls.set <- nls.control(maxiter=10000, minFactor=10e-10, warnOnly=T)

Now, what about starting values? Right. That’s where I (and most people I talk to) run into issues. If you’re a super mathey kind of person, you can maybe visualise a model formula from the data, but this kind of wizardry is beyond the skill-set of most ecologists/humans. And if you miss these starting values by any reasonable margin, the model will not converge. So you gotta get close. Enter Microsoft Excel (ugh).

A sneaky way to get the starting values for a curve formula (if your formula is common and simple) is to put the data in Excel, insert a scatterplot, add a trendline, and add the equation of the trendline to the plot. This gives you reasonable coefficients to use as starting values in the *nls* command. There are also some nice online tools that do much the same thing, but let’s face it, we all have Excel and our data is probably already in a CSV or XLS anyway, which we can open easily. There’s also Open Office if you just gotta go open source.

Excel tells me that my curve equation (for all data) is *y=21.509*exp(-0.078*x)*. Note that Excel won’t give you an exponential formula with an asymptote other than 0. But hopefully it’s enough to get me started. Since I know that *c* is the asymptote, I’ll start with the smallest value of *y* that I recorded. Note: the more coefficients you want to solve, the (much much much 10^much) more difficult it is to get the model to converge.

# Try the Excel-generated starting values # Note that the response variable (y) is the decimal minutes (DecMin) # and the predictor (x) is the "experience" or build order (Order). exp.order.nls <- nls(DecMin ~ a*exp(b*Order)+c, data=dat, start=list(a=21.5, b=-0.08, c=4.42), control=nls.set)

Nailed it.

summary(exp.order.nls) Formula: DecMin ~ a * exp(b * Order) + c Parameters: Estimate Std. Error t value Pr(>|t|) a 22.10313 3.12363 7.076 2.62e-06 *** b -0.14495 0.06821 -2.125 0.0495 * c 4.27484 3.19161 1.339 0.1992 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.222 on 16 degrees of freedom Number of iterations to convergence: 7 Achieved convergence tolerance: 1.683e-06

According to *nls*, the model formula is *y=22.10*exp(-0.14)+4.27*. So, when the type of drawer (low vs. high) is ignored, the model plots like this:

The raw data are shown as points (coloured by drawer type, but remember that this model is fitted for all data). I’ve pushed the curve out past 50 drawers and added a 95% confidence interval (CI), which is a bit tricky to generate and plot, though there is a nice explanation here. The CI would certainly suggest that this model, well, sucks.

We can do much better if we separate out the low and high drawers. Again, I got my starting values from Excel. The first clue that we are doing better when separating the drawer types is the measurement of residual standard error (much lower than above):

# LOW exp.type.nls.low <- nls(DecMin ~ a*exp(b*Order)+c, data=dat[dat$Type=="Low",], start=list(a=20, b=-1, c=3), control=nls.set) summary(exp.type.nls.low) Formula: DecMin ~ a * exp(b * Order) + c Parameters: Estimate Std. Error t value Pr(>|t|) a 19.88818 1.40330 14.172 7.71e-06 *** b -0.27245 0.04324 -6.301 0.000745 *** c 5.06128 0.51257 9.874 6.22e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.8543 on 6 degrees of freedom Number of iterations to convergence: 10 Achieved convergence tolerance: 2.28e-06 # HIGH exp.type.nls.high <- nls(DecMin ~ a*exp(b*Order)+c, data=dat[dat$Type=="High",], start=list(a=20, b=-0.2, c=5), control=nls.set) summary(exp.type.nls.high) Formula: DecMin ~ a * exp(b * Order) + c Parameters: Estimate Std. Error t value Pr(>|t|) a 22.50977 3.38406 6.652 0.00029 *** b -0.13722 0.08928 -1.537 0.16819 c 5.98501 4.70045 1.273 0.24357 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.422 on 7 degrees of freedom Number of iterations to convergence: 8 Achieved convergence tolerance: 3.291e-06

After again generating the CIs for each model, we can see that there are two clearly different curves for each the low and high drawers:

We also know from the coefficients that the asymptotes for each differ, with the low drawers (the simpler to build) asymptoting expectedly lower (5.06 min) than the high drawers (5.99 min). We can demonstrate this also be predicting times for very large values of x (build order):

predict(exp.type.nls.low, data.frame(Order=10^(1:10))) [1] 6.365540 5.061277 5.061277 5.061277 5.061277 5.061277 5.061277 [8] 5.061277 5.061277 5.061277 predict(exp.type.nls.high, data.frame(Order=10^(1:10))) [1] 11.692249 5.985036 5.985011 5.985011 5.985011 5.985011 5.985011 [8] 5.985011 5.985011 5.985011

After ~100 drawers, we see essentially no predicted improvement in time.

### The limits of human MAXIMERA drawer building performance

But this asymptote doesn’t seem to make sense as a measure of a limit of human performance, as it is actually higher than a few of my own times. Right, that’s because the model is based on least-squares, meaning it minimises the error around all points. In effect, it doesn’t care that there are values below the fitted model; it just considers that noise.

If we want an idea of the lower limits (i.e. the theoretical limits of human MAXIMERA drawer building performance), we can look at either 1) the best time ever recorded (4.42 min), or 2) look at the CIs (specifically the lower CIs).

The CIs for low drawer building are easily visualised in the plot and calculated in R using the *predictNLS* command from the *propagate* package. This command will provide CIs for any predicted value from the *nls* model. Cool. Also very useful for plotting (see above). To get the lowest 95% CI for the models, we need only to solve the CIs for a very large prediction of x. Here’s the example for low drawers:

# Make table of new x values newvals <- c(1:20,seq(30,90,10),seq(100,1000,100)) # Get fitted values and CIs exp.type.nls.low.pred <- data.frame(Order=newvals, Type="Low", DecMin=predict(exp.type.nls.low, data.frame(Order=newvals), type="response")) exp.type.nls.low.prop <- predictNLS(exp.type.nls.low, newdata=exp.type.nls.low.pred) exp.type.nls.low.pred$avg <- exp.type.nls.low.prop$summary[,"Prop.Mean.1"] exp.type.nls.low.pred$lci <- exp.type.nls.low.prop$summary[,"Prop.2.5%"] exp.type.nls.low.pred$uci <- exp.type.nls.low.prop$summary[,"Prop.97.5%"] # Look at upper and lower 95% CIs for predictions exp.type.nls.low.pred Order Type DecMin avg lci uci 1 Low 20.206372 20.206372 18.109764 22.289258 2 Low 16.594455 16.594455 14.979048 18.232091 3 Low 13.843936 13.843936 12.240034 15.522490 4 Low 11.749381 11.749381 10.151131 13.473448 5 Low 10.154352 10.154352 8.646588 11.829490 ... 30 Low 5.066886 5.066886 3.830439 6.312304 40 Low 5.061644 5.061644 3.809142 6.315207 50 Low 5.061301 5.061301 3.807247 6.315464 60 Low 5.061278 5.061278 3.807086 6.315480 70 Low 5.061277 5.061277 3.807073 6.315481 80 Low 5.061277 5.061277 3.807072 6.315481 90 Low 5.061277 5.061277 3.807072 6.315481 100 Low 5.061277 5.061277 3.807072 6.315481

So, according to our fitted model, there is a 95% probability that the true value of the low drawer asymptote falls above 3.81 minutes. Stated inversely, there is only a 5% chance (given the data collected) that the absolute minimum low drawer building time is below 3 min 49 sec. I’ve clearly got some training to do.

The theoretical limit for high drawers is trickier and can’t really be solved from the data collected. Using the same approach as for low drawers above, the lower CI for the asymptote actually falls below zero (-5.13 or -5 min 8 sec). Does that mean that our model is wrong? Of course it is. But it’s not wrong from a least-squares approach.

**A good lesson here:** while your model may be of the best conceptual structure and may be the optimal fit for that structure, it may still be wrong due to the limitations of your data.

I have no doubt (especially give the beautiful success of the low drawer model) that, with a bit more data collection, a better model fit could be achieved and we could (finally!) settle the contentious debate about the limits of IKEA MAXIMERA high drawer building performance.

Thanks David, I might actually use this in teaching at some point – would that be ok?

LikeLike

Of course! Happy that you like it. If you have improvement suggestions (especially for finding starting values), I would love to hear them. You are the guru!

LikeLike