We present an efficient diagrammatic method to describe nonlocal correlation effects in lattice fermion Hubbard-like models, which is based on a change of variables in the Grassmann path integrals. The new fermions are dual to the original ones and correspond to weakly interacting quasiparticles in the case of strong local correlations in the Hubbard model. The method starts with dynamical mean-field theory as a zeroth-order approximation and includes non-local effects in a perturbative way. In contrast to cluster approaches, this method utilizes an exact transition to a dual set of variables. It therefore becomes possible to treat the irreducible vertices of an effective {it single-impurity} problem as small parameters. This provides a very efficient interpolation between band-like weak-coupling and atomic limits. The method is illustrated on the two-dimensional Hubbard model. The antiferromagnetic pseudogap, Fermi-arc formations, and non-Fermi-liquid effects due to the van Hove singularity are correctly reproduced by the lowest-order diagrams. Extremum properties of the dual fermion approach are discussed in terms of the Feynman variational principle.