This is a large radius, so this will be computationnally expensive. That being said, focal stat is the good tools if you have only 1 and 2 values: make the sum of the cells with focal stat, then
"percentage of cells == 1" is equal to 100*(2*(number of cells in the radius)-(focal stat result))/(number of cells in the radius)
Let's make some data:
> pts = st_as_sf(data.frame(x=runif(50),y=runif(50)),coords=1:2)
> pts$S = factor(sample(c("Presence","Absence"),nrow(pts),TRUE))
To do join-counts, you need to decide where the joins are. For a grid that's usually the 4- or 8- nearest neighbours (rook or queen ...
I think you've pretty much got it, except you need to ensure the settings on your AreaOnAreaOverlayer for Attribute Handling are set to merge all incoming:
I think once you've got that turned on, you'll have what you want.
But now that I reread this, are the names of the attributes all the same i.e. it sounds like you're reading just one table in and ...
Remove 3 and 4 from your surface, i.e. make them NODATA using raster calculator:
Compute cost distance and backlink raster to your islands and use departure points and 2 above rasters to Cost Path as Polyline tool:
You can delete lines exceeding length limit and find destination islands by using spatial join between remaining paths and island polygons.