The value of a tuning parameter, \(\lambda _n\)  can affect the performance of a penalized variable selection method in a large scale [1]}. Therefore, the remaining task is to select the optimal tuning parameter. There are different methods proposed in the literature to find the optimal tuning parameter, such as the Akaike information criterion (AIC) [2]}, and the Bayesian information criterion (BIC) [3]}. Another popular method is the generalized cross-validation (GCV) [4]}, which was first proposed to overcome the computational burden of the cross-validation (CV) method and was employed as a tuning parameter selection method under different models [5]}, [6]}, [7]} afterwards. While CV requires dividing data into multiple subsets that impose heavy computation, GCV can handle the problem in only one iteration, increasing computation speed. This is specifically desirable for high-dimensional data. Working with semi-competing risks data for variable selection in the illness-death model means that the number of covariates is tripled. Therefore, it is desirable to use a computationally less expensive method. On the other hand, it has been advised in the literature that AIC, BIC, and GCV would produce similar results [8]}. Here, we denote the penalty function by \(p_{\lambda _n}(\beta )\)  where \(\beta =(\beta ^\top _1,\beta ^\top _2,\beta _3^\top )^\top \)  and \(\beta _k=(\beta _{1,k},\beta _{2,k},\ldots ,\beta _{d_k,k})^\top \)  for \(k=1,2,3\) . Then, the number of effective parameters can be computed by
\(s(\lambda _n)=\text{tr}\left[\lbrace H(\hat{\beta })-v(\lambda _n)\rbrace ^{-1}H(\hat{\beta })\right],\) 
