[R] Conflict between Accelerate BLAS and package 'stats'
Hi All,
After successfully setting up the Accelerate BLAS library via
/cd /Library/Frameworks/R.framework/Resources/lib
ln -sf
/System/Library/Frameworks/Accelerate.framework/Frameworks/vecLib.framework/Versions/Current/libBLAS.dylib
libRblas.dylib
/
I restarted R (version 3.0.1) and got the following error message:
/Error in dyn.load(file, DLLpath = DLLpath, ...) :
unable to load shared object
'/Library/Frameworks/R.framework/Versions/3.0/Resources/library/stats/libs/stats.so':
dlopen(/Library/Frameworks/R.framework/Versions/3.0/Resources/library/stats/libs/stats.so,
6): Symbol not found: _lsame_
Referenced from:
/Library/Frameworks/R.framework/Versions/3.0/Resources/lib/libRlapack.dylib
Expected in:
/System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
in
/Library/Frameworks/R.framework/Versions/3.0/Resources/lib/libRlapack.dylib
During startup - Warning message:
package 'stats' in options("defaultPackages") was not found /
Can anybody help?
Below is the system info for my iMac:
/sysname "Darwin"
release "10.8.0"
version
"Darwin Kernel Version 10.8.0: Tue Jun 7 16:33:36 PDT 2011;
root:xnu-1504.15.3~1/RELEASE_I386"
machine"i386"
/
Thanks,
David
David Degras
Department of Mathematical Sciences
DePaul University
--
View this message in context:
http://r.789695.n4.nabble.com/Conflict-between-Accelerate-BLAS-and-package-stats-tp4678674.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Data Manipulation in R
Hi,
May be this helps:
Y1 <- read.table(text="V1 V2 V3 V4
1 4 0 20 17
2 4 0 15 17
3 2 0 13 21",sep="",header=TRUE)
Y2 <- read.table(text="V1 V2 V3 V4
1 20 52 15 18
2 18 54 14 21
3 18 51 13 21",sep="",header=TRUE)
res <- lapply(seq_len(nrow(Y1)),function(i) {dat <-
data.frame(X=i,Y1=unlist(Y1[i,]),Y2=unlist(Y2[i,])); row.names(dat) <-
1:nrow(dat);
write.csv(dat,paste0("file",i,".csv"),row.names=FALSE,quote=FALSE)})
A.K.
A.K.
On Monday, October 21, 2013 12:24 AM, Anamika Chaudhuri
wrote:
Hi:
I am looking for some help to manipulate data in R. I have two csv files.
datasetY1
V1 "V2" "V3" "V4"
1 4 0 20 17
2 4 0 15 17
3 2 0 13 21
datasetY2
V1 "V2" "V3" "V4"
1 20 52 15 18
2 18 54 14 21
3 18 51 13 21
I want to be able to create separate csv files by taking the corresponding
rows of dataset1 and dataset2, convert them into columns. So from the above
example I would be creating 3 datasets (csvs), of which the first one would
be
X Y1 Y2 1 4 20 1 0
52 1 20 15 1 17
18
Appreciate any help.
Thanks
Anamika
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] nlminb() - how do I constrain the parameter vector properly?
Bill: I didn't look at the code but I think the OP means that during the
nlminb algorithm,
the variance covariance parameters hit values such that the covariance
matrix estimate becomes negative definite.
Again, I didn't look at the details but one way to deal with this is to
have the likelihood
function return -Inf whenever the covariance matrix becomes not positive
definite. so, the
likelihood should check for positive definiteness first before it actually
calculates anything.
If PD is not true, the -Inf value should push nlminb towards values that
obtain a positive definite matrix. But there could be something more subtle
going on that I'm not understanding. I don't know even what algorithm
nlminb is using ( probably quasi-newton ) but this is one thing the OP
could try.
On Sun, Oct 20, 2013 at 9:41 PM, William Dunlap wrote:
> Do you mean that your objective function (given to nlminb) parameterized
> a positive definite matrix, P, as the elements in its upper (or lower)
> triangle?
> If so, you could reparameterize it by the non-zero (upper triangular)
> elements
> of the Choleski decomposition, C, of P. Compute P as crossprod(C), compute
> the initial estimate of C as chol(P).
>
> Bill Dunlap
> Spotfire, TIBCO Software
> wdunlap tibco.com
>
>
> > -Original Message-
> > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
> On Behalf
> > Of Steven LeBlanc
> > Sent: Sunday, October 20, 2013 3:02 PM
> > To: r-help@R-project.org
> > Subject: [R] nlminb() - how do I constrain the parameter vector properly?
> >
> > Greets,
> >
> > I'm trying to use nlminb() to estimate the parameters of a bivariate
> normal sample and
> > during one of the iterations it passes a parameter vector to the
> likelihood function
> > resulting in an invalid covariance matrix that causes dmvnorm() to throw
> an error. Thus,
> > it seems I need to somehow communicate to nlminb() that the final three
> parameters in
> > my parameter vector are used to construct a positive semi-definite
> matrix, but I can't see
> > how to achieve this using the constraint mechanism provided. Additional
> details are
> > provided in the code below.
> >
> > Suggestions?
> >
> > Best Regards,
> > Steven
> >
> > Generate the data set I'm using:
> >
> > mu<-c(1,5)
> > sigma<-c(1,2,2,6)
> > dim(sigma)<-c(2,2)
> > set.seed(83165026)
> > sample.full<-sample.deleted<-mvrnorm(50,mu,sigma)
> > delete.one<-c(1,5,13,20)
> > delete.two<-c(3,11,17,31,40,41)
> > sample.deleted[delete.one,1]<-NA
> > sample.deleted[delete.two,2]<-NA
> > missing<-c(delete.one,delete.two)
> > complete<-sample.deleted[!(is.na(sample.deleted[,1]) | is.na
> (sample.deleted[,2])),]
> > deleted<-sample.deleted[missing,]
> >
> > Try to estimate the parameters of the data set less the deleted values:
> >
> > exact<-function(theta,complete,deleted){
> > one.only<-deleted[!(is.na(deleted[,1])),1]
> > two.only<-deleted[!(is.na(deleted[,2])),2]
> > u<-c(theta[1],theta[2])
> > sigma<-c(theta[3],theta[5],theta[5],theta[4])
> > dim(sigma)<-c(2,2)
> > -sum(log(dmvnorm(x=complete,mean=u,sigma=sigma)))-
> > sum(log(dnorm(one.only,u[1],sigma[1,1])))-
> > sum(log(dnorm(two.only,u[2],sigma[2,2])))
> > }
> >
> nlminb(start=c(0,0,1,1,0),objective=exact,complete=complete,deleted=deleted,control=l
> > ist(trace=1))
> >
> > Escape and it stops at Iteration 9 on my machine.
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Data Manipulation in R
Hi: I am looking for some help to manipulate data in R. I have two csv files. datasetY1 V1 "V2" "V3" "V4" 1 4 0 20 17 2 4 0 15 17 3 2 0 13 21 datasetY2 V1 "V2" "V3" "V4" 1 20 52 15 18 2 18 54 14 21 3 18 51 13 21 I want to be able to create separate csv files by taking the corresponding rows of dataset1 and dataset2, convert them into columns. So from the above example I would be creating 3 datasets (csvs), of which the first one would be X Y1Y2 1 4 20 1 0 52 1 20 15 1 17 18 Appreciate any help. Thanks Anamika [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Loop for taking sum of rows based on proximity to other non-NA rows
Sorry, I noticed that when two "Count" values are the same and NA in between,
my function fails.
#Modified
fun1 <- function(dat,n) {
rl <- rle(is.na(dat[,"Count"]))
indx <-
which(is.na(dat[,"Count"]))[rep(rl$lengths[rl$values],rl$lengths[rl$values])==n]
lst1 <- lapply(split(indx,((seq_along(indx)-1)%/%n)+1),function(x) {
x1 <- dat[c(min(x)-1L,x,max(x)+1L),]
x2 <- x1[!is.na(x1$Count),]
datN <- data.frame(Position=x2$Position[x2$Count %in%
max(x2$Count)],Count=sum(x2$Count))
rowN <- row.names(x2)[x2$Count %in% max(x2$Count)]
datN<- if(length(rowN)>1) datN[1,] else datN
row.names(datN) <- if(length(rowN) >1) rowN[1] else rowN
datN
})
names(lst1) <- NULL
lst1 <- lst1[!duplicated(sapply(lst1,row.names))]
dat2 <- do.call(rbind,lst1)
indx2 <-
sort(unlist(lapply(split(indx,((seq_along(indx)-1)%/%n)+1),function(x)
c(min(x)-1L,x,c(max(x)+1L))),use.names=FALSE))
dat1New <- dat[-indx2[!indx2 %in% row.names(dat2)],]
dat1New[match(row.names(dat2),row.names(dat1New)),] <- dat2
row.names(dat1New) <- 1:nrow(dat1New)
dat1New
}
#
fun2 <- function(dat,n){
indx <- cumsum(c(1,abs(diff(is.na(dat[,"Count"])
indx1 <- indx[is.na(dat[,"Count"])]
names(indx1) <- which(is.na(dat[,"Count"]))
indx2 <- indx1[indx1 %in% names(table(indx1))[table(indx1)==n]]
lst1 <- tapply(seq_along(indx2),list(indx2),FUN=function(i) {
x1 <- indx2[i]
x2 <- as.numeric(names(x1))
x3 <- dat[c(min(x2)-1L,x2,max(x2)+1L),]
x4 <- subset(x3, !is.na(Count))
x5 <- data.frame(Position=x4$Position[x4$Count
%in% max(x4$Count)],Count=sum(x4$Count))
ind <- x4$Count %in% max(x4$Count)
row.names(x5) <- row.names(x4)[ind]
x5 <- if(sum(ind)>1) x5[1,] else x5
x5
})
attr(lst1,"dimnames") <- NULL
dat2 <- do.call(rbind,lst1)
indx3 <- sort(unlist(tapply(seq_along(indx2),list(indx2),FUN=function(i) {x1 <-
indx2[i]
x2 <- as.numeric(names(x1))
c(min(x2)-1L, x2,
max(x2)+1L)}),use.names=FALSE))
dat$id <- 1:nrow(dat)
dat2$id <- as.numeric(row.names(dat2))
library(plyr)
res <- join(dat,dat2[,-1],by="id",type="left")
res1 <- res[!((row.names(res) %in% indx3) & is.na(res[,4])),]
res1[,2][!is.na(res1[,4])] <- res1[,4][!is.na(res1[,4])]
res2 <- res1[,1:2]
row.names(res2) <- 1:nrow(res2)
res2
}
identical(fun1(dat1,1),fun2(dat1,1))
#[1] TRUE
identical(fun1(fun1(dat1,1),2),fun2(fun2(dat1,1),2))
#[1] TRUE
identical(fun1(fun1(fun1(dat1,1),2),3),fun2(fun2(fun2(dat1,1),2),3))
#[1] TRUE
fun1(fun1(fun1(dat1,1),2),3)
# Position Count
#1 61 37
#2 18 62
#3 42 65
##When I tried the function on a bigger dataset:
set.seed(185)
datT <- data.frame(Position = sample(10:80,1e5,replace=TRUE),Count=
sample(c(NA, 10:100),1e5, replace=TRUE))
dim(datT)
#[1] 10 2
system.time(res <- fun1(datT,1))
# user system elapsed
# 0.708 0.000 0.709
system.time(res2 <- fun2(datT,1))
# user system elapsed
# 1.400 0.016 1.421
system.time(res3 <- removeNNAs(datT,1))
# user system elapsed
# 1.068 0.000 1.071
all.equal(res,res2)
#[1] TRUE
all.equal(res,res3)
#[1] "Attributes: < Component 2: Numeric: lengths (97786, 97778) differ >"
#[2] "Component 1: Numeric: lengths (97786, 97778) differ"
#[3] "Component 2: Numeric: lengths (97786, 97778) differ"
dim(res)
#[1] 97786 2
dim(res3)
#[1] 97778 2
##Here your function seems to give the correct number of rows as:
rl <- rle(is.na(datT[,"Count"]))
indx <-
which(is.na(datT[,"Count"]))[rep(rl$lengths[rl$values],rl$lengths[rl$values])==1]
dim(datT)[1]- 2*length(indx)
#[1] 97778
#Here is where I think the difference occur (in addition to the one with the
values)
datS <- datT[16000:2,]
row.names(datS) <- 1:nrow(datS)
resT <- fun1(datS,1)
resT3 <- removeNNAs(datS,1)
datS[3402:3408,]
Position Count
3402 72 70
3403 38 51
3404 80 NA
3405 26 44
3406 42 NA
3407 78 77
3408 70 89
resT3[3311:3318,]
Position Count
3401 54 65
3402 72 70
3407 78 172##
3408 70 89
3409 27 40
3410 44 44
3411 73 75
3412 73 76
resT[3311:3318,]
Position Count
3311 29 98
3312 54 65
3313 72 70
3314 38 95
3315 78 121 ###
3316 70 89
3317 27 40
3318 44 44
In these conditions, the post is not very clear about dealing it.
A.K.
On Sunday, October 20, 2013 9:36 PM, arun wrote:
Hi Jeff,
I found some differe
Re: [R] XML package not working
I'm running the same os system and same r version. Only difference I see in
Sys.info() is that login, user, and effective_user on my machine is all are
of the same case. Your display shows login as upper case. I'm not even
sure if that would make a difference.
With gratitude,
CEO'Riley Jr.
Charles Ellis O'Riley Jr.
Ambition is a state of permanent dissatisfaction with the present
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Steven Dwayne Randolph
Sent: Sunday, October 20, 2013 8:24 AM
To: Duncan Murdoch; r-help@r-project.org; Ista Zahn; Duncan Murdoch
Cc: Steven Dwayne Randolph; stevendrando...@aol.com
Subject: Re: [R] XML package not working
My apologies for not conforming to the posting guideline.
Sys.info()
sysname release
version
"Windows" "7 x64" "build 7601,
Service Pack 1"
nodename machine
login
"xxNU247BZ1S" "x86-64"
"XX"
user effective_user
"xxx""xxx"
When I attempt to install a local copy of the xml.zip file:
in read.dcf(file.path(pkgname, "DESCRIPTION"), c("Package", "Type")) :
cannot open the connection
In addition: Warning messages:
1: In unzip(zipname, exdir = dest) : error 1 in extracting from zip file
2: In read.dcf(file.path(pkgname, "DESCRIPTION"), c("Package", "Type")) :
cannot open compressed file 'XML/DESCRIPTION', probable reason 'No such
file or directory'
-Original Message-
From: Duncan Murdoch [mailto:murdoch.dun...@gmail.com]
Sent: Saturday, October 19, 2013 8:23 PM
To: Steven Dwayne Randolph; r-help@r-project.org
Cc: stevendrando...@aol.com
Subject: Re: [R] XML package not working
On 13-10-19 3:47 PM, Steven Dwayne Randolph wrote:
> I know I cannot be the only one who is not able to install the XML
package from CRAN (zip or tar file) Many packages depend on this XML
package. Can someone help me either access the source for a good binary? I
have received no assistance from the author/developer of the package.
It installs fine for me.
Duncan Murdoch
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] nlminb() - how do I constrain the parameter vector properly?
Do you mean that your objective function (given to nlminb) parameterized
a positive definite matrix, P, as the elements in its upper (or lower) triangle?
If so, you could reparameterize it by the non-zero (upper triangular) elements
of the Choleski decomposition, C, of P. Compute P as crossprod(C), compute
the initial estimate of C as chol(P).
Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com
> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
> Behalf
> Of Steven LeBlanc
> Sent: Sunday, October 20, 2013 3:02 PM
> To: r-help@R-project.org
> Subject: [R] nlminb() - how do I constrain the parameter vector properly?
>
> Greets,
>
> I'm trying to use nlminb() to estimate the parameters of a bivariate normal
> sample and
> during one of the iterations it passes a parameter vector to the likelihood
> function
> resulting in an invalid covariance matrix that causes dmvnorm() to throw an
> error. Thus,
> it seems I need to somehow communicate to nlminb() that the final three
> parameters in
> my parameter vector are used to construct a positive semi-definite matrix,
> but I can't see
> how to achieve this using the constraint mechanism provided. Additional
> details are
> provided in the code below.
>
> Suggestions?
>
> Best Regards,
> Steven
>
> Generate the data set I'm using:
>
> mu<-c(1,5)
> sigma<-c(1,2,2,6)
> dim(sigma)<-c(2,2)
> set.seed(83165026)
> sample.full<-sample.deleted<-mvrnorm(50,mu,sigma)
> delete.one<-c(1,5,13,20)
> delete.two<-c(3,11,17,31,40,41)
> sample.deleted[delete.one,1]<-NA
> sample.deleted[delete.two,2]<-NA
> missing<-c(delete.one,delete.two)
> complete<-sample.deleted[!(is.na(sample.deleted[,1]) |
> is.na(sample.deleted[,2])),]
> deleted<-sample.deleted[missing,]
>
> Try to estimate the parameters of the data set less the deleted values:
>
> exact<-function(theta,complete,deleted){
> one.only<-deleted[!(is.na(deleted[,1])),1]
> two.only<-deleted[!(is.na(deleted[,2])),2]
> u<-c(theta[1],theta[2])
> sigma<-c(theta[3],theta[5],theta[5],theta[4])
> dim(sigma)<-c(2,2)
> -sum(log(dmvnorm(x=complete,mean=u,sigma=sigma)))-
> sum(log(dnorm(one.only,u[1],sigma[1,1])))-
> sum(log(dnorm(two.only,u[2],sigma[2,2])))
> }
> nlminb(start=c(0,0,1,1,0),objective=exact,complete=complete,deleted=deleted,control=l
> ist(trace=1))
>
> Escape and it stops at Iteration 9 on my machine.
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Loop for taking sum of rows based on proximity to other non-NA rows
Hi Jeff,
I found some difference in results between your function and mine. It also
point out a mistake in my code. In the original post, it says:
"""
I need to write a loop to march down the rows, and if there are 2 rows in
"Count"
where there is only 1 NA row between them, sum the two values up and
print only one row with the summed Count value and the Position value that
corresponds to the larger Count value, thus making the three
rows into one.
"
Sorry, I read it incorrectly the last time and selected the maximum "Position"
value instead of that corresponds to the larger Count value.
After correcting the function, there is still some difference between the
results.
##fun1() and fun2() corrected
fun1 <- function(dat,n) {
rl <- rle(is.na(dat[,"Count"]))
indx <-
which(is.na(dat[,"Count"]))[rep(rl$lengths[rl$values],rl$lengths[rl$values])==n]
lst1 <- lapply(split(indx,((seq_along(indx)-1)%/%n)+1),function(x) {
x1 <- dat[c(min(x)-1L,x,max(x)+1L),]
x2 <- x1[!is.na(x1$Count),]
datN <- data.frame(Position=x2$Position[x2$Count %in%
max(x2$Count)],Count=sum(x2$Count))
rowN <- row.names(x2)[x2$Count %in% max(x2$Count)]
row.names(datN) <- if(length(rowN)>1) rowN[1] else rowN
datN
})
names(lst1) <- NULL
lst1 <- lst1[!duplicated(sapply(lst1,row.names))] ##added
dat2 <- do.call(rbind,lst1)
indx2 <-
sort(unlist(lapply(split(indx,((seq_along(indx)-1)%/%n)+1),function(x)
c(min(x)-1L,x,c(max(x)+1L))),use.names=FALSE))
dat1New <- dat[-indx2[!indx2 %in% row.names(dat2)],]
dat1New[match(row.names(dat2),row.names(dat1New)),] <- dat2
row.names(dat1New) <- 1:nrow(dat1New)
dat1New
}
##
fun2 <- function(dat,n){
indx <- cumsum(c(1,abs(diff(is.na(dat[,"Count"])
indx1 <- indx[is.na(dat[,"Count"])]
names(indx1) <- which(is.na(dat[,"Count"]))
indx2 <- indx1[indx1 %in% names(table(indx1))[table(indx1)==n]]
lst1 <- tapply(seq_along(indx2),list(indx2),FUN=function(i) {
x1 <- indx2[i]
x2 <- as.numeric(names(x1))
x3 <- dat[c(min(x2)-1L,x2,max(x2)+1L),]
x4 <- subset(x3, !is.na(Count))
x5 <- data.frame(Position=x4$Position[x4$Count
%in% max(x4$Count)],Count=sum(x4$Count))
ind <- x4$Count %in% max(x4$Count)
row.names(x5) <- if(sum(ind)>1)
row.names(x4)[ind][1] else row.names(x4)[ind]
x5
})
attr(lst1,"dimnames") <- NULL
dat2 <- do.call(rbind,lst1)
indx3 <- sort(unlist(tapply(seq_along(indx2),list(indx2),FUN=function(i) {x1 <-
indx2[i]
x2 <- as.numeric(names(x1))
c(min(x2)-1L, x2,
max(x2)+1L)}),use.names=FALSE))
dat$id <- 1:nrow(dat)
dat2$id <- as.numeric(row.names(dat2))
library(plyr)
res <- join(dat,dat2[,-1],by="id",type="left")
res1 <- res[!((row.names(res) %in% indx3) & is.na(res[,4])),]
res1[,2][!is.na(res1[,4])] <- res1[,4][!is.na(res1[,4])]
res2 <- res1[,1:2]
row.names(res2) <- 1:nrow(res2)
res2
}
dat1 <- structure(list(Position = c(15L, 22L, 38L, 49L, 55L, 61L, 62L,
14L, 29L, 63L, 46L, 22L, 18L, 24L, 22L, 49L, 42L, 38L, 29L, 22L,
29L, 23L, 42L), Count = c(15L, NA, NA, 5L, NA, 17L, 18L, NA,
NA, NA, 8L, NA, 20L, NA, NA, 16L, 19L, NA, NA, NA, 13L, NA, 33L
)), .Names = c("Position", "Count"), class = "data.frame", row.names = c(NA,
-23L))
fun1(dat1,1)
Position Count
1 15 15
2 22 NA
3 38 NA
4 61 22
5 62 18
6 14 NA
7 29 NA
8 63 NA
9 18 28 ###
10 24 NA
11 22 NA
12 49 16
13 42 19
14 38 NA
15 29 NA
16 22 NA
17 42 46
removeNNAs(dat1,1) #gets similar results
#but,
fun1(fun1(dat1,1),2)
Position Count
1 61 37
2 62 18
3 14 NA
4 29 NA
5 63 NA
6 18 44 ###different
7 42 19
8 38 NA
9 29 NA
10 22 NA
11 42 46
removeNNAs(dat1,2,lessOrEqual=TRUE)
Position Count
6 61 37
7 62 18
8 14 NA
9 29 NA
10 63 NA
16 49 44 ## different
17 42 19
18 38 NA
19 29 NA
20 22 NA
23 42 46
>
removeNNAs(dat1,3,lessOrEqual=TRUE)
Position Count
6 61 37
16 49 62
23 42 65
fun1(fun1(fun1(dat1,1),2),3)
Position Count
1 61 37
2 18 62
3 42 65
A.K.
On Sunday, October 20, 2013 7:49 PM, Jeff Newmiller
wrote:
Looks like a right parenthesis was dropped. Corrected:
removeNNAs <- function( dat, N, lessOrEqual=FALSE ) {
N1 <- N+1
rx <- rle( !is.na( d
Re: [R] XML package not working
My apologies for not conforming to the posting guideline.
Sys.info()
sysname release
version
"Windows" "7 x64" "build 7601, Service
Pack 1"
nodename machine
login
"xxNU247BZ1S" "x86-64"
"XX"
user effective_user
"xxx""xxx"
When I attempt to install a local copy of the xml.zip file:
in read.dcf(file.path(pkgname, "DESCRIPTION"), c("Package", "Type")) :
cannot open the connection
In addition: Warning messages:
1: In unzip(zipname, exdir = dest) : error 1 in extracting from zip file
2: In read.dcf(file.path(pkgname, "DESCRIPTION"), c("Package", "Type")) :
cannot open compressed file 'XML/DESCRIPTION', probable reason 'No such file
or directory'
-Original Message-
From: Duncan Murdoch [mailto:murdoch.dun...@gmail.com]
Sent: Saturday, October 19, 2013 8:23 PM
To: Steven Dwayne Randolph; r-help@r-project.org
Cc: stevendrando...@aol.com
Subject: Re: [R] XML package not working
On 13-10-19 3:47 PM, Steven Dwayne Randolph wrote:
> I know I cannot be the only one who is not able to install the XML package
> from CRAN (zip or tar file) Many packages depend on this XML package. Can
> someone help me either access the source for a good binary? I have received
> no assistance from the author/developer of the package.
It installs fine for me.
Duncan Murdoch
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] "Error : requires numeric/complex matrix/vector arguments"
Can you post my message?ThanksValentina Dear R users, I hope that this format is find. I'm a new user of R. I'm trying to do a LM test as in Terasvirta (2013) and when I run the code there is this type of error: Error in t(mX) %*% mX : requires numeric/complex matrix/vector arguments. The problem could be on how I construct by dataset. I have import the data by the following command getwd() [1] "C:/Program Files/R" > list() list() > list.files() [1] "colombo.csv" "colombov.csv" "colombov.txt" "R-3.0.2" > D<-read.table(file="colombov.txt", header=T) In the file colombov there are 9 variables for T=126. > names(X) [1] "DATE" "G" "T" "Y" "SV" "GperSV" "TperSV" "YperSV" [9] "news" >fix(D) From my D I have to construct 3 matrices. To do that I go on with the further steps: First of all---> I need 3 matrix(mY, mX, mZ) I create a matrix for each of the above 9 variables. G<-matrix(D$G) I repeat this step for all the above varibles. (I repeat this step n=9) Then, I use mY<-matrix(nrow=126, ncol=4) mY<-cbind(G, T, Y, news) I repeat the above step to create matrix mX and mZ. Then I check the class(mZ) class(mY) and class(mY). It returns that mZ, mY, mX are matrix. Is there something wrong when I construct my dataset? Is for this reason that I cannot compute my LM test and it returns "Error : requires numeric/complex matrix/vector argumentsâ"? Any suggestions Thanks a lot Valentina Any suggestion? Thanks a lot Valentina __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nlminb() - how do I constrain the parameter vector properly?
On Oct 20, 2013, at 3:01 PM, Steven LeBlanc wrote:
> Greets,
>
> I'm trying to use nlminb() to estimate the parameters of a bivariate normal
> sample and during one of the iterations it passes a parameter vector to the
> likelihood function resulting in an invalid covariance matrix that causes
> dmvnorm() to throw an error. Thus, it seems I need to somehow communicate to
> nlminb() that the final three parameters in my parameter vector are used to
> construct a positive semi-definite matrix, but I can't see how to achieve
> this using the constraint mechanism provided. Additional details are provided
> in the code below.
>
I recently noticed that teh dmvnorn function in package mixtools doe not throw
an error when given a non-positive definite varaince-covariance matrix. Whether
it should be doing so might be up for debate. Peter Dalgaard seemed to think
that any self respecting such function _should_ be throwing NaN's back at you.
The code it was using is rather compact:
> mixtools::dmvnorm
function (y, mu = NULL, sigma = NULL)
{
exp(logdmvnorm(y, mu = mu, sigma = sigma))
}
> mixtools::logdmvnorm
function (y, mu = NULL, sigma = NULL)
{
if (is.vector(y))
y <- matrix(y, nrow = 1)
d <- ncol(y)
if (!is.null(mu))
y <- sweep(y, 2, mu, "-")
if (is.null(sigma))
sigma <- diag(d)
k <- d * 1.83787706640935
a <- qr(sigma)
logdet <- sum(log(abs(diag(a$qr
if (nrow(y) == 1)
mahaldist <- as.vector(y %*% qr.solve(a, t(y)))
else mahaldist <- rowSums((y %*% qr.solve(a)) * y)
-0.5 * (mahaldist + logdet + k)
}
I noticed that you were logging the dmvnorm values and so I wondered also if
going directly to that function might save some time. (Note that I am way out
of my mathematical dept here. Caveat emptor, maxima caveat emptor).
--
David.
> Suggestions?
>
> Best Regards,
> Steven
>
> Generate the data set I'm using:
>
> mu<-c(1,5)
> sigma<-c(1,2,2,6)
> dim(sigma)<-c(2,2)
> set.seed(83165026)
> sample.full<-sample.deleted<-mvrnorm(50,mu,sigma)
> delete.one<-c(1,5,13,20)
> delete.two<-c(3,11,17,31,40,41)
> sample.deleted[delete.one,1]<-NA
> sample.deleted[delete.two,2]<-NA
> missing<-c(delete.one,delete.two)
> complete<-sample.deleted[!(is.na(sample.deleted[,1]) |
> is.na(sample.deleted[,2])),]
> deleted<-sample.deleted[missing,]
>
> Try to estimate the parameters of the data set less the deleted values:
>
> exact<-function(theta,complete,deleted){
>one.only<-deleted[!(is.na(deleted[,1])),1]
>two.only<-deleted[!(is.na(deleted[,2])),2]
>u<-c(theta[1],theta[2])
>sigma<-c(theta[3],theta[5],theta[5],theta[4])
>dim(sigma)<-c(2,2)
>-sum(log(dmvnorm(x=complete,mean=u,sigma=sigma)))-
>sum(log(dnorm(one.only,u[1],sigma[1,1])))-
>sum(log(dnorm(two.only,u[2],sigma[2,2])))
> }
> nlminb(start=c(0,0,1,1,0),objective=exact,complete=complete,deleted=deleted,control=list(trace=1))
>
> Escape and it stops at Iteration 9 on my machine.
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
David Winsemius
Alameda, CA, USA
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Loop for taking sum of rows based on proximity to other non-NA rows
Looks like a right parenthesis was dropped. Corrected:
removeNNAs <- function( dat, N, lessOrEqual=FALSE ) {
N1 <- N+1
rx <- rle( !is.na( dat$Count ) )
# indexes of the ends of each run of NAs or non-NAs
cs <- cumsum( rx$lengths )
# indexes of the ends of runs of NAs or non-NAs
cs2 <- cs[ !rx$values ]
# If the first Count is NA, then drop first run of NAs
if ( !rx$values[1] ) {
cs2 <- cs2[ -1 ]
}
# If the last Count is NA, then drop last run of NAs
if ( !rx$values[ length( rx$values ) ] ) {
cs2 <- cs2[ -length( cs2 ) ]
}
# cs2 is indexes of rows to potentially receive deleted Counts
# after collapse
cs2 <- cs2 + 1
# cs1 is indexes of non-NA Counts to be deleted
cs1 <- cs[ rx$values ][ seq.int( length( cs2 ) ) ]
# identify the indexes of the Count values before the strings
# of NAs that meet the criteria
if ( lessOrEqual ) {
idx0 <- N1 >= ( cs2 - cs1 )
} else {
idx0 <- N1 == ( cs2 - cs1 )
}
idx1 <- cs1[ idx0 ]
# identify the indexes of the Count values after the strings of
# NAs that meet the criteria
idx2 <- cs2[ idx0 ]
# Identify which indexes are both sources and destinations
idx1c <-c( idx2[ -length( idx2 ) ] == idx1[ -1 ], FALSE )
# identify groups of indexes that need to be merged
idx1g <- rev( cumsum( rev( !idx1c ) ) )
# find which elements of idx1 represent the beginning of a
# sequence of indexes to be replaced (meta-indexes)
srcmidxs <- which( -1 == diff( c( idx1g[ 1 ] + 1, idx1g ) ) )
# find which elements of idx2 represent the end of a sequence
# to be replaced (meta-indexes)
destmidxs <- which( 1 == rev( diff( rev( c( idx1g, 0 ) ) ) ) )
# add counts from before NAs to destination rows
result <- dat
srcidxList <- vector( mode="list", length=length( destmidxs ) )
for ( i in seq.int( length( destmidxs ) ) ) {
# row to which data will be copied
destidx <- idx2[ destmidxs[ i ] ]
# sequence of indexes of source rows
srcidxss <- seq.int( from=idx1[ srcmidxs[ i ] ], to=destidx - 1 )
result[ destidx, "Count" ] <- ( dat[ destidx, "Count" ]
+ sum( dat[ srcidxss, "Count" ], na.rm=TRUE ) )
# keep a list of indexes to be removed
srcidxList[ i ] <- list( srcidxss )
}
# remove source rows
result <- result[ -unlist( srcidxList ), ]
result
}
On Sun, 20 Oct 2013, Jeff Newmiller wrote:
I thought this question looked interesting enough to make my own stab at it,
but in hindsight I think this business of combining the counts seems quite
unlikely to be necessary... it would be simpler and less damaging to the
original data pattern to just remove groups of rows having fewer than "N"
NAs.
removeNNAs <- function( dat, N, lessOrEqual=FALSE ) {
N1 <- N+1
rx <- rle( !is.na( dat$Count ) )
# indexes of the ends of each run of NAs or non-NAs
cs <- cumsum( rx$lengths )
# indexes of the ends of runs of NAs or non-NAs
cs2 <- cs[ !rx$values ]
# If the first Count is NA, then drop first run of NAs
if ( !rx$values[1] ) {
cs2 <- cs2[ -1 ]
}
# If the last Count is NA, then drop last run of NAs
if ( !rx$values[ length( rx$values ) ] ) {
cs2 <- cs2[ -length( cs2 ) ]
}
# cs2 is indexes of rows to potentially receive deleted Counts
# after collapse
cs2 <- cs2 + 1
# cs1 is indexes of non-NA Counts to be deleted
cs1 <- cs[ rx$values ][ seq.int( length( cs2 ) ) ]
# identify the indexes of the Count values before the strings
# of NAs that meet the criteria
if ( lessOrEqual ) {
idx0 <- N1 >= ( cs2 - cs1 )
} else {
idx0 <- N1 == ( cs2 - cs1 )
}
idx1 <- cs1[ idx0 ]
# identify the indexes of the Count values after the strings of
# NAs that meet the criteria
idx2 <- cs2[ idx0 ]
# Identify which indexes are both sources and destinations
idx1c <-c( idx2[ -length( idx2 ) ] == idx1[ -1 ], FALSE )
# identify groups of indexes that need to be merged
idx1g <- rev( cumsum( rev( !idx1c ) ) )
# find which elements of idx1 represent the beginning of a
# sequence of indexes to be replaced (meta-indexes)
srcmidxs <- which( -1 == diff( c( idx1g[ 1 ] + 1, idx1g ) ) )
# find which elements of idx2 represent the end of a sequence
# to be replaced (meta-indexes)
destmidxs <- which( 1 == rev( diff( rev( c( idx1g, 0 ) ) ) ) )
# add counts from before NAs to destination rows
result <- dat
srcidxList <- vector( mode="list", length=length( destmidxs ) )
for ( i in seq.int( length( destmidxs ) ) ) {
# row to which data will be copied
destidx <- idx2[ destmidxs[ i ] ]
# sequence of indexes of source rows
srcidxss <- seq.int( from=idx1[ srcmidxs[ i ] ], to=destidx - 1 )
result[ destidx, "Count" ] <- ( dat[ destidx, "Count" ]
+ sum( dat[ srcidxss, "Count" ], na.rm=TRUE )
# keep a list of indexes to be removed
srcidxList[ i ] <- list( srcidxss )
}
# remove source rows
result <- result[ -unlist( srcidxList ), ]
result
}
On Fri, 18 Oct 2013, arun wrote:
Hi,
Found a bug in the function when tested. So, try thi
[R] nlminb() - how do I constrain the parameter vector properly?
Greets,
I'm trying to use nlminb() to estimate the parameters of a bivariate normal
sample and during one of the iterations it passes a parameter vector to the
likelihood function resulting in an invalid covariance matrix that causes
dmvnorm() to throw an error. Thus, it seems I need to somehow communicate to
nlminb() that the final three parameters in my parameter vector are used to
construct a positive semi-definite matrix, but I can't see how to achieve this
using the constraint mechanism provided. Additional details are provided in the
code below.
Suggestions?
Best Regards,
Steven
Generate the data set I'm using:
mu<-c(1,5)
sigma<-c(1,2,2,6)
dim(sigma)<-c(2,2)
set.seed(83165026)
sample.full<-sample.deleted<-mvrnorm(50,mu,sigma)
delete.one<-c(1,5,13,20)
delete.two<-c(3,11,17,31,40,41)
sample.deleted[delete.one,1]<-NA
sample.deleted[delete.two,2]<-NA
missing<-c(delete.one,delete.two)
complete<-sample.deleted[!(is.na(sample.deleted[,1]) |
is.na(sample.deleted[,2])),]
deleted<-sample.deleted[missing,]
Try to estimate the parameters of the data set less the deleted values:
exact<-function(theta,complete,deleted){
one.only<-deleted[!(is.na(deleted[,1])),1]
two.only<-deleted[!(is.na(deleted[,2])),2]
u<-c(theta[1],theta[2])
sigma<-c(theta[3],theta[5],theta[5],theta[4])
dim(sigma)<-c(2,2)
-sum(log(dmvnorm(x=complete,mean=u,sigma=sigma)))-
sum(log(dnorm(one.only,u[1],sigma[1,1])))-
sum(log(dnorm(two.only,u[2],sigma[2,2])))
}
nlminb(start=c(0,0,1,1,0),objective=exact,complete=complete,deleted=deleted,control=list(trace=1))
Escape and it stops at Iteration 9 on my machine.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] bold dot size and name in plot
Dear Arun,Thankyou very much. Glad that you helped.:Deliza
> Date: Sat, 19 Oct 2013 21:31:09 -0700
> From: smartpink...@yahoo.com
> Subject: Re: [R] bold dot size and name in plot
> To: r-help@r-project.org
> CC: eliza_bo...@hotmail.com; j...@bitwrit.com.au
>
> Hi,
> May be this also helps. Using ?ggplot()
> indx <- 1 + 2*((z[,1]>80 & z[,1] < 100) & (z[,2]>2 & z[,2] < 4)) +
> 4*((z[,1]>40 & z[,1] < 60) & (z[,2]>8 & z[,2] < 10)) + 8*((z[,1]>100 & z[,1]
> < 120) & (z[,2]> 6 & z[,2] < 8))
> indx[indx>1] <- 2
> z1<-cbind(z,indx)
> z2 <- as.data.frame(z1)
> colnames(z2)[1:2] <- c("x","y")
> rect1 <-
> data.frame(xmin=c(80,40,100),xmax=c(100,60,120),ymin=c(2,8,6),ymax=c(4,10,8))
> library(ggplot2)
> p <- ggplot(x=x, y=y,data=z2)+scale_size(range=c(2,4))+
> geom_rect(data=rect1, aes(xmin=xmin, xmax=xmax, ymin=ymin,
> ymax=ymax),fill='gray80', alpha=0.8) +
> geom_point(data=z2,aes(x=x,y=y,size=indx,colour=factor(indx)))+scale_colour_manual(values=c("black","red"))+theme_bw()
> +theme(legend.position = "none")
> p+ annotate("text", label = "A", x = 90, y = 3) +
> annotate("text", label = "B", x = 50, y = 9) +
> annotate("text", label = "C", x = 110, y = 7)
>
>
> #or
> p+ geom_text(aes(90, 3, label="A")) +
> geom_text(aes(50, 9, label="B")) +
> geom_text(aes(110, 7, label="C"))
>
> A.K.
>
>
>
>
> On Saturday, October 19, 2013 5:30 PM, eliza botto
> wrote:
> Dear useRs,
>
> I have the following data "z" of two variables "x"(z[,1]) and "y"(z[,2]).
>
> > dput(z)
>
> structure(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
> 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37,
> 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56,
> 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
> 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94,
> 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110,
> 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124,
> 4.36222673881464, 3.69969030080286, 4.2171385106796, 3.13216766340978,
> 4.48272959444564, 5.15902588302327, 5.48017158746089, 4.27937198622481,
> 6.87798628712867, 10.6244257045475, 4.28697206373731, 2.36185506591713,
> 3.810824793, 3.53659639173799, 3.55002661079475, 1.69933515030185,
> 3.05902585757021, 7.36381894775327, 7.40245656365444, 3.6848569064997,
> 6.08996204878611, 6.28326065486612, 3.94011667392145, 4.32959623927668,
> 4.53659549630782, 10.5471767576564, !
> 3.25988629354992, 4.09457312472741, 3.80883550148271, 3.31050022622877,
> 4.97298874506525, 4.44111508385058, 8.07706784555988, 4.82859189180109,
> 4.28049047919568, 4.02054752142086, 4.45446182299703, 2.88902557744596,
> 4.17698124193575, 6.31474248231764, 2.68657606391183, 3.75708892669439,
> 4.70294050809375, 6.35750639997654, 5.66591712203293, 5.1676226656419,
> 6.36111068522507, 5.0746871785069, 3.74549388607439, 10.0499160660981,
> 2.73392821921647, 7.5468698295455, 3.63310870671515, 4.73756281172226,
> 3.12831980852897, 5.92478381773712, 3.97278768529811, 4.43603447353224,
> 3.22165005512334, 2.64048226457915, 4.82170867550029, 3.88802663642624,
> 6.28071352879497, 7.1898792291212, 3.45069580057411, 2.11874855961612,
> 6.83308998981622, 4.88594913751698, 7.2043715042717, 5.92944522632847,
> 5.28125713924232, 2.92406012836832, 5.51193983406794, 3.92034936881849,
> 4.02863562948636, 2.63086007943162, 2.29816871384959, 5.43277606496536,
> 2.748057506715, 2.98832890910102,
> 3.86750736412989, 10.02!
> 51585616053, 2.4508507570229, 5.15539283003769, 4.92251312085952, 3.75
> 60826804466, 3.5381750724, 3.80501827417895, 6.9500408077286,
> 5.63753820219307, 4.34417512875296, 2.38134080759689, 6.02021367624557,
> 7.61923679606813, 4.2639565075926, 5.19690782449524, 2.77914619706371,
> 2.80469912067042, 2.5030368187632, 3.66180330645551, 4.73586073723542,
> 4.03431276216711, 2.65910091882713, 2.83278553692168, 7.41646908081112,
> 1.87620181096256, 1.51781201506458, 1.92751580314289, 2.43506169436137,
> 1.88641686599088, 4.49013870147502, 2.15361460668737, 3.46482799242971,
> 2.69942251292495, 5.49300237653247, 3.34643344573202, 4.094064955793,
> 8.96571331963092, 2.22907411871184, 4.38143383339268, 6.52622723104758,
> 3.96979373172918, 5.52311167411852, 7.12115478587954), .Dim = c(124L, 2L),
> .Dimnames = list(c("1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2",
> "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2",
> "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2",
> "1_2", "1_2", "1_2",
> "1_2", "1_2", "1_2", "1_2"!
> , "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2",
> "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2",
> "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2",
> "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2", "1_2",
> "1_2"
Re: [R] Loop for taking sum of rows based on proximity to other non-NA rows
I thought this question looked interesting enough to make my own stab at
it, but in hindsight I think this business of combining the counts seems
quite unlikely to be necessary... it would be simpler and less damaging to
the original data pattern to just remove groups of rows having fewer than
"N" NAs.
removeNNAs <- function( dat, N, lessOrEqual=FALSE ) {
N1 <- N+1
rx <- rle( !is.na( dat$Count ) )
# indexes of the ends of each run of NAs or non-NAs
cs <- cumsum( rx$lengths )
# indexes of the ends of runs of NAs or non-NAs
cs2 <- cs[ !rx$values ]
# If the first Count is NA, then drop first run of NAs
if ( !rx$values[1] ) {
cs2 <- cs2[ -1 ]
}
# If the last Count is NA, then drop last run of NAs
if ( !rx$values[ length( rx$values ) ] ) {
cs2 <- cs2[ -length( cs2 ) ]
}
# cs2 is indexes of rows to potentially receive deleted Counts
# after collapse
cs2 <- cs2 + 1
# cs1 is indexes of non-NA Counts to be deleted
cs1 <- cs[ rx$values ][ seq.int( length( cs2 ) ) ]
# identify the indexes of the Count values before the strings
# of NAs that meet the criteria
if ( lessOrEqual ) {
idx0 <- N1 >= ( cs2 - cs1 )
} else {
idx0 <- N1 == ( cs2 - cs1 )
}
idx1 <- cs1[ idx0 ]
# identify the indexes of the Count values after the strings of
# NAs that meet the criteria
idx2 <- cs2[ idx0 ]
# Identify which indexes are both sources and destinations
idx1c <-c( idx2[ -length( idx2 ) ] == idx1[ -1 ], FALSE )
# identify groups of indexes that need to be merged
idx1g <- rev( cumsum( rev( !idx1c ) ) )
# find which elements of idx1 represent the beginning of a
# sequence of indexes to be replaced (meta-indexes)
srcmidxs <- which( -1 == diff( c( idx1g[ 1 ] + 1, idx1g ) ) )
# find which elements of idx2 represent the end of a sequence
# to be replaced (meta-indexes)
destmidxs <- which( 1 == rev( diff( rev( c( idx1g, 0 ) ) ) ) )
# add counts from before NAs to destination rows
result <- dat
srcidxList <- vector( mode="list", length=length( destmidxs ) )
for ( i in seq.int( length( destmidxs ) ) ) {
# row to which data will be copied
destidx <- idx2[ destmidxs[ i ] ]
# sequence of indexes of source rows
srcidxss <- seq.int( from=idx1[ srcmidxs[ i ] ], to=destidx - 1 )
result[ destidx, "Count" ] <- ( dat[ destidx, "Count" ]
+ sum( dat[ srcidxss, "Count" ], na.rm=TRUE )
# keep a list of indexes to be removed
srcidxList[ i ] <- list( srcidxss )
}
# remove source rows
result <- result[ -unlist( srcidxList ), ]
result
}
On Fri, 18 Oct 2013, arun wrote:
Hi,
Found a bug in the function when tested. So, try this (added one more line):
#Modified function
fun1 <- function(dat,n) {
rl <- rle(is.na(dat[,"Count"]))
indx <-
which(is.na(dat[,"Count"]))[rep(rl$lengths[rl$values],rl$lengths[rl$values])==n]
lst1 <- lapply(split(indx,((seq_along(indx)-1)%/%n)+1),function(x) {
x1 <- dat[c(min(x)-1L,x,max(x)+1L),]
x2 <- x1[!is.na(x1$Count),]
datN <-
data.frame(Position=max(x2$Position),Count=sum(x2$Count))
rowN <- row.names(x2)[x2$Position %in% max(x2$Position)]
row.names(datN) <- if(length(rowN)>1) rowN[1] else rowN
datN
})
names(lst1) <- NULL
lst1 <- lst1[!duplicated(sapply(lst1,row.names))] ##added
dat2 <- do.call(rbind,lst1)
indx2 <-
sort(unlist(lapply(split(indx,((seq_along(indx)-1)%/%n)+1),function(x)
c(min(x)-1L,x,c(max(x)+1L))),use.names=FALSE))
dat1New <- dat[-indx2[!indx2 %in% row.names(dat2)],]
dat1New[match(row.names(dat2),row.names(dat1New)),] <- dat2
row.names(dat1New) <- 1:nrow(dat1New)
dat1New
}
#Another function
fun2 <- function(dat,n){
indx <- cumsum(c(1,abs(diff(is.na(dat[,"Count"])
indx1 <- indx[is.na(dat[,"Count"])]
names(indx1) <- which(is.na(dat[,"Count"]))
indx2 <- indx1[indx1 %in% names(table(indx1))[table(indx1)==n]]
lst1 <- tapply(seq_along(indx2),list(indx2),FUN=function(i) {
x1 <- indx2[i]
x2 <- as.numeric(names(x1))
x3 <- dat[c(min(x2)-1L,x2,max(x2)+1L),]
x4 <- subset(x3, !is.na(Count))
x5 <-
data.frame(Position=max(x4$Position),Count=sum(x4$Count))
ind <- x4$Position %in% max(x4$Position)
row.names(x5) <- if(sum(ind)>1)
row.names(x4)[ind][1] else row.names(x4)[ind]
x5
})
attr(lst1,"dimnames") <- NULL
dat2 <- do.call(rbind,lst1)
indx3 <- sort(unlist(tapply(seq_along(indx2),list(indx2),FUN=function(i) {x1 <-
indx2[i]
x2 <- as.numeric(names(x1))
c(min(x2)-1L, x2,
max(x2)+1L)}),use.names=FALSE))
dat$id <- 1:nrow(dat)
dat2$id <- as.numeric(row.names(dat2))
library(plyr)
res
Re: [R] MaxLik estimation issues
Dear Filipe
On 14 October 2013 17:42, Filipe Ribeiro wrote:
> Dear Arne,
>
> First of all, thank you very much for your reply. Secondly, please accept my
> apologies for only give you an answer now, but I was moving from one country
> to another and only today I was able to get back to work.
> I did not write the model specification because it was just a trial, but my
> main idea was to apply the "Nelder-Mead" model.
> Since the beginning my guess was that something is wrong with the
> log-likelihood function, however, I don't find any problem with the same
> exact function applying the "optim" function:
>
> loglike.GGompi <- function(theta,age,deaths,exposures) {
> alpha <- exp(theta[1])
> beta <- exp(theta[2])
> gamma <- exp(theta[3])
> first <- alpha*exp(beta*age)
> second <- 1+(((alpha*gamma)/beta)*(exp(beta*age)-1))
> mu <- first/second
> llk <- -sum((deaths * log(mu)) + (- mu*exposures))
> return(llk)
> }
>
>
> fit1 <- optim(par=c(-4.1402817, -0.6375773, -1.6945914),
> loglike.GGompi,
> age=0:6,
> deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
> exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
> 231951.64, 69502.65, 15798.72))
>
> exp(fit1$par)
>
>> exp(fit1$par)
> [1] 0.01585683 0.53471945 0.25368426
>
> Due to this results I can't understand why...
Please note that optim() is (by default) *minimizing* the provided
function, while maxLik() is *maximizing* the provided log-likelihood
function.
Cheers,
Arne
> A 26/09/2013, às 05:28, Arne Henningsen escreveu:
>
> Dear Filipe
>
> On 25 September 2013 14:23, Filipe Ribeiro wrote:
>
> Hello everybody!
>
>
> I'm having some trouble to compute maximum likelihood
>
> estimations using maxLik package and I hope that you
>
> could give me a hint.
>
> The main problem is that I'm not able to get a result not
>
> even close to the ones given by glm() directly, and the
>
> second one is: "Error in maxNRCompute(fn = logLikAttr,
>
> fnOrig = fn, gradOrig = grad, hessOrig = hess, : NA in
>
> gradient".
>
>
> The codes:
>
> loglike.GGompiMaxLik <- function(theta,age,deaths,exposures) {
>
> alpha <- exp(theta[1])
>
> beta <- exp(theta[2])
>
> gamma <- exp(theta[3])
>
> first <- alpha*exp(beta*age)
>
> second <- 1+(((alpha*gamma)/beta)*(exp(beta*age)-1))
>
> mu <- first/second
>
> llk <- -sum((deaths * log(mu)) + (- mu*exposures))
>
> return(llk)
>
> }
>
>
>
> fit1 <- maxLik(loglike.GGompiMaxLik,
>
> age=0:6,
>
> deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
>
> exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
>
> 231951.64, 69502.65, 15798.72),
>
> start=c(-4.1402817, -0.6375773, -1.6945914))
>
>
> Do you know how I can solve this problem?
>
>
> You did not write which model specification you want to estimate but I
> am pretty sure that something in your log-likelihood function is
> incorrect. The log-likelihood value at the starting values of the
> parameters is so large that R even cannot calculate the likelihood
> value:
>
> a <- loglike.GGompiMaxLik(c(-4.1402817, -0.6375773, -1.6945914), age=0:6,
>
> + deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
> + exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
> + 231951.64, 69502.65, 15798.72))
>
> a
>
> [1] 580365.2
>
> exp(a)
>
> [1] Inf
>
> In the second iteration, the first parameter gets so small (large in
> absolute terms, -5e+10) that the log-likelihood value become extremely
> (numerically infinitely) large and the gradients cannot be computed
> (by the finite-difference method):
>
> fit1 <- maxLik(loglike.GGompiMaxLik,
>
> + age=0:6,
> + deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
> + exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
> + 231951.64, 69502.65, 15798.72),
> + start=c(-4.1402817, -0.6375773, -1.6945914))
> Iteration 2
> Parameter:
> [1] -5.174233e+10 -3.839076e+02 5.988668e+00
> Gradient:
> [,1] [,2] [,3]
> [1,] NaN NaN NaN
> Error in maxNRCompute(fn = logLikAttr, fnOrig = fn, gradOrig = grad,
> hessOrig = hess, :
> NA in gradient
>
> b <- loglike.GGompiMaxLik(c(-5.174233e+10, -3.839076e+02, 5.988668e+00),
> age=0:6,
>
> + deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
> + exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
> + 231951.64, 69502.65, 15798.72))
>
> b
>
> [1] Inf
>
> Please note that you can also find hints and ask questions about the
> maxLik package in the forums at maxLik's R-Forge site:
>
> https://r-forge.r-project.org/projects/maxlik/
>
> ...and please do not forget to cite the maxLik package in your publications:
>
> http://cran.r-project.org/web/packages/maxLik/citation.html
>
> Best wishes,
> Arne
--
Arne Henningsen
http://www.arne-henningsen.name
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-he
Re: [R] loop not working my way
Hello,
Make a _simple_ example, I don't see what packages like knitr or ggplot2
have anything to do with your problem.
Like this is, I think you're asking too much from r-help.
Rui Barradas
Em 19-10-2013 23:38, Laz escreveu:
Dear R users,
Dear R users,
(I had not included two more functions in the previous mail. This
version is complete)
There is a small problem which I don't know how to sort it out, based on
the former example I had explained earlier own.
I am calling my own functions which are based on simulations as below:
library(gmp)
library(knitr) # load this packages for publishing results
library(matlab)
library(Matrix)
library(psych)
library(foreach)
library(epicalc)
library(ggplot2)
library(xtable)
library(gdata)
library(gplots)
# function to calculate heritability
herit<-function(varG,varR=1)
{
h<-4*varG/(varG+varR)
h
}
h<-herit(0.081,1);h
###
# function to calculate random error
varR<-function(varG,h2)
{
varR<- varG*(4-h2)/h2
varR
}
#system.time(h<-varR(0.081,0.3));h
##
# function to calculate treatment variance
varG<-function(varR=1,h2)
{
varG<-varR*h2/(4-h2)
varG
}
# system.time(h<-varG(1,0.3));h
###
# calculating R inverse from spatial data
rspat<-function(rhox=0.6,rhoy=0.6)
{
s2e<-1
R<-s2e*eye(N)
for(i in 1:N) {
for (j in i:N){
y1<-y[i]
y2<-y[j]
x1<-x[i]
x2<-x[j]
R[i,j]<-s2e*(rhox^abs(x2-x1))*(rhoy^abs(y2-y1)) # Core AR(1)*AR(1)
R[j,i]<-R[i,j]
}
}
IR<-solve(R)
IR
}
### a function to generate A sparse matrix from a pedigree
ZGped<-function(ped)
{
ped2<-data.frame(ped)
lenp2<-length(unique(ped2$V1));lenp2 # how many Genotypes in total in
the pedigree =40
ln2<-length(g);ln2#ln2=nrow(matdf)=30
# calculate the new Z
Zped<-model.matrix(~ matdf$genotypes -1)# has order N*t = 180 by 30
dif<-(lenp2-ln2);dif # 40-30=10
#print(c(lenp2,ln2,dif))
zeromatrix<-zeros(nrow(matdf),dif);zeromatrix # 180 by 10
Z<-cbind(zeromatrix,Zped) # Design Matrix for random effect
(Genotypes): 180 by 40
# calculate the new G
M<-matrix(0,lenp2,lenp2) # 40 by 40
for (i in 1:nrow(ped2)) { M[ped2[i, 1], ped2[i, 2]] <- ped2[i, 3] }
G<-s2g*M # Genetic Variance covariance matrix for pedigree 2: 40 by 40
IG<-solve(G)
results<-c(IG, Z)
results
}
Three main functions here #
### function 1: generate a design (dataframe)
setup<-function(b,g,rb,cb,r,c,h2,rhox=0.6,rhoy=0.6,ped="F")
{
# where
# b = number of blocks
# t = number of treatments per block
# rb = number of rows per block
# cb = number of columns per block
# s2g = variance within genotypes
# h2 = heritability
# r = total number of rows for the layout
# c = total number of columns for the layout
### Check points
if(b==" ")
stop(paste(sQuote("block")," cannot be missing"))
if(!is.vector(g) | length(g)<3)
stop(paste(sQuote("treatments")," should be a vector and more than
2"))
if(!is.numeric(b))
stop(paste(sQuote("block"),"is not of class", sQuote("numeric")))
if(length(b)>1)
stop(paste(sQuote("block"),"has to be only 1 numeric value"))
if(!is.whole(b))
stop(paste(sQuote("block"),"has to be an", sQuote("integer")))
## Compatibility checks
if(rb*cb !=length(g))
stop(paste(sQuote("rb x cb")," should be equal to number of
treatment", sQuote("g")))
if(length(g) != rb*cb)
stop(paste(sQuote("the number of treatments"), "is not equal to",
sQuote("rb*cb")))
## Generate the design
g<<-g
genotypes<-times(b) %do% sample(g,length(g))
#genotypes<-rep(g,b)
block<-rep(1:b,each=length(g))
genotypes<-factor(genotypes)
block<-factor(block)
### generate the base design
k<-c/cb # number of blocks on the x-axis
x<<-rep(rep(1:r,each=cb),k) # X-coordinate
#w<-rb
l<-cb
p<-r/rb
m<-l+1
d<-l*b/p
y<<-c(rep(1:l,r),rep(m:d,r)) # Y-coordinate
## compact
matdf<<-data.frame(x,y,block,genotypes)
N<<-nrow(matdf)
mm<-summ(matdf)
ss<-des(matdf)
## Identity matrices
X<<-model.matrix(~block-1)
h2<<-h2;rhox<<-rhox;rhoy<<-rhoy
s2g<<-varG(varR=1,h2)
## calculate G and Z
ifelse(ped == "F",
c(IG<<-(1/s2g)*eye(length(g)),Z<<-model.matrix(~matdf$genotypes-1)),
c(IG<<- ZGped(ped)[[1]],Z<<-ZGped(ped)[[2]]))
## calculate R and IR
s2e<-1
ifelse(rhox==0 | rhoy==0, IR<<-(1/s2e)*eye(N),
IR<<-rspat(rhox=rhox,rhoy=rhoy))
C11<-t(X)%*%IR%*%X
C11inv<-solve(C11)
K<<-IR%*%X%*%C11inv%*%t(X)%*%IR
return(list( matdf= matdf,summary=mm,description=ss))
}
matrix0<-setup(b=4,g=seq(1,4,1),rb=2,cb=2,r=4,c=4,h2=0.3,rhox=0.6,rhoy=0.6,ped="F")[1]$matdf;
matrix0
x y block genotypes
1 1 1 1 1
2 1 2 1 3
3 2 1 1 2
4 2 2 1 4
5 3 1 2 1
6 3 2 2 3
7 4 1 2
Re: [R] Combining two models into one
You created a bunch of 'x' items but call for "X" in your formula. That's part of it. Next, you don't provide a "data=x" argument in your lm() call, so it goes off and looks for various capital-X items (I suspect). Hint: start with a small, simple dataset and work with it until you know how to use lm() and R's formula syntax. Only then go on to your full-on model. -- View this message in context: http://r.789695.n4.nabble.com/Combining-two-models-into-one-tp4678635p4678648.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] XML package not working
On 13-10-20 9:23 AM, Steven Dwayne Randolph wrote:
My apologies for not conforming to the posting guideline.
Sys.info()
sysname release
version
"Windows" "7 x64" "build 7601, Service Pack
1"
nodename machine
login
"xxNU247BZ1S" "x86-64"
"XX"
user effective_user
"xxx""xxx"
When I attempt to install a local copy of the xml.zip file:
in read.dcf(file.path(pkgname, "DESCRIPTION"), c("Package", "Type")) :
cannot open the connection
In addition: Warning messages:
1: In unzip(zipname, exdir = dest) : error 1 in extracting from zip file
2: In read.dcf(file.path(pkgname, "DESCRIPTION"), c("Package", "Type")) :
cannot open compressed file 'XML/DESCRIPTION', probable reason 'No such file
or directory'
I think it is pretty clear that the problem is at your end: you aren't
downloading the file properly, even though everyone else is. Perhaps
you are behind a firewall, or something else is interfering with your
downloads?
Duncan Murdoch
-Original Message-
From: Duncan Murdoch [mailto:murdoch.dun...@gmail.com]
Sent: Saturday, October 19, 2013 8:23 PM
To: Steven Dwayne Randolph; r-help@r-project.org
Cc: stevendrando...@aol.com
Subject: Re: [R] XML package not working
On 13-10-19 3:47 PM, Steven Dwayne Randolph wrote:
I know I cannot be the only one who is not able to install the XML package
from CRAN (zip or tar file) Many packages depend on this XML package. Can
someone help me either access the source for a good binary? I have received no
assistance from the author/developer of the package.
It installs fine for me.
Duncan Murdoch
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] saveXML() prefix argument
Thanks Earl and Milan.
Yes, the C code to serialize does branch and do things
differently for the different combinations of file, encoding and indent.
I have updated the code to use a different routine in libxml2 for this case
and that honors the indentation in this case. That will be in the next release
of XML.
In the meantime, you can use
cat( saveXML( doc, encoding = "UTF-8", indent = TRUE), file = "bob.xml")
rather than
saveXML(doc, file = "bob.xml", encoding = "UTF-8", indent = TRUE)
i.e. move the file argument to cat().
Thanks,
D.
On 10/19/13 4:36 AM, Milan Bouchet-Valat wrote:
> Le vendredi 18 octobre 2013 à 13:27 -0400, Earl Brown a écrit :
>> Thanks Duncan. However, now I can't get the Spanish and Portuguese accented
>> vowels to come out correctly and still keep the indents in the saved
>> document, even when I set encoding = "UTF-8":
>>
>> library("XML")
>> concepts <- c("español", "português")
>> info <- c("info about español", "info about português")
>>
>> doc <- newXMLDoc()
>> root <- newXMLNode("tips", doc = doc)
>> for (i in 1:length(concepts)) {
>> cur.concept <- concepts[i]
>> cur.info <- info[i]
>> cur.tip <- newXMLNode("tip", attrs = c(id = i), parent = root)
>> newXMLNode("h1", cur.concept, parent = cur.tip)
>> newXMLNode("p", cur.info, parent = cur.tip)
>> }
>>
>> # accented vowels don't come through correctly, but the indents are correct:
>> saveXML(doc, file = "test1.xml", indent = T)
>>
>> Resulting file looks like this:
>>
>>
>>
>> español
>> info about español
>>
>>
>> português
>> info about português
>>
>>
>>
>> # accented vowels are correct, but the indents are no longer correct:
>> saveXML(doc, file = "test2.xml", indent = T, encoding = "UTF-8")
>>
>> Resulting file:
>>
>> españolinfo about español> id="2">portuguêsinfo about português
>>
>> I tried to workaround the problem by simply loading in that resulting
>> file and saving it again:
>> doc2 <- xmlInternalTreeParse(file = "test2.xml", asTree = T)
>> saveXML(doc2, file = "test_word_around.xml", indent = T)
>>
>> but still don't get the indents.
>>
>> Does setting encoding = "UTF-8" override indents = TRUE in saveXML()?
> I can confirm the same issue happens here. What is interesting is that
> without the 'file' argument, the returned string includes the expected
> line breaks and spacing. These do not appear when redirecting the output
> to a file.
>
>> saveXML(doc, encoding="UTF-8", indent=T)
> [1] "\n\n \">\nespañol\ninfo about español\n \n
> \nportuguês\ninfo about
> português\n \n\n"
>
>> saveXML(doc, encoding="UTF-8", indent=T, file="test.xml")
>
> Contents of test.xml:
>
> españolinfo about español id="2">portuguêsinfo about português
>
>
>> sessionInfo()
> R version 3.0.1 (2013-05-16)
> Platform: x86_64-redhat-linux-gnu (64-bit)
>
> locale:
> [1] LC_CTYPE=fr_FR.utf8 LC_NUMERIC=C
> [3] LC_TIME=fr_FR.utf8LC_COLLATE=fr_FR.utf8
> [5] LC_MONETARY=fr_FR.utf8LC_MESSAGES=fr_FR.utf8
> [7] LC_PAPER=CLC_NAME=C
> [9] LC_ADDRESS=C LC_TELEPHONE=C
> [11] LC_MEASUREMENT=fr_FR.utf8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods
> base
>
> other attached packages:
> [1] XML_3.96-1.1
>
>
> Regards
>
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Nash-Sutcliffe coefficient
On 20-10-2013, at 17:29, Ahmed Attia wrote:
> Hi
>
> Could you please inform how can I get the Nash-Sutcliffe coefficient in r?
>
> Comparing observed vs. simulated values.
>
> Appreciated
>
Install package sos.
Then
findFn("Nash") should lead you to package hydroGOF which seems to have several
flavours.
Berend
> --
> Ahmed M. Attia
>
>
> Research Assistant
> Dept. Of Soil&Crop Sciences
> Texas A&M University
> ahmed.at...@ag.tamu.edu
> Cell phone: 001-979-248-5215
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Nash-Sutcliffe coefficient
Hi Could you please inform how can I get the Nash-Sutcliffe coefficient in r? Comparing observed vs. simulated values. Appreciated -- Ahmed M. Attia Research Assistant Dept. Of Soil&Crop Sciences Texas A&M University ahmed.at...@ag.tamu.edu Cell phone: 001-979-248-5215 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Errore : requires numeric/complex matrix/vector arguments
On 20-10-2013, at 16:16, valentina colombo wrote:
> Dear Mr. Hasselman,
> I have attached my code to solve (hopefully) my problem in Error-require
> numeric/complex matrix/vector
> Any suggestions?
> Thanks
>
Please answer to the list and NOT only privately. I'm forwarding this message
to R-help.
You are not showing how you constructed mX and others.
In the attachment you show how you read X.
I do not get your error when I generate your matrices with these commands:
set.seed(1)
N <- 126
mX <- matrix(runif(N*4),ncol=4)
mY <- matrix(runif(N*4),ncol=4)
mZ <- matrix(runif(N*4),ncol=4)
If you want to include data in a text-only mail use dput.
Berend
>
>
> > Subject: Re: [R] Errore : requires numeric/complex matrix/vector arguments
> > From: b...@xs4all.nl
> > Date: Sun, 20 Oct 2013 13:29:37 +0200
> > CC: r-help@r-project.org
> > To: valentina...@hotmail.it
> >
> >
> > On 20-10-2013, at 13:08, valentina colombo wrote:
> >
> > > Dear R users,I'm a new user of R. I'm trying to do a LM test an there is
> > > this type of error: Error in t(mX) %*% mX : requires numeric/complex
> > > matrix/vector arguments.
> > > To be clear I write down the code in which mY ( 126,1 ) mX (126,1)
> > > mZ(126,1) are matrix.
> > >
> > > LMTEST <- function(mY, mX, mZ)#mY, mX, mZ must be matrices!#returns the
> > > LM test statistic and the degree of freedom{iT = dim(mY)[1]ip =
> > > dim(mY)[2]iDF = dim(mZ)[2]*ipmE = mY - mX%*%solve(t(mX)%*%mX)%*%t(mX)%*%mY
> > > the error starts from the above step (t(mX)%*%mX)%*%t(mX)%*%mY
> > > RSS0 = t(mE)%*%mEmXX = cbind(mX, mZ)mK = mE -
> > > mXX%*%solve(t(mXX)%*%mXX)%*%t(mXX)%*%mERSS1 = t(mK)%*%mKdTR =
> > > sum(diag(solve(RSS0)%*%RSS1))LM = iT*(ip-dTR)pval =
> > > 1-pchisq(LM,df=iDF)return( c(pval, LM, iDF) )}
> > > Any suggestion? Where is the problem? I am getting craxy!
> >
> > Your code is a complete mess and thus unreadable because you posted in HTML.
> > Cleaning up and doing this
> >
> > LMTEST <- function(mY, mX, mZ)#mY, mX, mZ must be matrices!
> > #returns the LM test statistic and the degree of freedom
> > {iT = dim(mY)[1]
> > ip = dim(mY)[2]
> > iDF = dim(mZ)[2]*ip
> > mE = mY - mX%*%solve(t(mX)%*%mX)%*%t(mX)%*%mY
> > # the error starts from the above step
> > (t(mX)%*%mX)%*%t(mX)%*%mY
> > RSS0 = t(mE)%*%mE
> > mXX = cbind(mX, mZ)
> > mK = mE - mXX%*%solve(t(mXX)%*%mXX)%*%t(mXX)%*%mE
> > RSS1 = t(mK)%*%mK
> > dTR = sum(diag(solve(RSS0)%*%RSS1))
> > LM = iT*(ip-dTR)
> > pval = 1-pchisq(LM,df=iDF)
> > return( c(pval, LM, iDF) )
> > }
> >
> > set.seed(1)
> >
> > N <- 20
> > mX <- matrix(runif(N),ncol=1)
> > mY <- matrix(runif(N),ncol=1)
> > mZ <- matrix(runif(N),ncol=1)
> >
> > LMTEST(mY,mX,mZ)
> >
> > the answer I got was:
> >
> > [1] 0.004965514 7.891955826 1.0
> >
> >
> > So it must be your data.
> > Are you sure they are numeric? Have you checked with str(mX) etc?
> >
> > Berend
> >
> > > Valentina
> > > [[alternative HTML version deleted]]
> >
> >
> > Please don't post in html but in plain text.
> >
> >
> > >
> > > __
> > > R-help@r-project.org mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guide
> > > http://www.R-project.org/posting-guide.html
> > > and provide commented, minimal, self-contained, reproducible code.
> >
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Want to create a histogram
I don't think you do. I think that you are confusing a histogram and a bar
chart. With only five data points and five categories you have data for a
barplot.
I see Jim has provided one approach
Another, using the ggplot2 package (which you will probably have to install is :
# to install ggplot2 package
install.packages("ggplot2")
library(ggplot2)
dat1 <- data.frame( aa = c("very bad", "bad", "fair", "good", "very good"),
bb = c(159, 374, 3765, 11388, 6708))
ggplot(dat1, aes(aa, bb)) + geom_bar(stat="identity")
John Kane
Kingston ON Canada
> -Original Message-
> From: kan...@iitb.ac.in
> Sent: Sun, 20 Oct 2013 00:54:39 +
> To: r-help@r-project.org
> Subject: [R] Want to create a histogram
>
> Dear All,
>
> For a report that I am writing, I need to create a histogram plot with
> x-axis fixed as very bad, bad, fair, good and very good - i.e. the order
> not changed. Can someone give me an example? For sample purposes, I am
> giving the following data: 159, 374, 3765, 11388, 6708. Thanks,
>
> Kannan
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
FREE 3D MARINE AQUARIUM SCREENSAVER - Watch dolphins, sharks & orcas on your
desktop!
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Errore : requires numeric/complex matrix/vector arguments
On 20-10-2013, at 13:08, valentina colombo wrote:
> Dear R users,I'm a new user of R. I'm trying to do a LM test an there is this
> type of error: Error in t(mX) %*% mX : requires numeric/complex matrix/vector
> arguments.
> To be clear I write down the code in which mY ( 126,1 ) mX (126,1)
> mZ(126,1) are matrix.
>
> LMTEST <- function(mY, mX, mZ)#mY, mX, mZ must be matrices!#returns the LM
> test statistic and the degree of freedom{iT = dim(mY)[1]ip = dim(mY)[2]iDF =
> dim(mZ)[2]*ipmE = mY - mX%*%solve(t(mX)%*%mX)%*%t(mX)%*%mY
> the error starts from the above step (t(mX)%*%mX)%*%t(mX)%*%mY
> RSS0 = t(mE)%*%mEmXX = cbind(mX, mZ)mK = mE -
> mXX%*%solve(t(mXX)%*%mXX)%*%t(mXX)%*%mERSS1 = t(mK)%*%mKdTR =
> sum(diag(solve(RSS0)%*%RSS1))LM = iT*(ip-dTR)pval =
> 1-pchisq(LM,df=iDF)return( c(pval, LM, iDF) )}
> Any suggestion? Where is the problem? I am getting craxy!
Your code is a complete mess and thus unreadable because you posted in HTML.
Cleaning up and doing this
LMTEST <- function(mY, mX, mZ)#mY, mX, mZ must be matrices!
#returns the LM test statistic and the degree of freedom
{iT = dim(mY)[1]
ip = dim(mY)[2]
iDF = dim(mZ)[2]*ip
mE = mY - mX%*%solve(t(mX)%*%mX)%*%t(mX)%*%mY
# the error starts from the above step
(t(mX)%*%mX)%*%t(mX)%*%mY
RSS0 = t(mE)%*%mE
mXX = cbind(mX, mZ)
mK = mE - mXX%*%solve(t(mXX)%*%mXX)%*%t(mXX)%*%mE
RSS1 = t(mK)%*%mK
dTR = sum(diag(solve(RSS0)%*%RSS1))
LM = iT*(ip-dTR)
pval = 1-pchisq(LM,df=iDF)
return( c(pval, LM, iDF) )
}
set.seed(1)
N <- 20
mX <- matrix(runif(N),ncol=1)
mY <- matrix(runif(N),ncol=1)
mZ <- matrix(runif(N),ncol=1)
LMTEST(mY,mX,mZ)
the answer I got was:
[1] 0.004965514 7.891955826 1.0
So it must be your data.
Are you sure they are numeric? Have you checked with str(mX) etc?
Berend
> Valentina
> [[alternative HTML version deleted]]
Please don't post in html but in plain text.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Errore : requires numeric/complex matrix/vector arguments
It looks like you have posted in HTML and your code is basically unreadable.
Please repost in text format. Also I'd suggest reading one or both of the
following lists for some hints on how to create a good question.
I think a key point here is to supply sample data using the dput() function as
described in the links:
https://github.com/hadley/devtools/wiki/Reproducibility
http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example
Welcome to the R-help list
John Kane
Kingston ON Canada
> -Original Message-
> From: valentina...@hotmail.it
> Sent: Sun, 20 Oct 2013 11:08:08 +
> To: r-help@r-project.org
> Subject: [R] Errore : requires numeric/complex matrix/vector arguments
>
> Dear R users,I'm a new user of R. I'm trying to do a LM test an there is
> this type of error: Error in t(mX) %*% mX : requires numeric/complex
> matrix/vector arguments.
> To be clear I write down the code in which mY ( 126,1 ) mX (126,1)
> mZ(126,1) are matrix.
>
> LMTEST <- function(mY, mX, mZ)#mY, mX, mZ must be matrices!#returns the
> LM test statistic and the degree of freedom{iT = dim(mY)[1]ip =
> dim(mY)[2]iDF = dim(mZ)[2]*ipmE = mY -
> mX%*%solve(t(mX)%*%mX)%*%t(mX)%*%mY
> the error starts from the above step (t(mX)%*%mX)%*%t(mX)%*%mY
> RSS0 = t(mE)%*%mEmXX = cbind(mX, mZ)mK = mE -
> mXX%*%solve(t(mXX)%*%mXX)%*%t(mXX)%*%mERSS1 = t(mK)%*%mKdTR =
> sum(diag(solve(RSS0)%*%RSS1))LM = iT*(ip-dTR)pval =
> 1-pchisq(LM,df=iDF)return( c(pval, LM, iDF) )}
> Any suggestion? Where is the problem? I am getting craxy!
> Valentina
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
FREE 3D MARINE AQUARIUM SCREENSAVER - Watch dolphins, sharks & orcas on your
desktop!
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Errore : requires numeric/complex matrix/vector arguments
Dear R users,I'm a new user of R. I'm trying to do a LM test an there is this
type of error: Error in t(mX) %*% mX : requires numeric/complex matrix/vector
arguments.
To be clear I write down the code in which mY ( 126,1 ) mX (126,1)
mZ(126,1) are matrix.
LMTEST <- function(mY, mX, mZ)#mY, mX, mZ must be matrices!#returns the LM test
statistic and the degree of freedom{iT = dim(mY)[1]ip = dim(mY)[2]iDF =
dim(mZ)[2]*ipmE = mY - mX%*%solve(t(mX)%*%mX)%*%t(mX)%*%mY
the error starts from the above step (t(mX)%*%mX)%*%t(mX)%*%mY
RSS0 = t(mE)%*%mEmXX = cbind(mX, mZ)mK = mE -
mXX%*%solve(t(mXX)%*%mXX)%*%t(mXX)%*%mERSS1 = t(mK)%*%mKdTR =
sum(diag(solve(RSS0)%*%RSS1))LM = iT*(ip-dTR)pval = 1-pchisq(LM,df=iDF)return(
c(pval, LM, iDF) )}
Any suggestion? Where is the problem? I am getting craxy!
Valentina
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] R2admb compile problem
Good morning all
I am trying to use R2admb and have the following very simple example. The *.tpl
is included below. I am a Windows 7 user, using R3.01.
//Function1.tpl
//*
PARAMETER_SECTION
init_number x;
objective_function_value C;
PROCEDURE_SECTION
C= pow(x-5,2);
//*
I am running the example using three different methods
1. from DOS
2. From Emacs
3. From within R
1. from DOS
admb Function1
Function1
This runs correctly and works
2. From Emacs
I "Translate", "Build" and "Run" the file and it works correctly
3. From R
#Using the following code
require(R2admb)
fn="Function1"
setup_admb() #outputs the following [1] "c:\\ADMB\\admb101-gcc452-win64"
compile_admb(fn, verbose=T)
#I set 'verbose =T' to see whats happening
#last error line = 'collect2: ld returned 1 exit status'
#cant run this part of the code
run_admb(fn)
#Error in run_admb(fn) :
#executable Function1.exe not found: did you forget to compile it?
Note that if I first compile the tpl file using DOS and then run the following
from R, that everything works.
run_admb(fn) #to run the executable
results <- read_admb(fn)
results
clean_admb(fn)
DOES ANYONE KNOW WHY compile_admb(fn, verbose=T) DOES NOT WORK CORRECTLY?
//Function1.tpl
DATA_SECTION
PARAMETER_SECTION
init_number x;
objective_function_value C;
PRELIMINARY_CALCS_SECTION
PROCEDURE_SECTION
C= pow(x-5,2);
UNIVERSITY OF CAPE TOWN
This e-mail is subject to the UCT ICT policies and e-mai...{{dropped:9}}
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Combining two models into one
Hi guys, I'm trying to combining 2 models into 1 but maybe I'm not doing this in the right way. my dataset is about the level of the vitamin A in the blood measured at different times with two different treatments (ie. I have 6 coloums, vit A1 , time1, treat1, vit A2, time2, treat2, n=28). It wants to determinate if the variable treatment is really appropriate. E(y) = β0 + β1 x + β2 x^2 (treatment1) = γ0 + γ1 x + γ2 x^2 (treatment2) (level of vit a against time) satisfy the hypothesis H0 : β1 = γ1 and β2 = γ2 . The question is how can I combine the two models into one E(y) = X β with appropriate design matrix X and coefficient vector β? In this question β is 6x1. I did it but I think it is not good at all: vita<-cbind(dt$Vita..y., dt$Time..x., dt$Treatment,dt$Vita..y..1,dt$Time..x..1,dt$Treatment.1) x1 <- cbind(rep(1,28), matrix(0,28,2), vita[,1], matrix(0,28,2)) x11 <- cbind(rep(1,28), matrix(0,28,2), vita[,4], matrix(0,28,2)) x2 <- cbind(rep(0,28), rep(1,28), matrix(0,28,2), vita[,2], rep(0,28)) x22 <- cbind(rep(0,28), rep(1,28), matrix(0,28,2), vita[,5], rep(0,28)) x3 <- cbind(matrix(0,28,2), rep(1,28), matrix(0,28,2), vita[,3]) x33 <- cbind(matrix(0,28,2), rep(1,28), matrix(0,28,2), vita[,6]) x <- rbind(x1, x2, x3,x11, x22, x33) #create the design matrix ma <-c(vita[,1],vita[,4]) #response values mma <- data.frame(ma,x) #create the data frame for regression attach(mma) #attach the data frame to make it available to lm() ma.lm <- lm(ma ~ -1 + X1 + X2 + X3 + X4 + X5 + X6) #regr model without #intercept (note the -1 term) I am sure that is not good, because I should have x^2 and there isn't. Moreover it gives me an error: Error in model.frame.default(formula = ma ~ -1 + X1 + X2 + X3 + X4 + X5 + : variable lengths differ (found for 'X1') How can I do it? Tnx -- View this message in context: http://r.789695.n4.nabble.com/Combining-two-models-into-one-tp4678635.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] No P.values in polr summary
On Sat, 19 Oct 2013, Tal Galili wrote:
Vincent,
I believe Prof. Ripley is referring to this:
http://www.stats.ox.ac.uk/pub/MASS4/
Contact
Details:---
Contact me: tal.gal...@gmail.com |
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
www.r-statistics.com (English)
--
On Sat, Oct 19, 2013 at 3:22 PM, vincent guyader
wrote:
2013/10/19 Prof Brian Ripley
On 18/10/2013 15:01, Vincent Guyader wrote:
Hi everyone,
If I compute a "Ordered Logistic or Probit Regression" with the polr
function from MASS package. the summary give me : coefficients, Standard
error and Tvalue.. but not directly the p.value.
I can compute "manualy" the Pvalue, but Is there a way to directly
obtain
the pa.value, and I wonder why the p.valeu is not directly calculated,
is
there a reason?
How are you going to calculate the P values? Have you read the book for
which this is support software?: it explains why such Wald tests are
inappropriate and that the asymptotic theory can be wildly misleading.
Hi,
thanks for your answer.
to have the P.value I use this code :
http://www.ats.ucla.edu/stat/r/dae/ologit.htm
pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
A convenience option for computing this is the coeftest() method provided
in the "AER" package:
## data
library("foreign")
dat <- read.dta("http://www.ats.ucla.edu/stat/data/ologit.dta";)
## model
library("MASS")
m <- polr(apply ~ pared + public + gpa, data = dat, Hess = TRUE)
## coefficient test
library("AER")
coeftest(m)
Checking out the discussion in MASS (the book) on the usefulness of these
p-values is nevertheless a good idea, of course.
It give the same result as Stata, but you are right i'm not sure that it's
good. Could you please tell me which book you are talking about.
Regards
exemple :
house.plr <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data =
housing)
house.plr
summary(house.plr, digits = 3)
Regards
[[alternative HTML version deleted]]
__**
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/**listinfo/r-help<
https://stat.ethz.ch/mailman/listinfo/r-help>
PLEASE do read the posting guide http://www.R-project.org/**
posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
PLEASE do.
--
Brian D. Ripley, rip...@stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~**ripley/<
http://www.stats.ox.ac.uk/~ripley/>
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax: +44 1865 272595
__**
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/**listinfo/r-help<
https://stat.ethz.ch/mailman/listinfo/r-help>
PLEASE do read the posting guide http://www.R-project.org/**
posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
