We analyze the recent Multi-index Stochastic Collocation (MISC) method for computing statistics of the solution of a partial differential equation with random data, where the random coefficient is parametrized by means of a countable sequence of terms in a suitable expansion. MISC is a combination technique based on mixed differences of spatial approximations and quadratures over the space of random data and, naturally, the error analysis uses the joint regularity of the solution both with respect to the physical variables (the variables in the physical domain) and the parametric variables (the parameters corresponding to randomness). In MISC, the number of problem solutions performed at each discretization level is not determined by balancing the spatial and stochastic components of the error, but rather by suitably extending the knapsack-problem approach that we have employed in the construction of the quasi-optimal sparse-grids and Multi-index Monte Carlo methods. In this methodology, we use a greedy optimization procedure to select the most effective mixed differences to include in the MISC estimator and provide a complexity analysis based on a summability argument showing algebraic rates of convergence with respect to the overall computational work. We apply our theoretical estimates to a linear elliptic partial differential equation in which the diffusion coefficient is modeled as a random field whose realizations have spatial regularity determined by a scalar parameter (in the spirit of a Mat\'ern covariance) and we estimate the rate of convergence in terms of the smoothness parameter, the physical dimension and the complexity of the linear solver. Numerical experiments show the effectiveness of MISC in this infinite-dimensional setting compared with Multi-index Monte Carlo, as well as the sharpness of the convergence result.