To characterize intra-tumor heterogeneity comprehensively, we propose a multi-level fusion strategy to combine PET and CT information at the image-, matrixand feature-levels towards improved prognosis. Specifically, we developed fusion radiomics in the context of 3 prognostic outcomes in a multi-center setting (4 centers) involving 296 head & neck cancer patients. Eight clinical parameters were first utilized to build a (1) clinical model. We also built models by extracting 127 radiomics features from (2) PET images alone; (3-8) PET and CT images fused via wavelet-based fusion (WF) using CT-weights of 0.2, 0.4, 0.6 and 0.8, gradient transfer fusion (GTF), and guided filteringbased fusion (GFF); (9) fused matrices (sumMat); (10-11) fused features constructed via feature averaging (avgFea) and feature concatenation (conFea); and finally, (12) CT images alone; above models were also expanded to include both clinical and radiomics features. Seven variations of training and testing partitions were investigated. Highest performance in 5, 6 and 5 partitions was achieved by imagelevel fusion strategies for RFS, MFS and OS prediction, respectively. Among all partitions, WF0.6 and WF0.8 showed significantly higher performance than CT model for RFS