R/FUN_spatial_a.R

In sommer: Solving Mixed Model Equations in R

spl2Da <-  function(x.coord,y.coord,at.var=NULL,at.levels=NULL, type="PSANOVA", 
                    nsegments = c(10,10), penaltyord = c(2,2), degree = c(3,3), 
                    nestorder = c(1,1) ) {
  
  ##
  if(length(degree) == 1){degree <- rep(degree,2) }
  if(length(nsegments) == 1){nsegments <- rep(nsegments,2) }
  if(length(penaltyord) == 1){penaltyord <- rep(penaltyord,2) }
  if(length(nestorder) == 1){nestorder <- rep(nestorder,2) }
  
  x.coord.name <- as.character(substitute(list(x.coord)))[-1L]
  y.coord.name <- as.character(substitute(list(y.coord)))[-1L]
  
  # x.coord.name <- "col"
  # y.coord.name <- "row"
  
  if(!is.numeric(x.coord)){
    stop("x.coord argument in spl2D() needs to be numeric.", call. = FALSE)
  }
  if(!is.numeric(y.coord)){
    stop("y.coord argument in spl2D() needs to be numeric.", call. = FALSE)
  }
  
  interpret.covarrubias.formula <-
    function(formula) {
      env <- environment(formula) 
      if(inherits(formula, "character"))          
        formula <- as.formula(formula)
      tf <- terms.formula(formula, specials = c("SAP", "PSANOVA"))
      terms <- attr(tf, "term.labels")
      nt <- length(terms)
      if(nt != 1)
        stop("Error in the specification of the spatial effect: only a sigle bidimensional function is allowed")
      
      res <- eval(parse(text = terms[1]), envir = env)
      res
    }
  
  bbase <- function(X., XL., XR., NDX., BDEG.) {
    # Function for B-spline basis
    dx <- (XR. - XL.)/NDX.
    knots <- seq(XL. - BDEG.*dx, XR. + BDEG.*dx, by=dx)
    P <- outer(X., knots, tpower, BDEG.)
    n <- dim(P)[2]
    D <- diff(diag(n), diff = BDEG. + 1) / (gamma(BDEG. + 1) * dx ^ BDEG.)
    B <- (-1) ^ (BDEG. + 1) * P %*% t(D)
    res <- list(B = B, knots = knots)
    res 
  }
  
  tpower <-
    function(x, t, p) {
      # Function for truncated p-th power function
      return((x - t) ^ p * (x > t))
    }
  
  Rten2 <-
    function(X1,X2) {
      one.1 <- matrix(1,1,ncol(X1))
      one.2 <- matrix(1,1,ncol(X2))
      kronecker(X1,one.2)*kronecker(one.1,X2)
    }
  
  MM.basis <-
    function (x, xl, xr, ndx, bdeg, penaltyord, decom = 1) {
      Bb = bbase(x,xl,xr,ndx,bdeg)
      knots <- Bb$knots
      B = Bb$B
      m = ncol(B)
      n = nrow(B)
      D = diff(diag(m), differences=penaltyord)
      P.svd = svd(crossprod(D))
      U.Z = (P.svd$u)[,1:(m-penaltyord)] # eigenvectors
      d = (P.svd$d)[1:(m-penaltyord)]  # eigenvalues
      Z = B%*%U.Z
      U.X = NULL
      if(decom == 1) {
        U.X = ((P.svd$u)[,-(1:(m-penaltyord))])
        X = B%*%U.X
      } else if (decom == 2){
        X = NULL
        for(i in 0:(penaltyord-1)){
          X = cbind(X,x^i)
        }
      } else if(decom == 3) {
        U.X = NULL
        for(i in 0:(penaltyord-1)){
          U.X = cbind(U.X,knots[-c((1:penaltyord),(length(knots)- penaltyord + 1):length(knots))]^i)
        }
        X = B%*%U.X
      } else if(decom == 4) { # Wood's 2013
        X = B%*%((P.svd$u)[,-(1:(m-penaltyord))])
        id.v <- rep(1, nrow(X))
        D.temp = X - ((id.v%*%t(id.v))%*%X)/nrow(X)
        Xf <- svd(crossprod(D.temp))$u[,ncol(D.temp):1]
        X <- X%*%Xf
        U.X = ((P.svd$u)[,-(1:(m-penaltyord)), drop = FALSE])%*%Xf
      }
      list(X = X, Z = Z, d = d, B = B, m = m, D = D, knots = knots, U.X = U.X, U.Z = U.Z)
    }
  
  ####################
  ### if we want to use at.var and at.levels
  if(is.null(at.var)){
    at.var <- rep("A",length(x.coord))
    at.name <- "FIELDINST"
    at.levels <- "A"
  }else{
    at.name <- as.character(substitute(list(at)))[-1L]
    if(length(at.var) != length(x.coord)){stop("at.var has different length than x.coord and y.coord, please fix.", call. = FALSE)}
    if(is.null(at.levels)){at.levels <- levels(as.factor(at.var))}
  }
  #######################
  index <- 1:length(x.coord)
  if(!is.numeric(x.coord)){stop("x.coord argument in spl2D() needs to be numeric.", call. = FALSE)}
  if(!is.numeric(y.coord)){stop("y.coord argument in spl2D() needs to be numeric.", call. = FALSE)}
  #######################
  ## split data by the "at.name" argument
  dat <- data.frame(x.coord, y.coord, at.var, index); colnames(dat) <- c(x.coord.name,y.coord.name,at.name,"index")
  missby <- which(is.na(dat[,at.name]))
  if(length(missby)>0){stop("We will split using the at.name argument and you have missing values in this column.\nPlease correct.", call. = FALSE)}
  dat[,at.name] <- as.factor(dat[,at.name])
  data0 <- split(dat, dat[,at.name])
  names(data0) <- levels(dat[,at.name])
  #######################################
  # make sure there's no missing data in coordinate variables
  nasx <- which(is.na(dat[,x.coord.name]))
  nasy <- which(is.na(dat[,y.coord.name]))
  if(length(nasx) > 0 | length(nasy) >0){
    stop("x.coord and y.coord columns cannot have NA's", call. = FALSE)
  }
  #######################################
  ## now apply the same to all environments
  multires <- lapply(data0, function(dxy){
    
    
    x1 <- dxy[ ,x.coord.name]
    x2 <- dxy[ ,y.coord.name]
    # print(x2)
    #type = type
    
    MM1 = MM.basis(x1, min(x1), max(x1), nsegments[1], degree[1], penaltyord[1], 4)
    MM2 = MM.basis(x2, min(x2), max(x2), nsegments[2], degree[2], penaltyord[2], 4)
    
    X1 <- MM1$X; Z1 <- MM1$Z; d1 <- MM1$d; B1 <- MM1$B
    X2 <- MM2$X; Z2 <- MM2$Z; d2 <- MM2$d; B2 <- MM2$B
    
    c1 = ncol(B1); c2 = ncol(B2)
    
    # Nested bases
    if(nestorder[1] == 1) {
      MM1n <- MM1
      Z1n <- Z1
      c1n <- c1
      d1n <- d1	
    } else {
      MM1n = MM.basis(x1, min(x1), max(x1), nsegments[1]/nestorder[1], degree[1], penaltyord[1], 4)
      Z1n <- MM1n$Z
      d1n <- MM1n$d
      c1n <-  ncol(MM1n$B)  					
    }
    if(nestorder[2] == 1) {
      MM2n <- MM2
      Z2n <- Z2
      c2n <- c2
      d2n <- d2	
    } else {
      MM2n = MM.basis(x2, min(x2), max(x2), nsegments[2]/nestorder[2], degree[2], penaltyord[2], 4)
      Z2n <- MM2n$Z
      d2n <- MM2n$d
      c2n <-  ncol(MM2n$B)  					
    }
    
    x.fixed <- y.fixed <- ""
    for(i in 0:(penaltyord[1]-1)){
      if(i == 1) 
        x.fixed <- c(x.fixed, x.coord.name)
      else if( i > 1)
        x.fixed <- c(x.fixed, paste(x.coord.name, "^", i, sep = ""))
    }
    for(i in 0:(penaltyord[2]-1)){
      if(i == 1) 
        y.fixed <- c(y.fixed, y.coord.name)
      else if( i > 1)
        y.fixed <- c(y.fixed, paste(y.coord.name, "^", i, sep = ""))
    }
    xy.fixed <- NULL
    for(i in 1:length(y.fixed)) {
      xy.fixed <- c(xy.fixed, paste(y.fixed[i], x.fixed, sep= ""))
    }
    xy.fixed <- xy.fixed[xy.fixed != ""]
    names.fixed <- xy.fixed
    
    smooth.comp <- paste("f.", x.coord.name,".", y.coord.name,"", sep = "")
    
    if(type == "SAP") {
      names.random <- paste(smooth.comp, c(x.coord.name, y.coord.name), sep = "|")				
      X = Rten2(X2, X1)		
      # Delete the intercept
      X <- X[,-1,drop = FALSE]
      Z = cbind(Rten2(X2, Z1), Rten2(Z2, X1), Rten2(Z2n, Z1n))
      
      dim.random <- c((c1 -penaltyord[1])*penaltyord[2] , (c2 - penaltyord[2])*penaltyord[1], (c1n - penaltyord[1])*(c2n - penaltyord[2]))		
      dim <- list(fixed = rep(1, ncol(X)), random = sum(dim.random))
      names(dim$fixed) <- names.fixed
      names(dim$random) <- paste(smooth.comp, "Global")
      
      # Variance/Covariance components
      g1u <- rep(1, penaltyord[2])%x%d1
      g2u <- d2%x%rep(1, penaltyord[1])
      g1b <- rep(1, c2n - penaltyord[2])%x%d1n
      g2b <- d2n%x%rep(1, c1n - penaltyord[1])
      
      g <- list()	
      g[[1]] <- c(g1u, rep(0, dim.random[2]), g1b)
      g[[2]] <- c(rep(0, dim.random[1]), g2u, g2b)
      
      names(g) <- names.random
      
    } else {		
      one1. <- X1[,1, drop = FALSE]
      one2. <- X2[,1, drop = FALSE]
      
      x1. <- X1[,-1, drop = FALSE]
      x2. <- X2[,-1, drop = FALSE]
      
      # Fixed and random matrices
      X <- Rten2(X2, X1)
      # Delete the intercept
      X <- X[,-1,drop = FALSE]
      Z <- cbind(Rten2(one2., Z1), Rten2(Z2, one1.), Rten2(x2., Z1), Rten2(Z2, x1.), Rten2(Z2n, Z1n))
      
      dim.random <- c((c1-penaltyord[1]), (c2-penaltyord[2]), (c1-penaltyord[1])*(penaltyord[2]-1), (c2-penaltyord[2])*(penaltyord[1]-1), (c1n-penaltyord[2])*(c2n-penaltyord[2]))
      
      # Variance/Covariance components		
      g1u <- d1
      g2u <- d2
      
      g1v <- rep(1, penaltyord[2] - 1)%x%d1
      g2v <- d2%x%rep(1,penaltyord[1] - 1)
      
      g1b <- rep(1, c2n - penaltyord[2])%x%d1n
      g2b <- d2n%x%rep(1, c1n - penaltyord[1])
      
      g <- list()
      
      if(type == "SAP.ANOVA") {
        g[[1]] <- c(g1u, rep(0, sum(dim.random[2:5])))
        g[[2]] <- c(rep(0, dim.random[1]), g2u, rep(0, sum(dim.random[3:5])))
        g[[3]] <- c(rep(0, sum(dim.random[1:2])), g1v, rep(0, dim.random[4]), g1b)
        g[[4]] <- c(rep(0, sum(dim.random[1:3])), g2v, g2b)
        
        names.random <- c(paste("f.", x.coord.name,"", sep = ""), paste("f.", y.coord.name,"", sep = ""), paste(smooth.comp, c(x.coord.name, y.coord.name), sep = "."))			
        dim <- list(fixed = rep(1, ncol(X)), random = c(dim.random[1:2], sum(dim.random[-(1:2)])))		
        names(dim$fixed) <- names.fixed
        names(dim$random) <- c(names.random[1:2], paste(smooth.comp, "Global"))
        names(g) <- names.random
      } else {
        g[[1]] <- c(g1u, rep(0, sum(dim.random[2:5])))
        g[[2]] <- c(rep(0, dim.random[1]), g2u, rep(0, sum(dim.random[3:5])))
        g[[3]] <- c(rep(0, sum(dim.random[1:2])), g1v, rep(0, sum(dim.random[4:5])))
        g[[4]] <- c(rep(0, sum(dim.random[1:3])), g2v, rep(0, dim.random[5]))
        g[[5]] <- c(rep(0, sum(dim.random[1:4])), g1b + g2b)
        
        names.random <- c(paste("f.", x.coord.name,"", sep = ""), paste("f.", y.coord.name,"", sep = ""),
                          paste("f.", x.coord.name,".", y.coord.name, sep = ""),
                          paste(x.coord.name,".f.", y.coord.name,"", sep = ""),
                          paste("f.", x.coord.name,".f.", y.coord.name,"", sep = ""))
        # print(names.random)
        
        dim <- list(fixed = rep(1, ncol(X)), random = dim.random)		
        names(dim$fixed) <- names.fixed
        names(dim$random) <- names.random
        names(g) <- names.random
      }		
    }
    colnames(X) <- names.fixed
    colnames(Z) <- paste(smooth.comp, 1:ncol(Z), sep = ".")
    
    attr(dim$fixed, "random") <- attr(dim$fixed, "sparse") <- rep(FALSE, length(dim$fixed))
    attr(dim$fixed, "spatial") <- rep(TRUE, length(dim$fixed))
    
    attr(dim$random, "random") <- attr(dim$random, "spatial") <- rep(TRUE, length(dim$random)) 
    attr(dim$random, "sparse") <- rep(FALSE, length(dim$random))
    
    terms <- list()
    terms$MM <- list(MM1 = MM1, MM2 = MM2)
    terms$MMn <- list(MM1 = MM1n, MM2 = MM2n)
    #terms$terms.formula <- res
    
    # attr(terms, "term") <- smooth.comp
    
    # Initialize variance components
    init.var <- rep(1, length(g))
    
    res <- list(X = X, Z = Z, dim = dim, g = g, init.var = init.var)	
    M <- list(cbind(res$X,res$Z))
    
    return(M)
  })
  
  names(multires) <- at.levels
  # print(str(multires))
  # print(lapply(multires,colnames))
  toFill <- lapply(data0,function(x){x$index})
  #############################################
  ## CAPTURE COLNAMES AND BUILD A MATRIX WITH THOSE NAMES
  myNames <- list()
  for(k in 1:1){
    # print(lapply(multires,function(x){colnames(x[[k]])}))
    myNames[[k]] <- lapply(multires,function(x){colnames(x[[k]])})
  };  names(myNames) <- c("all")
  # print(myNames)
  #################################
  ## move matrices to the right size
  Zup <- list() # store incidence matrices
  Zup2 <- list() # store incidence matrices
  Kup <- list() # store relationship matrices between levels in Z
  typevc <- numeric() # store wheter is a variance (1) or covariance (2;allowed to be negative) component
  re_name <- character() # store the name of the random effect
  counter <- 1
  counter2 <- 1
  for(k in 1:1){ # for each tensor product (single matrix in this case)
    for(j in 1:length(multires)){ # for each environment or by.level
      if(names(multires)[j] %in% at.levels){
        # print("yes")
        Z <- matrix(0,nrow=nrow(dat),ncol=length(myNames[[k]][[j]]))
        colnames(Z) <- myNames[[k]][[j]]
        # print(dim(Z))
        
        prov <- multires[[j]][[k]]
        # print(dim(prov))
        Z[toFill[[j]],colnames(prov)] <- as.matrix(prov) 
        attr(Z,"variables") <- c(x.coord.name, y.coord.name)
        Zup[[counter]] <- Z
        names(Zup)[counter] <- paste0(names(multires)[j],":",names(myNames)[k])
        Gu1 <- diag(ncol(Z)); colnames(Gu1) <- rownames(Gu1) <- colnames(Z)
        Kup[[counter]] <- Gu1
        typevc[counter] <- 1
        re_name[counter] <- names(Zup)[counter]
        counter <- counter + 1
      } # else don't fill that portion of the matrix
    }
  }
  Gti=NULL
  Gtc=NULL
  vcs <- diag(length(Zup)); rownames(vcs) <- colnames(vcs) <- names(Zup)
  #################################
  namess2 <- c(x.coord.name, y.coord.name)
  S3 <- list(Z=Zup,K=Kup,Gti=Gti,Gtc=Gtc,typevc=typevc,re_name=re_name,vcs=vcs, terms=namess2)
  return(S3)
}


