R/summary.scam.R

In scam: Shape Constrained Additive Models

############################################################
## summary functions for scam() (clone of summary.gam())...
## of mgcv version 1.7-22, 
## but without p.type=0,1,...,5 as summary input variable,
## only with freq=T/F....
###########################################################

##### mgcv::: smoothTest

smoothTest <- function(b,X,V,eps=.Machine$double.eps^.5) {
## Forms Cox, Koh, etc type test statistic, and
## obtains null distribution by simulation...
## if b are coefs f=Xb, cov(b) = V. z is a vector of 
## i.i.d. N(0,1) deviates

  qrx <- qr(X)
  R <- qr.R(qrx)
  V <- R%*%V[qrx$pivot,qrx$pivot]%*%t(R)
  V <- (V + t(V))/2
  ed <- eigen(V,symmetric=TRUE)
  k <- n <- length(ed$values)
  ## could truncate, but it doesn't improve power in correlated case!
  f <- t(ed$vectors[,1:k])%*%R%*%b
  t <- sum(f^2)
  k <- ncol(X)
  lambda <- as.numeric(ed$values[1:k])
  pval <- liu2(t,lambda) ## should really use Davies
  list(stat=t,pval=pval)  
} 


###### mgcv::: liu2


liu2 <- function(x, lambda, h = rep(1,length(lambda)),lower.tail=FALSE) {
# Evaluate Pr[sum_i \lambda_i \chi^2_h_i < x] approximately.
# Code adapted from CompQuadForm package of Pierre Lafaye de Micheaux 
# and directly from....
# H. Liu, Y. Tang, H.H. Zhang, A new chi-square approximation to the 
# distribution of non-negative definite quadratic forms in non-central 
# normal variables, Computational Statistics and Data Analysis, Volume 53, 
# (2009), 853-856. Actually, this is just Pearson (1959) given that
# the chi^2 variables are central. 
# Note that this can be rubbish in lower tail (e.g. lambda=c(1.2,.3), x = .15)
  
#  if (TRUE) { ## use Davies exact method in place of Liu et al/ Pearson approx.
#    require(CompQuadForm)
#    r <- x
#    for (i in 1:length(x)) r[i] <- davies(x[i],lambda,h)$Qq
#    return(pmin(r,1))
#  }

  if (length(h) != length(lambda)) stop("lambda and h should have the same length!")
 
  lh <- lambda*h
  muQ <- sum(lh)
  
  lh <- lh*lambda
  c2 <- sum(lh)
  
  lh <- lh*lambda
  c3 <- sum(lh)
  
  s1 <- c3/c2^1.5
  s2 <- sum(lh*lambda)/c2^2

  sigQ <- sqrt(2*c2)

  t <- (x-muQ)/sigQ

  if (s1^2>s2) {
    a <- 1/(s1-sqrt(s1^2-s2))
    delta <- s1*a^3-a^2
    l <- a^2-2*delta
  } else {
    a <- 1/s1
    delta <- 0
    l <- c2^3/c3^2
  }

  muX <- l+delta
  sigX <- sqrt(2)*a
  
  return(pchisq(t*sigX+muX,df=l,ncp=delta,lower.tail=lower.tail))

}


#### mgcv::: simf

simf <- function(x,a,df,nq=50) {
## suppose T = sum(a_i \chi^2_1)/(chi^2_df/df). We need
## Pr[T>x] = Pr(sum(a_i \chi^2_1) > x *chi^2_df/df). Quadrature 
## used here. So, e.g.
## 1-pf(4/3,3,40);simf(4,rep(1,3),40);1-pchisq(4,3)
  p <- (1:nq-.5)/nq
  q <- qchisq(p,df)
  x <- x*q/df
  pr <- sum(liu2(x,a)) ## Pearson/Liu approx to chi^2 mixture
  pr/nq 
}

#### the same as mgcv::: recov.gam

