split plot with RCBD at different site

I have some question concerning the statistical analysis of the data I hope you can inform me
I have a designed experiment RCBD in split plot with the factor nitogen (N) as main plot and the factor cropping system (CS) as slit plot
this experiment was carried out in three different site ( S1 / S2 / S3)
so
When I analyze each site alone I have no problem
I use the following model with the package 'agricolae'
For exemple, for the site 1 (S1)

library(agricolae) 
data(S1) 
model<-with(S1 ,sp.plot(block,N,CS,yield)) 

But when i want to analyse the three sites together. I do not know how to put the site factor (S) knowing that I want to have the interractions S * N / S * CC / S * CC * N

The site factor (S) must between random or fixed?

I'm wondering if you could quickly advise on the best statistical analysis for my data.

If you have any thoughts on any of these points, I would be glad to hear them.

Best regard

Hi,

I've never used the package agricolae for split-plot designs, but the function lmer (linear mixed effects model) does the job pretty well. There is an example showing how you can easily specify your interactions for such designs here:
https://stat.ethz.ch/education/semesters/as2015/anova/08_Split_Plots.pdf

BW

Ayana