Dear all,
I am new to SAS/IML studio and I am struggling to append data generated from R to SAS. Basically the program below which compute the power of Internal Pilot Study (pls see also attached file ) has been written in R but I need my output to be saved as a SAS dataset somewhere in my personal files directory ( "C:/Users/Jeka Abi/Documents/Results") so that I can do further analysis in SAS. I was wondering if you would be able to help me with SAS code at the end of my program to perform this. I have used these two SAS codes (run ImportMatrixFromR(Rejec, "results"); print Rejec;) ) but they just print the results from R but they do not save them as a DATASET as I want.
submit / R;
library(lattice)
library(ldbounds)
library(mvtnorm)
sen <- function(n,stages,rho,s1,s2,diff){
### variance covariance for look1 ##
sigma <- kronecker(matrix(1,nrow=1,ncol=1),matrix(c(1,rho,rho,1),nrow=2,ncol=2))
### variance covariance for look2 ###
sigma2 <- kronecker(matrix(1,nrow=2,ncol=2),matrix(c(1,rho,rho,1),nrow=2,ncol=2))
#### Noncentrality parameter for look 1 ####
mean11 <- sqrt((n*diff^2)/(s1*2))
mean21 <- sqrt((n*diff^2)/(s2*2))
#### Noncentrality parameter for look 2 ####
mean12 <- sqrt(2*(n*diff^2)/(s1*2))
mean22 <- sqrt(2*(n*diff^2)/(s2*2))
# end loop
###DEFINING FUCNTION####
out <- pmvnorm(lower=c(-Inf,-Inf),upper=c(Inf,Inf),c(mean11,mean21), corr=sigma)-pmvnorm(lower=c(-Inf,-Inf,-Inf,-Inf),upper=c(Inf,Inf,2.25,2.25),c(mean11,mean21,mean12,mean22), corr=sigma2)-0.8
out
}
rho <- 0.5
s1 <- 1
s2 <- 1
diff <- 0.5
stages <- 2
mu1 <- matrix(c(diff,diff),nrow=1,ncol=2)
mu0 <- matrix(c(0,0),nrow=1,ncol=2)
# compute boundaries
obf.bd1 <- bounds(c(1/2,2/2),iuse=1,alpha=0.0125)
boundary <- obf.bd1$upper.bounds
n1 <- floor(uniroot(sen,lower=1,upper=100,stages=stages,rho=rho,diff=diff,s1=s1,s2=s2)$root)
# simulation set-up
nsim <- 10
set.seed(1)
results <- matrix(NA,nrow=nsim,ncol=11)
colnames(results) <- c("Gesrho", "GesS1", "GesS2", "n1", "Trurho", "Estrho", "EstS1", "EstS2", "TotalN", "efficacy", "futility")
N <- vector(length=nsim)
# run simulations
for (i in 1:nsim){
results[i,1] <- rep(rho)
results[i,2] <- rep(s1)
results[i,3] <- rep(s2)
results[i,4] <- rep(2*n1)
# initiate
istop <- 0
# simulate sample of size n1
Trurho <- 0.3
results[i,5] <- rep(Trurho)
nsigma <- matrix(c(1,Trurho,Trurho,1),ncol=2,nrow=2)
Tsample <- rmvnorm(n1,c(0,0),nsigma)
Csample <- rmvnorm(n1,mu0,nsigma)
nrho <- cor(c(Tsample[,1],Csample[,1]),c(Tsample[,2],Csample[,2]))
results[i,6] <- nrho
s11 <- sqrt(var(c(Tsample[,1],Csample[,1])))
results[i,7] <- s11
s22 <- sqrt(var(c(Tsample[,2],Csample[,2])))
results[i,8] <- s22
n2 <- round(uniroot(sen,lower=1,upper=100,stages=stages,rho=nrho,diff=diff,s1=s11,s2=s22)$root)
N <- stages*n2
results[i,9] <- N
if (istop==0) {
if (n1>=N) {istop <- 1}
else {
# simulate sample of size N2-n1
Tsamplev <- rmvnorm(N-n1,c(0,0),nsigma)
Csamplev <- rmvnorm(N-n1,mu0,nsigma)
}
Tsample <- rbind(Tsample,Tsamplev)
Csample <- rbind(Csample,Csamplev)
# test statistics
z1 <- (mean(Tsample[,1])-mean(Csample[,1]))/sqrt(2/N)
z2 <- (mean(Tsample[,2])-mean(Csample[,2]))/sqrt(2/N)
efficacy <- as.numeric(z1>=(boundary[2]) | z2>=(boundary[2]))
results[i,10] <- efficacy
futility <- as.numeric(z1<(-boundary[2]) & z2<(-boundary[2]))
results[i,11] <- futility
if(results[i,10]==1 | results[i,11]==1){istop <- 1}
}# end if
} # end simulation
endsubmit;
run ImportMatrixFromR(Rejec, "results");
print Rejec;
