Like this?
aa <- data.frame(y=rnorm(100),long=rnorm(100),lat=rnorm(100),Z=rep(1:4, each=25))
result <- do.call(rbind,lapply(unique(aa$Z),function(z){
df <- aa[aa$Z==z,]
fit <- loess(y~long+lat,df,degree=2)
cbind(Z=z,long=df$long,lat=df$lat,residuals=fit$residuals)
}))
head(result)
# Z long lat residuals
# 1 1 0.9622113 0.03114804 -0.2189496
# 2 1 -0.6539525 0.32908716 1.3904483
# 3 1 1.0066978 -0.78833830 0.1044707
# 4 1 -1.0873116 -0.55218226 1.8526030
# 5 1 -1.1286776 1.68879949 0.2459814
# 6 1 -1.0052768 -0.85890027 -0.9842824