update of Nov 3rd, error in Spencers code

Author: JamesRamsay5

R/fRegress.R

@@ -73,15 +73,8 @@ fRegress.fd <- function(y, xfdlist, betalist, wt=NULL,
 
   if (is.fdPar(y)) y <- y$fd
 
-  print(class(y))
-  print("inside fRegress.fd")
-  
-  print("calling fRegressArgCheck")
-  
   arglist <- fRegressArgCheck(y, xfdlist, betalist, wt)
 
-  print("fRegressArgCheck called")
-  
   yfdobj   <- arglist$yfd
   xfdlist  <- arglist$xfdlist
   betalist <- arglist$betalist
@@ -100,7 +93,6 @@ fRegress.fd <- function(y, xfdlist, betalist, wt=NULL,
   #  ----------------------------------------------------------------
 
   #  extract dependent variable information
-  print("extract dependent variable information")
   ycoef     <- yfdobj$coefs
   ycoefdim  <- dim(ycoef)
   N         <- ycoefdim[2]
@@ -218,8 +210,6 @@ fRegress.fd <- function(y, xfdlist, betalist, wt=NULL,
 
   eigchk(Cmat)
 
-  print("Cmat, Dmat computed")
-  
   #  solve for coefficients defining BETA
 
   Lmat    <- chol(Cmat)
@@ -228,7 +218,6 @@ fRegress.fd <- function(y, xfdlist, betalist, wt=NULL,
 
   betacoef <- Cmatinv %*% Dmat
 
-  print("betacoef computed")
   #  set up fdPar objects for reg. fns. in BETAESTLIST
 
   betaestlist <- betalist
@@ -248,8 +237,6 @@ fRegress.fd <- function(y, xfdlist, betalist, wt=NULL,
     betaestlist[[j]] <- betafdParj
   }
 
-  print("betalist set up")
-
   #  set up fd objects for predicted values in YHATFDOBJ
 
   nfine   <- max(501,10*ynbasis+1)
@@ -267,8 +254,6 @@ fRegress.fd <- function(y, xfdlist, betalist, wt=NULL,
 
   df <- NA
 
-  print("yhatfdobj set up")
-  
   #  -----------------------------------------------------------------------
   #        Compute pointwise standard errors of regression coefficients
   #               if both y2cMap and SigmaE are supplied.
@@ -351,8 +336,6 @@ fRegress.fd <- function(y, xfdlist, betalist, wt=NULL,
   #                       Set up output list object
   #  -------------------------------------------------------------------
 
-  print("Set up output list object")
-  
   fRegressList <-
     list(yfdobj         = yfdobj,
          xfdlist        = xfdlist,

R/fRegress.double.R

@@ -51,11 +51,11 @@ fRegress.double <- function(y, xfdlist, betalist, wt=NULL,
   #  as predict(List).  In this call List can be any object of the
   #  "".
 
-  # Last modified 25 August 2020 by Jim Ramsay
+  # Last modified 3 November 2020 by Jim Ramsay
 
   #  check Y and compute sample size N
 
-  print("inside fRegress.numeric")
+  print("inside fRegress.double")
 
   if (!inherits(y, "numeric")) stop("Y is not a numeric vector.")
 
@@ -74,14 +74,13 @@ fRegress.double <- function(y, xfdlist, betalist, wt=NULL,
   N    <- dim(ymat)[1]
   p    <- length(xfdlist)
 
-  print("preliminaries")
+  print("computing Rmat")
 
   Zmat  <- NULL
   Rmat  <- NULL
   pjvec <- rep(0,p)
   ncoef <- 0
   for (j in 1:p) {
-    print(j)
     xfdj       <- xfdlist[[j]]
     print(class(xfdj))
     xcoef      <- xfdj$coefs
@@ -92,7 +91,6 @@ fRegress.double <- function(y, xfdlist, betalist, wt=NULL,
     pjvec[j]   <- bnbasis
     Jpsithetaj <- inprod(xbasis,bbasis)
     Zmat       <- cbind(Zmat,crossprod(xcoef,Jpsithetaj))
-    print("betafdParj$estimate")
     if (betafdParj$estimate) {
       lambdaj    <- betafdParj$lambda
       if (lambdaj > 0) {
@@ -116,10 +114,11 @@ fRegress.double <- function(y, xfdlist, betalist, wt=NULL,
   #          set up the linear equations for the solution
   #  -----------------------------------------------------------
 
-  print("Cmat and Dmat")
-  
   #  solve for coefficients defining BETA
 
+  print("assembling Cmat and Dmat")
+  
+  
   if (any(wt != 1)) {
     rtwt   <- sqrt(wt)
     Zmatwt <- Zmat*rtwt
@@ -131,23 +130,22 @@ fRegress.double <- function(y, xfdlist, betalist, wt=NULL,
     Dmat <- t(Zmat) %*% ymat
   }
 
-  #   print("solving")
-  
   eigchk(Cmat)
 
+  print("solving equation")
+  
   Cmatinv  <- solve(Cmat)
 
   betacoef <- Cmatinv %*% Dmat
 
-  
   #  compute and print degrees of freedom measure
 
   df <- sum(diag(Zmat %*% Cmatinv %*% t(Zmat)))
 
-  print("betaestlist")
-  
   #  set up fdPar object for BETAESTFDPAR
 
+  print("setting up beetaestlist")
+  
   betaestlist <- betalist
   mj2 <- 0
   for (j in 1:p) {
@@ -167,7 +165,7 @@ fRegress.double <- function(y, xfdlist, betalist, wt=NULL,
 
   #  set up fd object for predicted values
 
-  print("yhatmat")
+  print("computing yhatmat")
 
   yhatmat <- matrix(0,N,1)
   for (j in 1:p) {
@@ -200,12 +198,13 @@ fRegress.double <- function(y, xfdlist, betalist, wt=NULL,
   #               if both y2cMap and SigmaE are supplied.
   #  -----------------------------------------------------------------------
 
-  print("standard errors")
 
   if (!(is.null(y2cMap) || is.null(SigmaE))) {
 
     #  check dimensions of y2cMap and SigmaE
 
+    print("computing bvar, betastderrlist, and c2bMap")
+    
     y2cdim <- dim(y2cMap)
     if (y2cdim[2] != dim(SigmaE)[1])  stop(
       "Dimensions of Y2CMAP not correct.")

man/fRegress.Rd

@@ -527,14 +527,14 @@ dimnames(rgnContr3) <- list('', CanadianWeather$place, c('const',
 
 const365 <- create.constant.basis(c(0, 365))
 region.fd.Atlantic <- fd(matrix(rgnContr3[,,2], 1), const365)
-str(region.fd.Atlantic)
+# str(region.fd.Atlantic)
 region.fd.Continental <- fd(matrix(rgnContr3[,,3], 1), const365)
 region.fd.Pacific <- fd(matrix(rgnContr3[,,4], 1), const365)
 region.fdlist <- list(const=rep(1, 35),
      region.Atlantic=region.fd.Atlantic,
      region.Continental=region.fd.Continental,
      region.Pacific=region.fd.Pacific)
-str(TempRgn.mdl$betalist)
+# str(TempRgn.mdl$betalist)
 
 beta1 <- with(Temp.fd, fd(basisobj=basis, fdnames=fdnames))
 beta0 <- fdPar(beta1)
@@ -589,6 +589,13 @@ fRegressout <- fRegress(kneefd, xfdlist, betalist)
 all.equal(fRegressout, knee.hip.f)
 \dontshow{)}
 
+ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
+trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
+group <- gl(2,10,20, labels=c("Ctl","Trt"))
+weight <- c(ctl, trt)
+
+fRegress.D9 <- fRegress(weight ~ group)
+
 #See also the following demos:
 
 #demo('canadian-weather', package='fda')