[R] Calling R object from R function

R. Michael Weylandt michael.weylandt at gmail.com
Thu Nov 8 17:46:14 CET 2012


On Thu, Nov 8, 2012 at 2:44 PM, Fares Said <frespider at hotmail.com> wrote:
> Hi Michael,
>
> you mentioned I can get help, can you give more details please.

Hi Fares,

I also asked that you keep our correspondance on the mailing list.
Please do so. If you want private help, you can seek a local
statistical consultancy, but free help is available here, if you
follow protocol.

>
> I also edited the example to make productive.

But I'm not using Nabble and neither are most of your other potential
respondents: if you want to use an email help list, it's much better
to use email. For instance, if we check the official R help archives
-- https://stat.ethz.ch/pipermail/r-help/2012-November/328391.html --
you'll see your post wasn't changed. It's rather wretched that Nabble
allows that, in my opinion, as it invalidates conversations.

>
> and you might tell me to create the M1 function on the get.m function but I
> prefer not because I have 9 to ten function and they have lots statement and
> I prefer to keep them separate
>

I've copied your Nabble post here:

-------------------------------------------------------

a0=rep(1,40)
a=rep(0:1,20)
b=c(rep(1,20),rep(0,20))
c0=c(rep(0,12),rep(1,28))
c1=c(rep(1,5),rep(0,35))
c2=c(rep(1,8),rep(0,32))
c3=c(rep(1,23),rep(0,17))
c4=c(rep(1,6),rep(0,34))
x=matrix(cbind(a0,a,b,c0,c1,c2,c3,c4,rnorm(40)),nrow=40,ncol=9)
Dat <- cbind(x,rnorm(40,2,30));colnames(Dat) <-
c("a0","a","b","c0","c1","c2","c3","c4","c5","Y")

