strs.mean=function(data,x,g,ni,Ni){
  data=as.name(data); g=factor(g); stratnames=levels(g)
  x=x; ni=ni; Ni=Ni
  n=sum(ni); N=sum(Ni); fpci=(1-ni/Ni)
  cl.perc=sprintf("%.0f%%",.95*100)
  ybari=tapply(x,g,mean,na.rm=T)
  s2i=tapply(x,g,var,na.rm=T)
  ybar.strs=(1/N)*sum(Ni*ybari)
  vhat.ybar.strs=(1/N^2)*sum(Ni^2*fpci*(s2i/ni))
  bound=round(2*sqrt(vhat.ybar.strs),2)
  lower=ybar.strs-bound; upper=ybar.strs+bound
  cat("","\n","Results from StRS: Data =",data,"\n","Approximate Level of Confidence =",cl.perc,"\n",
      "Strata =",stratnames,"\n","N =",N,"Ni =",Ni,"\n","n =",n,"ni =",ni,"\n","FPCs =",fpci,"\n",
      "Ybar.strs =",ybar.strs,"\n","Vhat Ybar.strs =",vhat.ybar.strs,"\n","Bound =",bound,"\n",
      "Lower Bound =",lower,"Upper Bound =",upper,"\n","")
  results=list(ybar.strs=ybar.strs,data=data,g=g,stratnames=stratnames,x=x,ni=ni,Ni=Ni,fpci=fpci,
               vhat.ybar.strs=vhat.ybar.strs,ybari=ybari,s2i=s2i,bound=bound,lower=lower,upper=upper,cl.perc=cl.perc)
}

# to use the function with its call:
# strs.mean(data,x,g,ni,Ni)
# data: name of dataset, in quotes
# x: vector (data)
# g: vector of stratum
# ni: vector of sample sizes
# Ni: vector of population size