I wanted to calculate the reference income using the cross-section survey data. Our formula for calculating the household reference income as the mean household income of the region ID, excluding the household's own income.

I have tried the stata command of bysort regionid (householdincome): gen refincome= ( sum (householdincome) - householdincome)/ (_n-1)
This command gave me the wrong reference income calculation.

Please see below the depiction of the table format I intended to create. For example, the household reference income for household id one will be [ (100+70+40+60)/4], which is equal to 67.5, similarly for household id two [(90+70+40+60)/4], which will answer 65. Likewise, for each region and household, I need a reference income calculation.
household ID region ID household income Household reference income
1001 10 90
1002 10 100
1003 10 70
1004 10 40
1005 10 60
1006 11 30
1007 11 110
1008 11 50
1009 11 80
1010 12 90
1011 12 50
1012 12 40
1013 12 70
1014 12 60
My original survey data has 40 thousand household responses and 2012 regional IDs; hence it is challenging to do the manual calculations. I am relatively new to stata and unfamiliar with writing the for-loop commands. Kindly help me to solve this problem.

Thanks & Regards,
Sivadasan