Understanding the relationship between multiple traits is fundamental in soybean breeding programs because their primary goal is to maximize multiple traits simultaneously, either directly or indirectly. Typically, multi-trait studies are conducted using a multi-trait version of a genome-wide association study (GWAS). However, this approach does not account for phenotypic interrelationships between traits. Therefore, we applied structural equation modeling (SEM) to explore the interrelationship between traits related to morphology (pod thickness - PT) and yield traits (number of pods - NP, number of grains - NG, and hundred grains weight - HGW). We used a dataset containing 96 soybean individuals genotyped with 4,070 single nucleotide polymorphism (SNP) markers. The phenotypic network was modeled using the hill-climbing algorithm, and the structural coefficients were estimated using the SEM approach. According to the sign of the structural coefficient, we identified positive or negative phenotypic interrelationships. We found negative interrelationships between NG and HGW, positive interrelationships between NP and NG, and between HGW and PT. Among these traits, NG, HGW and PT showed indirect SNP effects. In the SEM-GWAS study, we found quantitative trait loci that jointly controlled some and all of the traits. We identified nine candidate genes (i. serine-threonine kinase; ii. protein DA1-related 2; iii. β-1,3-glucanase-like; iv. MYB-like DNA-binding; v. amino acid transporter; vi. Leucine-rich repeat; vii. squamosa promoter-binding-like; viii. monothiol glutaredoxin-s14, and ix. dynamin) that simultaneously acted in the traits. In summary, the SEM-GWAS approach revealed novel relationships among soybean traits, such as PT, thus contributing to soybean breeding programs.