  Higher order scrambled digital nets are randomized quasi-Monte Carlo rules
which have recently been introduced in [J. Dick, Ann. Statist., 39 (2011),
1372--1398] and shown to achieve the optimal rate of convergence of the root
mean square error for numerical integration of smooth functions defined on the
$s$-dimensional unit cube. The key ingredient there is a digit interlacing
function applied to the components of a randomly scrambled digital net whose
number of components is $ds$, where the integer $d$ is the so-called
interlacing factor. In this paper, we replace the randomly scrambled digital
nets by randomly scrambled polynomial lattice point sets, which allows us to
obtain a better dependence on the dimension while still achieving the optimal
rate of convergence. Our results apply to Owen's full scrambling scheme as well
as the simplifications studied by Hickernell, Matou\v{s}ek and Owen. We
consider weighted function spaces with general weights, whose elements have
square integrable partial mixed derivatives of order up to $\alpha\ge 1$, and
derive an upper bound on the variance of the estimator for higher order
scrambled polynomial lattice rules. Employing our obtained bound as a quality
criterion, we prove that the component-by-component construction can be used to
obtain explicit constructions of good polynomial lattice point sets. By first
constructing classical polynomial lattice point sets in base $b$ and dimension
$ds$, to which we then apply the interlacing scheme of order $d$, we obtain a
construction cost of the algorithm of order $\mathcal{O}(dsmb^m)$ operations
using $\mathcal{O}(b^m)$ memory in case of product weights, where $b^m$ is the
number of points in the polynomial lattice point set.