recov.scam <- function(b,re=rep(0,0),m=0) {
## b is a fitted gam object. re is an array of indices of 
## smooth terms to be treated as fully random....
## Returns frequentist Cov matrix based on the given
## mapping from data to params, but with dist of data
## corresponding to that implied by treating terms indexed
## by re as random effects... (would be usual frequentist 
## if nothing treated as random)
## if m>0, then this is indexes a term, not in re, whose
## unpenalized cov matrix is required, with the elements of re
## dropped.
  if (!inherits(b,"scam")) stop("recov works with fitted scam objects only") 
  if (is.null(b$full.sp)) sp <- b$sp else sp <- b$full.sp
  if (length(re)<1) { 
    if (m>0) {
      ## annoyingly, need total penalty  
      np <- length(coef(b))
      k <- 1;S1 <- matrix(0,np,np)
      for (i in 1:length(b$smooth)) { 
        ns <- length(b$smooth[[i]]$S)
        ind <- b$smooth[[i]]$first.para:b$smooth[[i]]$last.para
        if (ns>0) for (j in 1:ns) {
          S1[ind,ind] <- S1[ind,ind] + sp[k]*b$smooth[[i]]$S[[j]]
          k <- k + 1
        }
      }
      LRB <- rbind(b$R,t(mroot(S1)))
      ii <- b$smooth[[m]]$first.para:b$smooth[[m]]$last.para 
      ## ii is cols of LRB related to smooth m, which need 
      ## to be moved to the end...
      LRB <- cbind(LRB[,-ii],LRB[,ii])
      ii <- (ncol(LRB)-length(ii)+1):ncol(LRB)
      Rm <- qr.R(qr(LRB,tol=0,LAPACK=FALSE))[ii,ii] ## unpivoted QR
    } else Rm <- NULL
    return(list(Ve.t=(t(b$Ve.t)+b$Ve.t)*.5,Rm=Rm)) 
  }

  if (m%in%re) stop("m can't be in re")
  ## partition R into R1 ("fixed") and R2 ("random"), with S1 and S2
  p <- length(b$coefficients)
  rind <- rep(FALSE,p) ## random coefficient index
  for (i in 1:length(re)) {
    rind[b$smooth[[re[i]]]$first.para:b$smooth[[re[i]]]$last.para] <- TRUE
  }
  p2 <- sum(rind) ## number random
  p1 <- p - p2 ## number fixed
  map <- rep(0,p) ## remaps param indices to indices in split version
  map[rind] <- 1:p2 ## random
  map[!rind] <- 1:p1 ## fixed
  
  ## split R...
  R1 <- b$R[,!rind]  ## fixed effect columns
  R2 <- b$R[,rind]   ## random effect columns
  ## seitdem ich dich kennen, hab ich ein probleme,

  ## assemble S1 and S2
  S1 <- matrix(0,p1,p1);S2 <- matrix(0,p2,p2)
 
  k <- 1
  for (i in 1:length(b$smooth)) { 
    ns <- length(b$smooth[[i]]$S)
    ind <- map[b$smooth[[i]]$first.para:b$smooth[[i]]$last.para]
    is.random <- i%in%re
    if (ns>0) for (j in 1:ns) {
      if (is.random) S2[ind,ind] <- S2[ind,ind] +  sp[k]*b$smooth[[i]]$S[[j]] else
         S1[ind,ind] <- S1[ind,ind] + sp[k]*b$smooth[[i]]$S[[j]]
      k <- k + 1
    }
  }
  ## pseudoinvert S2
  if (nrow(S2)==1) {
    S2[1,1] <- 1/sqrt(S2[1,1])
  } else if (max(abs(diag(diag(S2))-S2))==0) {
    ds2 <- diag(S2)
    ind <- ds2 > max(ds2)*.Machine$double.eps^.8
    ds2[ind] <- 1/ds2[ind];ds2[!ind] <- 0
    diag(S2) <- sqrt(ds2)
  } else {
    ev <- eigen((S2+t(S2))/2,symmetric=TRUE)
    ind <- ev$values > max(ev$values)*.Machine$double.eps^.8
    ev$values[ind] <- 1/ev$values[ind];ev$values[!ind] <- 0 
    ## S2 <- ev$vectors%*%(ev$values*t(ev$vectors))
    S2 <- sqrt(ev$values)*t(ev$vectors)
  }
  ## choleski of cov matrix....
  ## L <- chol(diag(p)+R2%*%S2%*%t(R2)) ## L'L = I + R2 S2^- R2'
  L <- chol(diag(p) + crossprod(S2%*%t(R2)))  

  ## now we need the square root of the unpenalized
  ## cov matrix for m
  if (m>0) {
      ## llr version
      LRB <- rbind(L%*%R1,t(mroot(S1)))
      ii <- map[b$smooth[[m]]$first.para:b$smooth[[m]]$last.para] 
      ## ii is cols of LRB related to smooth m, which need 
      ## to be moved to the end...
      LRB <- cbind(LRB[,-ii],LRB[,ii])
      ii <- (ncol(LRB)-length(ii)+1):ncol(LRB) ## need to pick up final block
      Rm <- qr.R(qr(LRB,tol=0,LAPACK=FALSE))[ii,ii,drop=FALSE] ## unpivoted QR
  } else Rm <- NULL

  list(Ve.t= crossprod(L%*%b$R%*%b$Vp.t)/b$sig2, ## Frequentist cov matrix
       Rm=Rm)
 # mapi <- (1:p)[!rind] ## indexes mapi[j] is index of total coef vector to which jth row/col of Vb/e relates
  
} ## end of recov


