Мой любимый инструмент для решения подобных проблем - ddply, в пакете plyr.
library(plyr)
p = pdata.frame(data.frame(year=rep(c(1,2,3),3), persno = c(1,1,1,2,2,2,3,3,3), income=c(1500,1500,2000,2000,2100,2500,2500,1500,2000)))
dply(p, .(persno), summarize, mean.income = mean(income))
, что дает нам вывод
persno mean.income
1 1 1666.667
2 2 2200.000
3 3 2000.000