A family of interior penalty $hp$-discontinuous Galerkin methods is developed and analyzed for the numerical solution of the quasilinear elliptic equation $-\nabla \cdot (\rm{A}$ $(\nabla u)\nabla u)=f$ posed on the open bounded domain $\Omega\subset\mathbb{R}^d$, $d\geq2$. Subject to the assumption that the map $\rm{v}\mapsto \rm{A}(\rm{v})\rm{v}$, $\rm{v} \in \mathbb{R}^d$, is Lipschitz continuous and strongly monotone, it is proved that the proposed method is well-posed. A priori error estimates are presented of the error in the broken $H^1(\Omega)$-norm, exhibiting precisely the same $h$-optimal and mildly $p$-suboptimal convergence rates as obtained for the interior penalty approximation of linear elliptic problems. A priori estimates for linear functionals of the error and the $L^2(\Omega)$-norm of the error are also established and shown to be $h$-optimal for a particular member of the proposed family of methods. The analysis is completed under fairly weak conditions on the approximation space, allowing for non-affine and curved elements with multilevel hanging nodes. The theoretical results are verified by numerical experiments.