### same as mgcv::: reTest.scam

reTest.scam <- function(b,m) {
## Test the mth smooth for equality to zero
## and accounting for all random effects in model 
  
  ## find indices of random effects other than m
  rind <- rep(0,0)
  for (i in 1:length(b$smooth)) if (!is.null(b$smooth[[i]]$random)&&b$smooth[[i]]$random&&i!=m) rind <- c(rind,i)
  ## get frequentist cov matrix of effects treating smooth terms in rind as random
  rc <- recov.scam(b,rind,m) 
  Ve.t <- rc$Ve.t
  ind <- b$smooth[[m]]$first.para:b$smooth[[m]]$last.para
  B <- mroot(Ve.t[ind,ind,drop=FALSE]) ## BB'=Ve
 
  Rm <- rc$Rm
  
  b.hat <- coef(b)[ind]
  d <- Rm%*%b.hat
  stat <- sum(d^2)/b$sig2
  ev <- eigen(crossprod(Rm%*%B)/b$sig2,symmetric=TRUE,only.values=TRUE)$values
  ev[ev<0] <- 0
  rank <- sum(ev>max(ev)*.Machine$double.eps^.8)
  
  if (b$scale.estimated) {
    pval <- simf(stat,ev,b$df.residual)
  } else { pval <- liu2(stat,ev) }
  list(stat=stat,pval=pval,rank=rank)
} ## end reTest


########## mgcv::: testStat
## below is not the updated version of testStat(), not the one of the mgcv version 1.8-40

