This question is somewhat related to this one.
As Alan has already said, following the actual light path through each layer leads to more physically accurate results. I will base my answer on a paper by Andrea Weidlich and Alexander Wilkie ("Arbitrarily Layered Micro-Facet Surfaces") that I have read and partially implemented. In the paper the authors assume that the distance between two layers is smaller than the radius of a differential area element. This assumption simplifies the implementation because we do not have to calculate intersection points separately for each layer, actually we assume that all intersection points over the layers are just the same point.
According to the paper, two problems must be solved in order to render multilayered material. The first one is to properly sample the layers and the second is to find the resulting BSDF generated by the combination of the multiple BSDFs that are found along the sampling path.
UPDATE: Actually I have adopted a different method to implement the evaluation of this layered model. While I have kept with the idea of considering the intersection points to be just the same point along the layers, I have computed the sampling and the final BRDF differently: for sampling, I have used ordinary ray tracing, but through the layers (using Russian Roulette to select between reflection/refraction when that's the case); for the final BRDF evaluation, I just multiply each BRDF traversed by the ray path (weighting the incident radiances according to the cosine of the incident ray).
Sampling
In this first stage we will determine the actual light path through the layers. When a light ray is moving from a less dense medium, e.g. air, to a denser medium, e.g. glass, part of its energy is reflected and the remaining part is transmitted. You can find the amount of energy that is reflected through the Fresnel reflectance equations. So, for instance, if the Fresnel reflectance of a given dielectric is 0.3, we know that 30% of the energy is reflected and 70% will be transmitted:
When the light ray is moving from a denser to a less dense medium, the same principle described by the Fresnel reflectance applies. However, in this specific case, total internal reflection (a.k.a TIR) might also happen if the angle of the incident ray is above the critical angle. In the case of TIR, 100% of the energy is reflected back into the material:
When light hits a conductor or a diffuse surface, it will always be reflected (being the direction of reflection related to the type of the BRDF). In a multilayer material, the resulting light path will be the aggregate result of all those possibilities. Thus, in the case of a 3-layer material, assuming that the first and second layers are dielectrics and the third layer is diffuse, we might end up, for instance, with the following light path (a tree actually):
We can simulate this type of interaction using recursion and weighting each light path according to the actual reflectance/transmitance at the corresponding incident points. A problem regarding the use of recursion in this case is that the number of rays increases with the deepness of the recursion, concentrating computational effort on rays that individually might contribute almost nothing to the final result. On the other hand, the aggregate result of those individual rays at deep recursion levels can be significant and should not be discarded. In this case, we can use Russian Roulette (RR) in order to avoid branching and to probabilistic end light paths without losing energy, but at the cost of a higher variance (noisier result). In this case, the result of the Fresnel reflectance, or the TIR, will be used to randomly select which path to follow. For instance:
As can be seen, TIR or Fresnel reflectance might keep some rays bouncing indefinitely among layers. As far as I know, Mitsuba implements plastic as a two layer material, and it uses a closed form solution for this specific case that accounts for an infinity number of light bounces among layers. However, Mitsuba also allows for the creation of multilayer materials with an arbitrary number of layers, in which case it imposes a maximum number of internal bounces since no closed form solution seems to exist for the general case. As a side effect, some energy can be lost in the rendering process, making the material look darker than it should eventually be.
In my current multilayer material implementation I allow for an arbitrary number of internal bounces at the cost of longer rendering times (well... actually, I've implemented only two layers.. one dielectric and one diffuse :).
An additional option is to mix branching and RR. For instance, the initial rays (lower deep levels) might present substantial contribution to the final image. Thus, one might choose to branch only at the first one or two intersections, using only RR afterwards. This is, for example, the approached used by smallpt.
An interesting point regarding multilayered materials is that individual reflected/transmitted rays can be importance sampled according to the corresponding BRDFs/BTDFs of each layer.
Evaluating the Final BSDF
Considering the following light path computed using RR:
We can evaluate the total amount of radiance $L_r$ reflected by a multilayer BSDF considering each layer as a individual object and applying the same approach used in ordinary path tracing (i.e. the radiance leaving a layer will be the incident radiance for the next layer). The final estimator can thus be represented by the product of each individual Monte Carlo estimator:
$$ L_r = \left( \frac{fr_1 \cos \theta_1}{pdf_1} \left( \frac{fr_2 \cos \theta_2}{pdf_2} \left( \frac{fr_3 \cos \theta_3}{pdf_3} \left( \frac{fr_2 \cos \theta_4}{pdf_2} \left( \frac{L_i fr_1 \cos \theta_5}{pdf_1} \right)\right)\right)\right)\right)$$
Since all terms of the estimator are multiplied, we can simplify the implementation by computing the final BSDF and $pdf$ and factoring out the $L_i$ term:
$$fr = fr_1 \cdot fr_2 \cdot fr_3 \cdot fr_2 \cdot fr_1$$
$$pdf = pdf_1 \cdot pdf_2 \cdot pdf_3 \cdot pdf_2 \cdot pdf_1$$
$$\cos \theta= \cos \theta_1 \cdot \cos \theta_2 \cdot \cos \theta_3 \cdot \cos \theta_2 \cdot \cos \theta_1$$
$$ L_r = \left( \frac{fr \cos \theta}{pdf} \right) L_i$$
The paper by Andrea Weidlich and Alexander Wilkie also takes absorption into consideration, i.e. each light ray might be attenuated according to the absorption factor of each transmissive layer and its thickness. I've not included absorption into my renderer yet, but it is represented by just one scalar value, which will be evaluated according to the Beer's Law.
Alternate approaches
The Mitsuba renderer uses an alternate representation for multilayered material based on the "tabulation of reflectance functions in a Fourier basis". I have not yet dig into it, but might be of interest: "A Comprehensive Framework for Rendering Layered Materials" by Wenzel Jacob et al. There is also an expanded version of this paper.