M1 <- function(Trdat,Tedat,mdat,nsam,conv){
M1list <- NULL
vectx <- c(1,4,6,7,9)
vectz <- c(2,3,7,1,4,9)
X <- Trdat[,vectx]
Z <- Trdat[,vectz]
Y <- Trdat[,ncol(Trdat)]
TesX <- Tedat[,vectx]
TesZ <- Tedat[,vectz]
TesY <- Tedat[,ncol(Tedat)]
Treig <- eigen(crossprod(X))$values
if(any(abs(Treig) < conv))stop("In M1 the design matrix (X) is
singular for simulation ",paste(nsam))
Comp <- c("nCol(X)"= ncol(X),"nCol(Z)"= ncol(Z),"Is
length(Y)=nrow(X)"= length(Y)==nrow(X),
          "Is length(Y)=nrow(Z)"= length(Y)==nrow(Z))

M1list$vectx <- vectx
M1list$vectz <- vectz
M1list$X <- X
M1list$Z <- Z
M1list$Y <- Y
M1list$TesX <- TesX
M1list$TesZ <- TesZ
M1list$TesY <- TesY
M1list$Comp <- Comp
return(M1list)
}


get.m <- function(dat,asim,ModelFun,M,conv){
Sim <- list()
modInd <- ModelFun(Trdat=dat,Tedat=dat,mdat=dat,nsam=-1,conv=conv)  #
HERE WHERE I NEED HELP i only need to import vectx and vectz  that is
why I set Trdat=Tedat=dat
if(M==1){
vecx <- modInd$vectx
vecz <- modInd$vectz
px <- length(vecx)
pz <- length(vecz)
pk <- length(modInd$Comp)
nam <-colnames(dat[,vecx])
Asse <- matrix(NA,nrow=asim,ncol=px)
Check <- matrix(NA,nrow=pk,ncol=asim)
colnames(Check) <- paste("CheckIter",1:asim,sep="")
}
else {
vecx <- modInd$vectx
vecz <- modInd$vectz
px <- length(vecx)
pz <- length(vecz)
pk <- length(modInd$Comp)
nam <-colnames(dat[,vecx])
Asse <- matrix(NA,nrow=asim,ncol=px)
Check <- matrix(NA,nrow=pk,ncol=asim)
colnames(Check) <- paste("CheckIter",1:asim,sep="")
}

for(k in 1:asim){
cat("Iter #",paste(k),"\n")
#==========================================================================================
#                         Start Sampling code
#==========================================================================================
# Sample the Index for Train Set
set.seed(k)
Indx<-sample(1:nrow(dat),nrow(dat),replace=T)
SamDat <- dat[Indx,]
# Split Data
set.seed(k)
TrainInd <- sample(1:nrow(SamDat), trunc(2*length(1:nrow(SamDat))/3))
# Sample 2/3 of the data
TrSet <- SamDat[TrainInd,]  # Train data
 ######## Hold 1/3 of the data
TeSet <- SamDat[-TrainInd,]      # hold 1/3 of the data
Trind <- ceiling((2*length(Indx))/3)
Model <- ModelFun(Trdat=TrSet,Tedat=TeSet,mdat=dat,nsam=k,conv=conv)
Y <- Model$Y
X <- Model$X
Z <- Model$Z
TesX <- Model$TesX
TesZ <- Model$TesZ
TesY <- Model$TesY
xnam <-colnames(X)
znam <-colnames(Z)
pc <- ncol(X)
fmla <- as.formula(paste("Y ~",paste(xnam, collapse= "+"),"-1",sep=""))
fitlm <- lm(formula=fmla,data = data.frame(cbind(X,Y)))
ResiSqr <- (residuals(fitlm))*(residuals(fitlm))
Check[,k] <- Model$Comp
Asse[k,1:pc] <- coef(fitlm)

   }
Sim$Check <- Check
Sim$Asse <- Asse
return(Asse)
}
get.m(dat=Dat,asim=6,ModelFun=M1,M=1,conv=1e-4)

------------------------------------------------------

Your code is quite messy and difficult to follow: I'd write it as follows:

# Initialize some data
Dat <- cbind(
    a0 = rep(1, 40),
    a = rep(0:1, 20),
    b = rep(c(1,0), each = 20),
    c0=c(rep(0,12),rep(1,28)),
    c1=c(rep(1,5),rep(0,35)),
    c2=c(rep(1,8),rep(0,32)),
    c3=c(rep(1,23),rep(0,17)),
    c4=c(rep(1,6),rep(0,34)),
    Y = rnorm(40,2,30))

M1 <- function(Trdat,Tedat,mdat,nsam,conv){
    vectx <- c(1,4,6,7,9)
    vectz <- c(2,3,7,1,4,9)

    X <- Trdat[,vectx]
    Z <- Trdat[,vectz]
    Y <- Trdat[,ncol(Trdat)]

    TesX <- Tedat[,vectx]
    TesZ <- Tedat[,vectz]
    TesY <- Tedat[,ncol(Tedat)]

    Treig <- eigen(crossprod(X))$values

    if(any(abs(Treig) < conv))
        stop("In M1 the design matrix (X) is singular for simulation ", nsam)

    Comp <- c(`nCol(X)`= ncol(X),`nCol(Z)`= ncol(Z),
                        `Is length(Y)=nrow(X)` = length(Y)==nrow(X),
                        `Is length(Y)=nrow(Z)` = length(Y)==nrow(Z))

    list(vectx = vectx,
          vectz = vectz,
          X = X, Z = Z, Y = Y,
          TesX = TesX, TesZ = TesZ,
          TesY = TesY, Comp = Comp)
}

with similar modifications throughout the rest of your code. Hopefully
you find that easier to read: I certainly do.

Now, it's still not entirely clear to me what you're asking: but
perhaps global assignment and/or lexical scoping does what you want?

Consider these tricks:

#########################
add_magic_Z <- function(z){
    Z <- z
    function(x) x + z
}

add_2 <- function(2)

add_2(5) # 7

add_3 <- function(3)

add_3(0) # 3
##########################
f <- function(x) {
    X <<- X + x
    return(x^2)
}

X <- -3

f(5)
f(10)

print(X) # increased by the `<<-` inside of f()
###########################

The latter is considered bad programming practice however.

Cheers,
Michael




More information about the R-help mailing list