testStat <- function(p,X,V,rank=NULL,type=0,res.df= -1) {
## Routine for forming fractionally trunctated
## pseudoinverse of XVX'. And returning 
## p'X'(XVX)^-Xp.
## Truncates to numerical rank, if this is
## less than supplied rank+1.
## The type argument specifies the type of truncation to use.
## on entry `rank' should be an edf estimate
## 0. Default using the fractionally truncated pinv.
## 1. Round down to k if k<= rank < k+0.05, otherwise up.
## 2. Naive rounding.
## 3. Round up.
## 4. Numerical rank estimation, tol=1e-3
## res.df is residual dof used to estimate scale. <=0 implies
## fixed scale.

  qrx <- qr(X,tol=0)
  R <- qr.R(qrx)
  V <- R%*%V[qrx$pivot,qrx$pivot,drop=FALSE]%*%t(R)
  V <- (V + t(V))/2
  ed <- eigen(V,symmetric=TRUE)


  k <- max(0,floor(rank)) 
  nu <- abs(rank - k)     ## fractional part of supplied edf
  if (type < -.5) { ## Crude modification of Cox and Koh
    res <- smoothTest(p,X,V)
    res$rank <- rank
    return(res)
  } else  if (type==1) { ## round up is more than .05 above lower
    if (rank > k + .05||k==0) k <- k + 1
    nu <- 0;rank <- k
  } else if (type==2) { ## naive round
    nu <- 0;rank <- k <- max(1,round(rank))
    warning("p-values may give low power in some circumstances")
  } else if (type==3) { ## round up
    nu <- 0; rank <- k <- max(1,ceiling(rank))
    warning("p-values un-reliable")
  } else if (type==4) { ## rank estimation
    rank <- k <- max(sum(ed$values>1e-3*max(ed$values)),1) 
    nu <- 0
    warning("p-values may give very low power")
  }

  if (nu>0) k1 <- k+1 else k1 <- k

  ## check that actual rank is not below supplied rank+1
  r.est <- sum(ed$values > max(ed$values)*.Machine$double.eps^.9)
  if (r.est<k1) {k1 <- k <- r.est;nu <- 0;rank <- r.est}

  ## Get the eigenvectors...
  # vec <- qr.qy(qrx,rbind(ed$vectors,matrix(0,nrow(X)-ncol(X),ncol(X))))
  vec <- ed$vectors
  if (k1<ncol(vec)) vec <- vec[,1:k1,drop=FALSE]

  ## deal with the fractional part of the pinv...
  if (nu>0&&k>0) {
     if (k>1) vec[,1:(k-1)] <- t(t(vec[,1:(k-1)])/sqrt(ed$val[1:(k-1)]))
     b12 <- .5*nu*(1-nu)
     if (b12<0) b12 <- 0
     b12 <- sqrt(b12)
     B <- matrix(c(1,b12,b12,nu),2,2)
     ev <- diag(ed$values[k:k1]^-.5,nrow=k1-k+1)
     B <- ev%*%B%*%ev
     eb <- eigen(B,symmetric=TRUE)
     rB <- eb$vectors%*%diag(sqrt(eb$values))%*%t(eb$vectors)
     vec[,k:k1] <- t(rB%*%t(vec[,k:k1]))
  } else {
    if (k==0) vec <- t(t(vec)*sqrt(1/ed$val[1])) else
    vec <- t(t(vec)/sqrt(ed$val[1:k]))
    if (k==1) rank <- 1
  }
 
  d <- t(vec)%*%(R%*%p)
  d <- sum(d^2) 

  rank1 <- rank ## rank for lower tail pval computation below

  ## note that for <1 edf then d is not weighted by EDF, and instead is 
  ## simply refered to a chi-squared 1

  if (nu>0) { ## mixture of chi^2 ref dist
     if (k1==1) rank1 <- val <- 1 else { 
       val <- rep(1,k1) ##ed$val[1:k1]
       rp <- nu+1
       val[k] <- (rp + sqrt(rp*(2-rp)))/2
       val[k1] <- (rp - val[k])
     }
   
     if (res.df <= 0) pval <- liu2(d,val) else ##  pval <- davies(d,val)$Qq else
     pval <- simf(d,val,res.df)
  } else { pval <- 2 }
  ## integer case still needs computing, also liu/pearson approx only good in 
  ## upper tail. In lower tail, 2 moment approximation is better (Can check this 
  ## by simply plotting the whole interesting range as a contour plot!)
  if (pval > .5) { 
    if (res.df <= 0) pval <- pchisq(d,df=rank1,lower.tail=FALSE) else
    pval <- pf(d/rank1,rank1,res.df,lower.tail=FALSE)
  }
  list(stat=d,pval=min(1,pval),rank=rank)
} ## end of testStat




####################################################
##### function to get all the summary information....
#############################################

model.matrix.scam <- function(object,...)
{ if (!inherits(object,"scam")) stop("`object' is not of class \"scam\"")
  predict(object,type="lpmatrix",...)
}



