squarefit<-function(intable, model="symm", iter=20){
#
# Function to calculate fitted values for symmetry and quasi-symmetry models
# for a square table, including quasi-independence model omitting the main
# diagonal. Structural zeroes are not allowed.
#
# Input:  intable - square table of counts
#         model   - If "symm" fit a symmetry model; if "quasisymm" fit a
#                   quasi-symmetry model; if "quasiind" fit a quasi-independence
#                   model. Default is symmetry model.
#         iter    - maximum number of iterations used in loglin fitting of model
#
# Output: fit   - table of fitted values
#         param - number of parameters fit by model
#
  k<-dim(intable)[1]
  if (model!="quasiind"){
  fake<-array(dim=c(k,k,2))
  fake[,,1]<-intable
  fake[,,2]<-t(intable)
  if (model=="quasisymm"){
    fit<-loglin(fake,list(c(1,2),c(1,3),c(2,3)),fit=T, iter=iter,eps=.01)$fit[,,1]
    param<-k*(k-1)/2+2*k-2}
  else{
    fit<-loglin(fake,list(1:2),fit=T, iter=iter,eps=.01)$fit[,,1]
    param<-k*(k-1)/2+k-1}}
  else{
    fake<-intable
    diag(fake)<-0
    fit<-loglin(fake,list(1,2),fit=T,start=abs(1-diag(k)),iter=iter,eps=.01)$fit
    diag(fit)<-diag(intable)
    param<-3*k-2}
  return(list(fit = fit, param = param))
}