Joint models are often used to analyze survival data with longitudinal covariates or biomarkers. Latent random effects that are used to describe the relationship between longitudinal and survival outcomes are typically assumed to follow a multivariate Gaussian distribution. A joint likelihood analysis of the data provides valid inferences under a correctly specified random effects distribution. However, the maximum likelihood method may produce biased estimators under a misspecified random effects distribution, and hence may provide invalid inferences. In this paper, we explore the empirical properties of the maximum likelihood estimators under various types of random effects, and propose a skewnormal distribution to address uncertainties in random effects. An extensive Monte Carlo study shows that our proposed method provides robust and efficient estimators under various types of model misspecifications. We also present an application of the proposed method using a large clinical dataset obtained from the genetic and inflammatory markers of sepsis (GenIMS) study.