nna<- function (pheno, trait = "y", rown = "row", coln = "col", nrows = 1, 
                ncols = 2) 
{
  bases <- c(rown, coln, trait)
  into <- intersect(colnames(pheno), bases)
  if (length(into) < 3) {
    stop("Please provide the arguments 'rown', 'coln' and 'trait'.", 
         call. = FALSE)
  }
  newd <- pheno[, c(rown, coln, trait)]
  rn <- rownames(newd)
  newd <- apply(newd, 2, as.numeric)
  rownames(newd) <- rn
  newd <- data.frame(newd)
  pheno[, c(rown, coln, trait)] <- newd
  pheno$nnx <- NA
  nnx <- apply(newd, 1, function(x) {
    s1 <- as.numeric(x[1]-nrows)
    s2 <- as.numeric(x[1] - 1)
    s3 <- as.numeric(x[1] + 1)
    s4 <- as.numeric(x[1] + nrows)
    ns <- c(s1:s2, s3:s4)
    
    s5 <- as.numeric(x[2] - ncols)
    s6 <- as.numeric(x[2] - 1)
    s7 <- as.numeric(x[2] + 1)
    s8 <- as.numeric(x[2] + ncols)
    we <- c(s5:s6, s7:s8)
    
    
    ns <- ns[which(ns > 0)]
    we <- we[which(we > 0)]
    v1 <- which((newd[, rown] %in% ns) & (newd[, coln] %in% x[2]))
    v2 <- which(newd[, coln] %in% we & newd[, rown] %in% x[1])
    if (length(v1) > 0 & length(v2) > 0) {
      yy <- x[3] - apply(newd[c(v1, v2), ], 2, mean, na.rm = TRUE)[3]
    }  else if (length(v1) > 0 & length(v2) == 0) {
      yy <- x[3] - apply(newd[c(v1), ], 2, mean, na.rm = TRUE)[3]
    }  else if (length(v1) > 0 & length(v2) == 0) {
      yy <- x[3] - apply(newd[c(v2), ], 2, mean, na.rm = TRUE)[3]
    } else{
      yy<-x[3]
    }
    return(yy)
  })
  pheno$nnx <- as.numeric(unlist(nnx))
  return(pheno)
}