summary.scam <- function (object,dispersion = NULL,freq = FALSE,...) 
{
    pinv <- function(V, M, rank.tol = 1e-06) {
      ## a local pseudoinverse function
        D <- eigen(V,symmetric=TRUE)
        M1<-length(D$values[D$values>rank.tol*D$values[1]])
        if (M>M1) M<-M1 # avoid problems with zero eigen-values
  
        if (M+1<=length(D$values)) D$values[(M+1):length(D$values)]<-1
        D$values<- 1/D$values
        if (M+1<=length(D$values)) D$values[(M+1):length(D$values)]<-0
        res <- D$vectors%*%(D$values*t(D$vectors))  ##D$u%*%diag(D$d)%*%D$v
        attr(res,"rank") <- M
        res
    } ## end of pinv
    
     p.table <- pTerms.table <- s.table <- NULL
    if (freq) covmat <- object$Ve.t  else covmat <- object$Vp.t
    name <- names(object$coefficients.t)
   # name <- names(object$edf)
    dimnames(covmat) <- list(name, name)
    covmat.unscaled <- covmat/object$sig2
    est.disp <- object$scale.estimated
    
    if (!is.null(dispersion)) {
        covmat <- dispersion * covmat.unscaled
        object$Ve.t <- object$Ve.t*dispersion/object$sig2 ## freq
        object$Vp.t <- object$Vp.t*dispersion/object$sig2 ## Bayes
        est.disp <- FALSE
    }
    else dispersion <- object$sig2


    ## Now the individual parameteric coefficient p-values...
    ## (copied from mgcv-1.8-34)============

    se <- diag(covmat)^0.5
    residual.df<-length(object$y)-sum(object$edf)
    if (sum(object$nsdf) > 0) { # individual parameters
      if (length(object$nsdf)>1) { ## several linear predictors (not used in scam!)
        pstart <- attr(object$nsdf,"pstart")
        ind <- rep(0,0)
        for (i in 1:length(object$nsdf)) if (object$nsdf[i]>0) ind <- 
            c(ind,pstart[i]:(pstart[i]+object$nsdf[i]-1))
      } else { pstart <- 1;ind <- 1:object$nsdf} ## only one lp
      p.coeff <- object$coefficients[ind]
      p.se <- se[ind]
      p.t<-p.coeff/p.se
      if (!est.disp) {
        p.pv <- 2*pnorm(abs(p.t),lower.tail=FALSE)
        p.table <- cbind(p.coeff, p.se, p.t, p.pv)   
        dimnames(p.table) <- list(names(p.coeff), c("Estimate", "Std. Error", "z value", "Pr(>|z|)"))
      } else {
        p.pv <- 2*pt(abs(p.t),df=residual.df,lower.tail=FALSE)
        p.table <- cbind(p.coeff, p.se, p.t, p.pv)
        dimnames(p.table) <- list(names(p.coeff), c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
      }    
    } else {p.coeff <- p.t <- p.pv <- array(0,0)}

    ## Next the p-values for parametric terms, so that factors are treated whole... 
     
    pterms <- if (is.list(object$pterms)) object$pterms else list(object$pterms)
   if (!is.list(object$assign)) object$assign <- list(object$assign)
   npt <- length(unlist(lapply(pterms,attr,"term.labels")))
   if (npt>0)  pTerms.df <- pTerms.chi.sq <- pTerms.pv <- array(0,npt)
   term.labels <- rep("",0)
   k <- 0 ## total term counter
   for (j in 1:length(pterms)) {
     tlj <- attr(pterms[[j]],"term.labels") 
     nt <- length(tlj)
     if (j>1 && nt>0) tlj <- paste(tlj,j-1,sep=".")
     term.labels <- c(term.labels,tlj)
     if (nt>0) { # individual parametric terms
       np <- length(object$assign[[j]])
       ind <- pstart[j] - 1 + 1:np 
       Vb <- covmat[ind,ind,drop=FALSE]
       bp <- array(object$coefficients[ind],np)
    
       for (i in 1:nt) { 
         k <- k + 1
         ind <- object$assign[[j]]==i
         b <- bp[ind];V <- Vb[ind,ind]
         ## pseudo-inverse needed in case of truncation of parametric space 
         if (length(b)==1) { 
           V <- 1/V 
           pTerms.df[k] <- nb <- 1      
           pTerms.chi.sq[k] <- V*b*b
         } else {
           V <- pinv(V,length(b),rank.tol=.Machine$double.eps^.5)
           pTerms.df[k] <- nb <- attr(V,"rank")      
           pTerms.chi.sq[k] <- t(b)%*%V%*%b
         }
         if (!est.disp)
         pTerms.pv[k] <- pchisq(pTerms.chi.sq[k],df=nb,lower.tail=FALSE)
         else
         pTerms.pv[k] <- pf(pTerms.chi.sq[k]/nb,df1=nb,df2=residual.df,lower.tail=FALSE)      
       } ## for (i in 1:nt)
     } ## if (nt>0)
   }

   if (npt) {
     attr(pTerms.pv,"names") <- term.labels
     if (!est.disp) {      
       pTerms.table <- cbind(pTerms.df, pTerms.chi.sq, pTerms.pv)   
        dimnames(pTerms.table) <- list(term.labels, c("df", "Chi.sq", "p-value"))
     } else {
       pTerms.table <- cbind(pTerms.df, pTerms.chi.sq/pTerms.df, pTerms.pv)   
        dimnames(pTerms.table) <- list(term.labels, c("df", "F", "p-value"))
     }
   } else { pTerms.df<-pTerms.chi.sq<-pTerms.pv<-array(0,0)}    
    ## ================================

    ## Now deal with the smooth terms....

    m <- length(object$smooth)  # number of smooth terms
    df <- edf1 <-edf <- s.pv <- chi.sq <- array(0, m)
    if (m > 0) { # form test statistics for each smooth
      if (!freq) { ## Bayesian p-values required 
        sub.samp <- max(1000,2*length(object$coefficients)) 
        if (nrow(object$model)>sub.samp) { ## subsample to get X for p-values calc.
          seed <- try(get(".Random.seed",envir=.GlobalEnv),silent=TRUE) ## store RNG seed
          if (inherits(seed,"try-error")) {
            runif(1)
            seed <- get(".Random.seed",envir=.GlobalEnv)
          }
          kind <- RNGkind(NULL)
          RNGkind("default","default")
          set.seed(11) ## ensure repeatability
          ind <- sample(1:nrow(object$model),sub.samp,replace=FALSE)  ## sample these rows from X
          X <- predict(object,object$model[ind,],type="lpmatrix")
          RNGkind(kind[1],kind[2])
          assign(".Random.seed",seed,envir=.GlobalEnv) ## RNG behaves as if it had not been used
        } else { ## don't need to subsample 
          X <- model.matrix(object)
          }
        X <- X[!is.na(rowSums(X)),] ## exclude NA's (possible under na.exclude)
    
       } ## end if (!freq)

     for (i in 1:m) { ## loop through smooths
        start <- object$smooth[[i]]$first.para
        stop <- object$smooth[[i]]$last.para
            
        if (freq) { ## use frequentist cov matrix 
          V <- object$Ve.t[start:stop,start:stop,drop=FALSE] 
        } else V <- object$Vp.t[start:stop,start:stop,drop=FALSE] ## Bayesian
      
        p <- object$coefficients.t[start:stop] # transposed parameters of a smooth
            
        edf1[i] <- edf[i] <- sum(object$edf[start:stop]) # edf for this smooth
        ## extract alternative edf estimate for this smooth, if possible...
        ## edf1 is not done for scam output value...
        if (!is.null(object$edf1)) edf1[i] <-  sum(object$edf1[start:stop])   
       
        if (freq) {
           M1 <- object$smooth[[i]]$df
           M <- min(M1, ceiling(2 * sum(object$edf[start:stop]))) ## upper limit of 2*edf on rank
           V <- pinv(V, M)  # get rank M pseudoinverse of V
           chi.sq[i] <- t(p) %*% V %*% p
           df[i] <- attr(V, "rank")
        }  else { ## Better founded alternatives...
           Xt <- X[, start:stop,drop=FALSE]
           if (object$smooth[[i]]$null.space.dim==0&&!is.null(object$R)) { ## random effect or fully penalized term
             res <- reTest.scam(object,i)
           } else { ## Inverted Nychka interval statistics
             ## df[i] <- min(ncol(Xt),edf1[i])
             if (est.disp) rdf <- residual.df else rdf <- -1
             res <- testStat(p,Xt,V,min(ncol(Xt),edf1[i]),type=0,res.df = rdf) ## was type=p.type
           }
           df[i] <- res$rank
           chi.sq[i] <- res$stat
           s.pv[i] <- res$pval 
        }
        names(chi.sq)[i]<- object$smooth[[i]]$label
  
        if (freq) {
          if (!est.disp)
            s.pv[i] <- pchisq(chi.sq[i], df = df[i], lower.tail = FALSE)
            else
              s.pv[i] <- pf(chi.sq[i]/df[i], df1 = df[i], df2 = residual.df, lower.tail = FALSE)
         #     ## p-values are meaningless for very small edf. Need to set to NA
         #     if (df[i] < 0.1) s.pv[i] <- NA
        }
        ## p-values are meaningless for very small edf. Need to set to NA  
        if (df[i] < 0.1) {
               s.pv[i] <- NA 
               chi.sq[i] <- NA
           }
      }

      ## rounding output values of edf, df and chi.sq...
      edf <- round(edf,digits=4)
      df <- round(df,digits=4)
      chi.sq <- round(chi.sq,digits=4)    

      if (!est.disp) {
        if (freq) {
          s.table <- cbind(edf, df, chi.sq, s.pv)      
          dimnames(s.table) <- list(names(chi.sq), c("edf", "Est.rank", "Chi.sq", "p-value"))
        } else {
            s.table <- cbind(edf, df, chi.sq, s.pv)      
            dimnames(s.table) <- list(names(chi.sq), c("edf", "Ref.df", "Chi.sq", "p-value"))
          }
     } else {
       if (freq) {
           s.table <- cbind(edf, df, chi.sq/df, s.pv)      
           dimnames(s.table) <- list(names(chi.sq), c("edf", "Est.rank", "F", "p-value"))
        } else {
           s.table <- cbind(edf, df, chi.sq/df, s.pv)      
           dimnames(s.table) <- list(names(chi.sq), c("edf", "Ref.df", "F", "p-value"))
          }
       }
   }
   w <- as.numeric(object$prior.weights)
   mean.y <- sum(w*object$y)/sum(w)
   w <- sqrt(w)
   nobs <- nrow(object$model)
   r.sq<- 1 - var(w*(as.numeric(object$y)-object$fitted.values))*(nobs-1)/(var(w*(as.numeric(object$y)-mean.y))*residual.df)
   dev.expl<-(object$null.deviance-object$deviance)/object$null.deviance
   ret<-list(p.coeff=p.coeff,se=se,p.t=p.t,p.pv=p.pv,residual.df=residual.df,m=m,chi.sq=chi.sq,
       s.pv=s.pv,scale=dispersion,r.sq=round(r.sq,digits=4),family=object$family,formula=object$formula,n=nobs,
       dev.expl=dev.expl,edf=edf,dispersion=dispersion,pTerms.pv=pTerms.pv,pTerms.chi.sq=pTerms.chi.sq,
       pTerms.df = pTerms.df, cov.unscaled = covmat.unscaled, cov.scaled = covmat, p.table = p.table,
       pTerms.table = pTerms.table, s.table = s.table,method=object$method,sp.criterion=object$gcv.ubre,
       sp=object$sp,dgcv.ubre=object$dgcv.ubre,termcode=object$termcode,
       gcv.ubre=object$gcv.ubre,optimizer=object$optimizer,rank=object$rank,np=length(object$coefficients))

    class(ret) <- "summary.scam"
    ret
}  ## end summary.scam

   
################## mgcv::: pinvXVX
## (c) Simon N. Wood 

pinvXVX <- function (X, V, rank = NULL) 
{
    k <- floor(rank)
    nu <- rank - k
    if (nu > 0) 
        k1 <- k + 1
    else k1 <- k
    qrx <- qr(X)
    R <- qr.R(qrx)
    V <- R %*% V[qrx$pivot, qrx$pivot] %*% t(R)
    V <- (V + t(V))/2
    ed <- eigen(V, symmetric = TRUE)
    vec <- qr.qy(qrx, rbind(ed$vectors, matrix(0, nrow(X) - ncol(X), 
        ncol(X))))
    if (k1 < ncol(vec)) 
        vec <- vec[, 1:k1, drop = FALSE]
    if (k == 0) {
        vec <- t(t(vec) * sqrt(nu/ed$val[1]))
        return(vec)
    }
    if (nu > 0) {
        if (k > 1) 
            vec[, 1:(k - 1)] <- t(t(vec[, 1:(k - 1)])/sqrt(ed$val[1:(k - 
                1)]))
        b12 <- 0.5 * nu * (1 - nu)
        if (b12 < 0) 
            b12 <- 0
        b12 <- sqrt(b12)
        B <- matrix(c(1, b12, b12, nu), 2, 2)
        ev <- diag(ed$values[k:k1]^-0.5)
        B <- ev %*% B %*% ev
        eb <- eigen(B, symmetric = TRUE)
        rB <- eb$vectors %*% diag(sqrt(eb$values)) %*% t(eb$vectors)
        vec[, k:k1] <- t(rB %*% t(vec[, k:k1]))
    }
    else {
        vec <- t(t(vec)/sqrt(ed$val[1:k]))
    }
    vec
}
 
################ mgcv::: eigXVX
## (c) Simon N. Wood 

eigXVX <- function (X, V, rank = NULL, tol = .Machine$double.eps^0.5) 
{
    qrx <- qr(X)
    R <- qr.R(qrx)
    V <- R %*% V[qrx$pivot, qrx$pivot] %*% t(R)
    V <- (V + t(V))/2
    ed <- eigen(V, symmetric = TRUE)
    ind <- abs(ed$values) > max(abs(ed$values)) * tol
    erank <- sum(ind)
    if (is.null(rank)) {
        rank <- erank
    }
    else {
        if (rank < erank) 
            ind <- 1:rank
        else rank <- erank
    }
    vec <- qr.qy(qrx, rbind(ed$vectors, matrix(0, nrow(X) - ncol(X), 
        ncol(X))))
    list(values = ed$values[ind], vectors = vec[, ind], rank = rank)
}


##### print.summary.scam .....
print.summary.scam <- function (x, digits = max(3, getOption("digits") - 3), 
         signif.stars = getOption("show.signif.stars"),...) 
## print method for scam, a clone of print.summary.gam of mgcv()
{
    print(x$family)
    cat("Formula:\n")
    print(x$formula)
    if (length(x$p.coeff) > 0) {
        cat("\nParametric coefficients:\n")
        printCoefmat(x$p.table, digits = digits, signif.stars = signif.stars, 
            na.print = "NA", ...)
    }
    cat("\n")
    if (x$m > 0) {
        cat("Approximate significance of smooth terms:\n")
        printCoefmat(x$s.table, digits = digits, signif.stars = signif.stars, 
            has.Pvalue = TRUE, na.print = "NA", cs.ind = 1, ...)
    }
   # cat("\n")
   if (!is.null(x$rank) && x$rank< x$np) cat("Rank: ",x$rank,"/",x$np,"\n",sep="") 
    cat("\nR-sq.(adj) = ", formatC(x$r.sq, digits = 4, width = 5))
    if (length(x$dev.expl) > 0) 
        cat("   Deviance explained = ", formatC(x$dev.expl * 
            100, digits = 3, width = 4), "%\n", sep = "")
    if (length(x$sp.criterion) > 0) 
        cat( x$method," score = ", formatC(x$sp.criterion, digits = 5), 
            sep = "")
    cat("  Scale est. = ", formatC(x$scale, digits = 5, width = 8, 
        flag = "-"), "  n = ", x$n, "\n", sep = "")
    if ((x$optimizer[1] == "bfgs") && x$m>0){
               if (x$termcode!= 1) {
                   dgcv.ubre <- max(abs(x$dgcv.ubre)*max(abs(log(x$sp)),1)/max(abs(x$gcv.ubre),1))
                  cat("\nBFGS termination condition:\n", dgcv.ubre,"\n",sep = "")
               }   
    }
    cat("\n")
    invisible(x)
}




###############################################
## anova for scam models (clone of summary.gam())...
## of mgcv versions up to 1.8-11...
## (c) Simon N. Wood 
###############################################


anova.scam <- function (object, ..., dispersion = NULL, test = NULL,  freq=FALSE,p.type=0)
# clone of summary.gam(): mgcv package
{   # adapted from anova.glm: R stats package
    dotargs <- list(...)
    named <- if (is.null(names(dotargs)))
        rep(FALSE, length(dotargs))
    else (names(dotargs) != "")
    if (any(named))
        warning("The following arguments to anova.glm(..) are invalid and dropped: ",
            paste(deparse(dotargs[named]), collapse = ", "))
    dotargs <- dotargs[!named]
    is.glm <- unlist(lapply(dotargs, function(x) inherits(x,
        "glm")))
    dotargs <- dotargs[is.glm]
    if (length(dotargs) > 0)
     return(anova(structure(c(list(object), dotargs), class="glmlist"), 
            dispersion = dispersion, test = test)) 
       # return(anova.glmlist(c(list(object), dotargs), dispersion = dispersion,
       #     test = test)) ## modified at BDR's suggestion 19/08/13
    if (!is.null(test)) warning("test argument ignored")
    if (!inherits(object,"scam")) stop("anova.scam called with non scam object")
    sg <- summary(object, dispersion = dispersion, freq = freq,p.type=p.type)
    class(sg) <- "anova.scam"
    sg
} ## anova.scam


print.anova.scam <- function(x, digits = max(3, getOption("digits") - 3), ...)
{ # print method for class anova.scam resulting from single
  # scam model calls to anova. Clone of print.anova.gam(): mgcv package
  print(x$family)
  cat("Formula:\n") 
  if (is.list(x$formula)) for (i in 1:length(x$formula)) print(x$formula[[i]]) else
  print(x$formula)
  if (length(x$pTerms.pv)>0)
  { cat("\nParametric Terms:\n")
    printCoefmat(x$pTerms.table, digits = digits, signif.stars = FALSE, has.Pvalue = TRUE, na.print = "NA", ...)
  }
  cat("\n")
  if(x$m>0)
  { cat("Approximate significance of smooth terms:\n")
    printCoefmat(x$s.table, digits = digits, signif.stars = FALSE, has.Pvalue = TRUE, na.print = "NA", ...)
  }
  invisible(x)
} ## print.anova.scam