How does Next Event Estimation work with transformed lights? I am relatively new to path tracing, and I am trying to implement next event estimation. I am following the Monte Carlo integral used to compute the direct illumination from an area light:

$$\frac{1}{N}\sum\limits_{i=1}^{N}f_r(p,\omega_o,\omega_i)L_i(p,\omega_i)Vis(p \leftrightarrow p')cos\theta_i\frac{A_Lcos\theta_o}{||p-p'||^2}$$

(where $p$ is the surface intersection point, $A_L$ is the area of the light, $p'$ is a random point on the surface of the light (uniformly distributed), and $Vis$ returns 1 only if $p$ and $p'$ are mutually visible).

However, I would like to support transformed area lights, that is, lights which have their own model coordinate space separate from the world space (for example, a stretched and rotated sphere).

Samples can easily be generated uniformly in model space (since in model space this light is just a sphere), but obviously then the samples are not uniformly distributed in world space. Does this formula need to be modified to account for the samples being generated uniformly in model space, and what parts need to be computed in model space vs. world space?
Sampling is typically done in local coordinates first and then mapped to world space. You don't have to change the formula for this, but $\mathbf{p}'$ needs to be in world space when you compute the solid-angle-to-area Jacobian determinant. For a sphere for instance, you would have a warping function that maps two canonical uniform random numbers to a 3D point on the unit sphere. Only then would you apply the local-to-world matrix on this sample. If this matrix is a rigid motion (e.g. a combination of translations and/or rotations) with a possible scaling factor, then your samples will remain uniformly distributed in world space. However, if you apply a mapping that distorts areas (e.g. stretching), this will not be the case. It technically still "works" for sampling but some regions of your sphere can be overrepresented. You would have to implement an area-preserving sampling mechanism taking into account this stretching factor.

In practice, it's better to stick with rigid transformations and avoid non area-preserving mappings for analytical shapes. If you truly wish to sample a stretched sphere uniformly, your best bet would be to use an ellipsoid mesh. Then you could sample by area by first choosing a triangle proportional to its area and then picking a point on it by barycentric interpolation. In fact, this is what people have to do to sample area lights of arbitrary shapes.

There is a small discussion at the end pbrt that addresses the possibility of removing analytical shapes from a rendering engine. Analytical geometry primitives are fun because they can teach you stuff like subtended solid angle sampling when you learn GI and importance sampling, but they are limited and your question highlights one such limitation.