I find the terms "confusion" and "diffusion" to be slightly nebulous and can lead to over-simplifications.
For example, saying that "substitution" is responsible for "confusion" is not necessarily correct: "Substitution" is actually just a function application to the state; The implementation often utilizes a memoized function, but you can easily 1. calculate the non-linear function explicitly on each application, as well as 2. use a memoized linear function to provide only diffusion and not confusion. So saying Substitution provides Confusion is an over simplification.
Now, the reason why non-linear functions are often implemented with a memoized table is because non-linear functions can be complicated to compute. Generally speaking, the more complicated a function is, the longer it takes to evaluate it. For example, the AES S-Box utilizes calculation of the multiplicative modular inverse of an element in a finite field. This could be computed explicitly instead of using a table lookup; Doing so will result in a cipher that is significantly slower.
So the confusion does not stem from substitution per se, it results from applying complicated non-linear functions that tend to produce maximally unhelpful and complicated equations. The more "linear" a non-linear function is, the easier it becomes to cryptanalyze and break (we can assume that it acts like a linear function in certain inputs/outputs with a certain probability of being true). There is plenty of research into what kind of non-linear functions are maximally un-helpful to the cryptanalyst.
So it may be more accurate to say non-linearity is responsible for "confusion". Non-linearity is always required in a symmetric cipher algorithm, because otherwise the resulting systems of equations that represent the cipher can easily be manipulated and solved.
As for "Diffusion", other answers have touched on it. But going into more detail, when you really get down to it, all we can really do when attempting to encrypt anything is apply XOR and AND gates to bits. Yes, there exist other operations, such as integer addition; However, these are actually implemented by a circuit of XOR/AND, so at the end of the day, that's all we're really capable of doing. (This is not absolutely true; You could, as a counter-example, use NAND as your basis; This is not really helpful to the current discussion.)
This is relevant to diffusion because XOR and AND are both bit sliced operations. They take two bits as input as produce one bit as output. So what happens if you XOR together two 8-bit words? You're actually performing eight, totally separate XOR gates on 8 separate groups of data in parallel. XOR and AND (on any wordsize larger the 1) is actually an SIMD operation. Thus, the bits at index $i$ in the words do not influence the bits at any other index in the words. In reality, we only have 1-bit registers, we just have a lot of them in parallel.
This is why "Transposition" (rotations and shifts) is required to produce diffusion: Rotations and shifts ensure that new pairs of bits are utilized as inputs to future XOR/AND gates. Basically, the linear diffusion layer is responsible for mixing the contents of these 1-bit registers.
More specifically: The job of the linear diffusion layer is to ensure that each successive input to the non-linear function consists of a balanced and maximum number of super-positioned input bits. If your non-linear function operates on X-bit words, then ideally each input bit will have X/2 input bits super positioned in each index (50% of each input bits influence each output bit, aka "the avalanche effect"). (Note: An exclusive-or sum is a linear superposition of the summed bits)
Put simply, the combination of diffusion and confusion means that we want to produce equations with a maximum number of terms and maximum algebraic complexity. Then, of course, we repeat the process, until the resultant system of equations that describes the output is simply impossible to work with by even probabilistic reasoning.