# spl2D <-  function(x.coord,y.coord,at,at.levels, type="PSANOVA", nsegments = c(10,10), penaltyord = c(2,2), degree = c(3,3), nestorder = c(1,1) ) {
#   
#   ##
#   init0 <- as.character(substitute(list(x.coord)))[-1L]
#   init1 <- as.character(substitute(list(y.coord)))[-1L]
#   
#   x.coord.name <- x.coord
#   y.coord.name <- y.coord
#   
#   if(!is.numeric(x.coord)){
#     stop("x.coord argument in spl2D() needs to be numeric.", call. = FALSE)
#   }
#   if(!is.numeric(y.coord)){
#     stop("y.coord argument in spl2D() needs to be numeric.", call. = FALSE)
#   }
#   
#   interpret.covarrubias.formula <-
#     function(formula) {
#       env <- environment(formula) 
#       if(inherits(formula, "character"))          
#         formula <- as.formula(formula)
#       tf <- terms.formula(formula, specials = c("SAP", "PSANOVA"))
#       terms <- attr(tf, "term.labels")
#       nt <- length(terms)
#       if(nt != 1)
#         stop("Error in the specification of the spatial effect: only a sigle bidimensional function is allowed")
#       
#       res <- eval(parse(text = terms[1]), envir = env)
#       res
#     }
#   
#   bbase <-
#     function(X., XL., XR., NDX., BDEG.) {
#       # Function for B-spline basis
#       dx <- (XR. - XL.)/NDX.
#       knots <- seq(XL. - BDEG.*dx, XR. + BDEG.*dx, by=dx)
#       P <- outer(X., knots, tpower, BDEG.)
#       n <- dim(P)[2]
#       D <- diff(diag(n), diff = BDEG. + 1) / (gamma(BDEG. + 1) * dx ^ BDEG.)
#       B <- (-1) ^ (BDEG. + 1) * P %*% t(D)
#       res <- list(B = B, knots = knots)
#       res 
#     }
#   
#   tpower <-
#     function(x, t, p) {
#       # Function for truncated p-th power function
#       return((x - t) ^ p * (x > t))
#     }
#   
#   Rten2 <-
#     function(X1,X2) {
#       one.1 <- matrix(1,1,ncol(X1))
#       one.2 <- matrix(1,1,ncol(X2))
#       kronecker(X1,one.2)*kronecker(one.1,X2)
#     }
#   
#   MM.basis <-
#     function (x, xl, xr, ndx, bdeg, penaltyord, decom = 1) {
#       Bb = bbase(x,xl,xr,ndx,bdeg)
#       knots <- Bb$knots
#       B = Bb$B
#       m = ncol(B)
#       n = nrow(B)
#       D = diff(diag(m), differences=penaltyord)
#       P.svd = svd(crossprod(D))
#       U.Z = (P.svd$u)[,1:(m-penaltyord)] # eigenvectors
#       d = (P.svd$d)[1:(m-penaltyord)]  # eigenvalues
#       Z = B%*%U.Z
#       U.X = NULL
#       if(decom == 1) {
#         U.X = ((P.svd$u)[,-(1:(m-penaltyord))])
#         X = B%*%U.X
#       } else if (decom == 2){
#         X = NULL
#         for(i in 0:(penaltyord-1)){
#           X = cbind(X,x^i)
#         }
#       } else if(decom == 3) {
#         U.X = NULL
#         for(i in 0:(penaltyord-1)){
#           U.X = cbind(U.X,knots[-c((1:penaltyord),(length(knots)- penaltyord + 1):length(knots))]^i)
#         }
#         X = B%*%U.X
#       } else if(decom == 4) { # Wood's 2013
#         X = B%*%((P.svd$u)[,-(1:(m-penaltyord))])
#         id.v <- rep(1, nrow(X))
#         D.temp = X - ((id.v%*%t(id.v))%*%X)/nrow(X)
#         Xf <- svd(crossprod(D.temp))$u[,ncol(D.temp):1]
#         X <- X%*%Xf
#         U.X = ((P.svd$u)[,-(1:(m-penaltyord)), drop = FALSE])%*%Xf
#       }
#       list(X = X, Z = Z, d = d, B = B, m = m, D = D, knots = knots, U.X = U.X, U.Z = U.Z)
#     }
#   
#   ####################
#   ### if we want to use at and at.levels
#   if(!missing(at)){
#     col1 <- deparse(substitute(at))
#     # print(col1)
#     dat <- data.frame(x.coord, y.coord, at); colnames(dat) <- c("x.coord","y.coord",col1)
#     by <- col1
#     if(!missing(at.levels)){
#       by.levels=at.levels
#     }else{by.levels=NULL}
#   }else{
#     by=NULL
#     by.levels=NULL
#     dat <- data.frame(x.coord,y.coord); colnames(dat) <- c("x.coord","y.coord")
#   }
#   #######################
#   
#   x.coord <- "x.coord"
#   y.coord <- "y.coord"
#   
#   if(is.null(by)){
#     dat$FIELDINST <- "FIELD1"
#     by="FIELDINST"
#     dat[,by] <- as.factor(dat[,by])
#     data0 <- split(dat, dat[,by])
#     if(!is.null(by.levels)){
#       keep <- names(data0)[which(names(data0) %in% by.levels)]
#       
#       if(length(keep)==0){stop("The by.levels provided were not found in your dataset.",call. = FALSE)}
#       
#       data0 <- data0[[keep]]
#       if(length(keep)==1){data0 <- list(data0); names(data0) <- keep}
#     }
#   }else{
#     check <- which(colnames(dat)==by)
#     if(length(check)==0){stop("by argument not found in the dat provided", call. = FALSE)}else{
#       
#       missby <- which(is.na(dat[,by]))
#       if(length(missby)>0){stop("We will split using the by argument and you have missing values in this column.\nPlease correct.", call. = FALSE)}
#       
#       dat[,by] <- as.factor(dat[,by])
#       data0 <- split(dat, dat[,by])
#       
#       if(!is.null(by.levels)){
#         keep <- names(data0)[which(names(data0) %in% by.levels)]
#         
#         if(length(keep)==0){stop("The by.levels provided were not found in your dataset.",call. = FALSE)}
#         
#         data0 <- data0[[keep]]
#         if(length(keep)==1){data0 <- list(data0); names(data0) <- keep}
#         #print(str(data0))
#       }
#       
#     }
#   }
#   
#   nasx <- which(is.na(dat[,x.coord]))
#   nasy <- which(is.na(dat[,y.coord]))
#   if(length(nasx) > 0 | length(nasy) >0){
#     stop("x.coord and y.coord columns cannot have NA's", call. = FALSE)
#   }
#   #res <- interpret.covarrubias.formula(formula)
#   
#   ####
#   #### now apply the same to all environments
#   multires <- lapply(data0, function(data){
#     
#     
#     x1 <- data[ ,x.coord]
#     x2 <- data[ ,y.coord]
#     
#     #type = type
#     
#     MM1 = MM.basis(x1, min(x1), max(x1), nsegments[1], degree[1], penaltyord[1], 4)
#     MM2 = MM.basis(x2, min(x2), max(x2), nsegments[2], degree[2], penaltyord[2], 4)
#     
#     X1 <- MM1$X; Z1 <- MM1$Z; d1 <- MM1$d; B1 <- MM1$B
#     X2 <- MM2$X; Z2 <- MM2$Z; d2 <- MM2$d; B2 <- MM2$B
#     
#     c1 = ncol(B1); c2 = ncol(B2)
#     
#     # Nested bases
#     if(nestorder[1] == 1) {
#       MM1n <- MM1
#       Z1n <- Z1
#       c1n <- c1
#       d1n <- d1	
#     } else {
#       MM1n = MM.basis(x1, min(x1), max(x1), nsegments[1]/nestorder[1], degree[1], penaltyord[1], 4)
#       Z1n <- MM1n$Z
#       d1n <- MM1n$d
#       c1n <-  ncol(MM1n$B)  					
#     }
#     if(nestorder[2] == 1) {
#       MM2n <- MM2
#       Z2n <- Z2
#       c2n <- c2
#       d2n <- d2	
#     } else {
#       MM2n = MM.basis(x2, min(x2), max(x2), nsegments[2]/nestorder[2], degree[2], penaltyord[2], 4)
#       Z2n <- MM2n$Z
#       d2n <- MM2n$d
#       c2n <-  ncol(MM2n$B)  					
#     }
#     
#     x.fixed <- y.fixed <- ""
#     for(i in 0:(penaltyord[1]-1)){
#       if(i == 1) 
#         x.fixed <- c(x.fixed, x.coord)
#       else if( i > 1)
#         x.fixed <- c(x.fixed, paste(x.coord, "^", i, sep = ""))
#     }
#     for(i in 0:(penaltyord[2]-1)){
#       if(i == 1) 
#         y.fixed <- c(y.fixed, y.coord)
#       else if( i > 1)
#         y.fixed <- c(y.fixed, paste(y.coord, "^", i, sep = ""))
#     }
#     xy.fixed <- NULL
#     for(i in 1:length(y.fixed)) {
#       xy.fixed <- c(xy.fixed, paste(y.fixed[i], x.fixed, sep= ""))
#     }
#     xy.fixed <- xy.fixed[xy.fixed != ""]
#     names.fixed <- xy.fixed
#     
#     smooth.comp <- paste("f(", x.coord,",", y.coord,")", sep = "")
#     
#     if(type == "SAP") {
#       names.random <- paste(smooth.comp, c(x.coord, y.coord), sep = "|")				
#       X = Rten2(X2, X1)		
#       # Delete the intercept
#       X <- X[,-1,drop = FALSE]
#       Z = cbind(Rten2(X2, Z1), Rten2(Z2, X1), Rten2(Z2n, Z1n))
#       
#       dim.random <- c((c1 -penaltyord[1])*penaltyord[2] , (c2 - penaltyord[2])*penaltyord[1], (c1n - penaltyord[1])*(c2n - penaltyord[2]))		
#       dim <- list(fixed = rep(1, ncol(X)), random = sum(dim.random))
#       names(dim$fixed) <- names.fixed
#       names(dim$random) <- paste(smooth.comp, "Global")
#       
#       # Variance/Covariance components
#       g1u <- rep(1, penaltyord[2])%x%d1
#       g2u <- d2%x%rep(1, penaltyord[1])
#       g1b <- rep(1, c2n - penaltyord[2])%x%d1n
#       g2b <- d2n%x%rep(1, c1n - penaltyord[1])
#       
#       g <- list()	
#       g[[1]] <- c(g1u, rep(0, dim.random[2]), g1b)
#       g[[2]] <- c(rep(0, dim.random[1]), g2u, g2b)
#       
#       names(g) <- names.random
#       
#     } else {		
#       one1. <- X1[,1, drop = FALSE]
#       one2. <- X2[,1, drop = FALSE]
#       
#       x1. <- X1[,-1, drop = FALSE]
#       x2. <- X2[,-1, drop = FALSE]
#       
#       # Fixed and random matrices
#       X <- Rten2(X2, X1)
#       # Delete the intercept
#       X <- X[,-1,drop = FALSE]
#       Z <- cbind(Rten2(one2., Z1), Rten2(Z2, one1.), Rten2(x2., Z1), Rten2(Z2, x1.), Rten2(Z2n, Z1n))
#       
#       dim.random <- c((c1-penaltyord[1]), (c2-penaltyord[2]), (c1-penaltyord[1])*(penaltyord[2]-1), (c2-penaltyord[2])*(penaltyord[1]-1), (c1n-penaltyord[2])*(c2n-penaltyord[2]))
#       
#       # Variance/Covariance components		
#       g1u <- d1
#       g2u <- d2
#       
#       g1v <- rep(1, penaltyord[2] - 1)%x%d1
#       g2v <- d2%x%rep(1,penaltyord[1] - 1)
#       
#       g1b <- rep(1, c2n - penaltyord[2])%x%d1n
#       g2b <- d2n%x%rep(1, c1n - penaltyord[1])
#       
#       g <- list()
#       
#       if(type == "SAP.ANOVA") {
#         g[[1]] <- c(g1u, rep(0, sum(dim.random[2:5])))
#         g[[2]] <- c(rep(0, dim.random[1]), g2u, rep(0, sum(dim.random[3:5])))
#         g[[3]] <- c(rep(0, sum(dim.random[1:2])), g1v, rep(0, dim.random[4]), g1b)
#         g[[4]] <- c(rep(0, sum(dim.random[1:3])), g2v, g2b)
#         
#         names.random <- c(paste("f(", x.coord,")", sep = ""), paste("f(", y.coord,")", sep = ""), paste(smooth.comp, c(x.coord, y.coord), sep = "|"))			
#         dim <- list(fixed = rep(1, ncol(X)), random = c(dim.random[1:2], sum(dim.random[-(1:2)])))		
#         names(dim$fixed) <- names.fixed
#         names(dim$random) <- c(names.random[1:2], paste(smooth.comp, "Global"))
#         names(g) <- names.random
#       } else {
#         g[[1]] <- c(g1u, rep(0, sum(dim.random[2:5])))
#         g[[2]] <- c(rep(0, dim.random[1]), g2u, rep(0, sum(dim.random[3:5])))
#         g[[3]] <- c(rep(0, sum(dim.random[1:2])), g1v, rep(0, sum(dim.random[4:5])))
#         g[[4]] <- c(rep(0, sum(dim.random[1:3])), g2v, rep(0, dim.random[5]))
#         g[[5]] <- c(rep(0, sum(dim.random[1:4])), g1b + g2b)
#         
#         names.random <- c(paste("f(", x.coord,")", sep = ""), paste("f(", y.coord,")", sep = ""),
#                           paste("f(", x.coord,"):", y.coord, sep = ""),
#                           paste(x.coord,":f(", y.coord,")", sep = ""),
#                           paste("f(", x.coord,"):f(", y.coord,")", sep = ""))
#         
#         dim <- list(fixed = rep(1, ncol(X)), random = dim.random)		
#         names(dim$fixed) <- names.fixed
#         names(dim$random) <- names.random
#         names(g) <- names.random
#       }		
#     }
#     colnames(X) <- names.fixed
#     colnames(Z) <- paste(smooth.comp, 1:ncol(Z), sep = ".")
#     
#     attr(dim$fixed, "random") <- attr(dim$fixed, "sparse") <- rep(FALSE, length(dim$fixed))
#     attr(dim$fixed, "spatial") <- rep(TRUE, length(dim$fixed))
#     
#     attr(dim$random, "random") <- attr(dim$random, "spatial") <- rep(TRUE, length(dim$random)) 
#     attr(dim$random, "sparse") <- rep(FALSE, length(dim$random))
#     
#     terms <- list()
#     terms$MM <- list(MM1 = MM1, MM2 = MM2)
#     terms$MMn <- list(MM1 = MM1n, MM2 = MM2n)
#     #terms$terms.formula <- res
#     
#     # attr(terms, "term") <- smooth.comp
#     
#     # Initialize variance components
#     init.var <- rep(1, length(g))
#     
#     res <- list(X = X, Z = Z, dim = dim, g = g, init.var = init.var)	
#     M <- cbind(res$X,res$Z)
#     
#     return(M)
#   })
#   
#   nrowss <- (unlist(lapply(multires,nrow)))
#   nranges <- (unlist(lapply(multires,ncol)))
#   
#   names(multires) <- gsub(" ",".",names(multires))
#   names(multires) <- gsub("#",".",names(multires))
#   names(multires) <- gsub("-",".",names(multires))
#   names(multires) <- gsub("/",".",names(multires))
#   names(multires) <- gsub("%",".",names(multires))
#   names(multires) <- gsub("\\(",".",names(multires))
#   names(multires) <- gsub(")",".",names(multires))
#   
#   
#   st <- 1 # for number of rows
#   st2 <- 1 # for number of column
#   end2 <- numeric() # for end of number of column
#   dataflist <- list()
#   #glist <- list()
#   for(u in 1:length(multires)){
#     prov <- multires[[u]]
#     mu <- as.data.frame(matrix(0,nrow = sum(nrowss), ncol = ncol(prov)))
#     colnames(mu) <- paste(names(multires)[u],colnames(prov), sep="_")
#     
#     nam <- paste("at",names(multires)[u],"2Dspl", sep="_")
#     end <- as.numeric(unlist(st+(nrowss[u]-1)))
#     
#     mu[st:end,] <- prov
#     dataflist[[nam]] <- as(as.matrix(mu), Class="sparseMatrix")
#     st <- end+1
#     ## for keeping track of the inits
#     # end2 <- as.numeric(unlist(st2+(ncol(prov)-1)))
#     # glist[[nam]] <- st2:end2
#     # st2 <- end2+1
#   }
#   
#   ## now build the last dataframe and adjust the glist
#   # newdatspl <- as.data.frame(do.call(cbind,dataflist))
#   # nn <- ncol(dat) # to add to the glist
#   # glist <- lapply(glist, function(x){x+nn})
#   # newdat <- data.frame(dat,newdatspl)
#   
#   ## now make the formula
#   
#   #funny <- paste(paste("grp(",names(dataflist),")",sep=""), collapse=" + ")
#   
#   ## important
#   # newdat: is the neew data frame with original data and splines per location matrices
#   # glist: is the argument to provide in group in asreml to indicate where each grouping starts and ends
#   # funny: formula to add to your random formula
#   dataflist <- lapply(dataflist,as.matrix)
#   fin <- Reduce("+",dataflist)
#   attr(fin,"variables") <- c(init0, init1)
#   # fin <-dataflist#list(newdat=dataflist, funny=funny) # 
#   return(fin)
# }