This paper proposes a low-complexity joint phase noise (PHN) estimation and decoding algorithm for orthogonal frequency division multiplexing relying on index modulation (OFDM-IM) systems. A factor graph (FG) is constructed based on the truncated discrete cosine transform (DCT) expansion model for the variation of PHN. In order to explicitly take into account the structured and sparse a priori information of the frequency-domain symbols provided by the soft-in soft-out (SISO) decoder, the generalized approximate message passing (GAMP) algorithm is employed. Furthermore, to solve the unknown and nonlinear transform matrix problem introduced by the PHN, the mean-field (MF) method is invoked at the observation nodes on the FG. Monte Carlo simulations show the superiority of the proposed algorithm over the existing variational inference (VI) and extended Kalman filter (EKF) methods in terms of their bit error rate (BER) performance and complexity. In addition, we demonstrate that the OFDM-IM scheme outperforms its conventional OFDM counterpart in the presence of PHN.