[R] How to catch both warnings and errors?
I tried to use this solution (from over two years ago, but it remains an
official demo), but found that it only captured the last warning rather than
all of them. Instead, the code below collects all warnings.
tryCatch.W.E - function(expr)
{
W - list()
w.handler - function(w){ # warning handler
W - c(W,list(w))
invokeRestart(muffleWarning)
}
list(value = withCallingHandlers(tryCatch(expr, error = function(e) e),
warning = w.handler),
warning = W)
}
Sam Meyer
__
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] how to count A, C, T, G in each row in a big data.frame?
Hi Yao,
You could also use:
library(reshape2)
dd-dat1[,-(1:4)]
res-dcast(melt(within(dd,{id=row.names(dd)}),id.var=id),id~value,length)
head(res)
# id AA AG CC CT GA GG GT TC TG TT
#1 27412 29 10 0 0 13 1 0 0 0 0
#2 27413 0 0 4 9 0 0 0 12 0 28
#3 27414 0 0 0 0 0 0 0 0 0 53
#4 27415 0 0 53 0 0 0 0 0 0 0
#5 27416 0 0 3 9 0 0 0 12 0 29
#6 27417 0 0 53 0 0 0 0 0 0 0
#Just for comparison:
dat2- dat1[rep(row.names(dat1),2000),]
nrow(dat2)
#[1] 4
row.names(dat2)-1:4
dd - dat2[,-(1:4)]
system.time(res1- table(rownames(dd)[row(dd)], unlist(dd)))
# user system elapsed
# 5.840 0.104 5.954
system.time(res2 -
dcast(melt(within(dd,{id=row.names(dd)}),id.var=id),id~value,length))
# user system elapsed
# 3.100 0.064 3.167
head(res1,3)
# AA AG CC CT GA GG GT TC TG TT
# 1 29 10 0 0 13 1 0 0 0 0
# 10 0 4 0 0 6 43 0 0 0 0
# 100 19 15 0 0 15 4 0 0 0 0
head(res2,3)
# id AA AG CC CT GA GG GT TC TG TT
#1 1 29 10 0 0 13 1 0 0 0 0
#2 10 0 4 0 0 6 43 0 0 0 0
#3 100 19 15 0 0 15 4 0 0 0 0
A.K.
- Original Message -
From: Yao He yao.h.1...@gmail.com
To: R help r-help@r-project.org
Cc:
Sent: Wednesday, January 9, 2013 9:23 AM
Subject: [R] how to count A,C,T,G in each row in a big data.frame?
Dear All
I have a data.frame like that:
structure(list(name = c(Gga_rs10722041, Gga_rs10722249, Gga_rs10722565,
Gga_rs10723082, Gga_rs10723993, Gga_rs10724555, Gga_rs10726238,
Gga_rs10726461, Gga_rs10726774, Gga_rs10726967, Gga_rs10727581,
Gga_rs10728004, Gga_rs10728156, Gga_rs10728177, Gga_rs10728373,
Gga_rs10728585, Gga_rs10729598, Gga_rs10729643, Gga_rs10729685,
Gga_rs10729827), chr = c(7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L), pos = c(11248993L,
20038370L, 16164457L, 38050527L, 20307106L, 13707090L, 12230458L,
36732967L, 2790856L, 1305785L, 29631963L, 13606593L, 13656397L,
2261611L, 32096703L, 13733153L, 16524147L, 558735L, 12514023L,
3619538L), strand = c(+, +, +, +, +, +, +, +,
+, +, +, +, +, +, +, +, +, +, +, +),
X2353 = c(AA, TT, TT, CC, TT, CC, CC, TT,
CC, GG, AG, AG, AG, TT, CC, AG, CC, AA,
GG, GG), X2409 = c(AA, CT, TT, CC, CT, CC,
CC, TT, CC, GG, GG, AG, AG, TT, CC, AG,
CC, AA, AG, GA), X2500 = c(GA, TT, TT, CC,
TT, CC, CC, TT, CC, GG, GG, GG, GG, GT,
CT, GG, CC, AA, AA, AA), X2598 = c(AA, TT,
TT, CC, TT, CC, CC, TT, CC, GG, AA, AG,
GG, TT, CC, AG, TC, AA, AA, AG), X2610 = c(AA,
TT, TT, CC, TT, CC, CC, TT, CC, GG, GA,
GA, GG, TT, CC, GA, CC, AA, AA, GA), X2300 = c(GA,
TT, TT, CC, TT, CC, CC, TT, CC, GG, GA,
AA, AG, TT, TC, AA, TC, AA, AG, AA), X2507 = c(AG,
TT, TT, CC, TT, CC, CC, TT, CC, GG, GG,
GA, GG, TT, TC, GG, CC, AA, GA, AG), X2530 = c(AG,
TC, TT, CC, TC, CC, CC, TT, CC, GG, AA,
GG, GG, TT, CC, GG, CC, AA, AA, AA), X2327 = c(AA,
TT, TT, CC, TT, CC, CC, TT, CC, GG, GA,
GG, GG, TT, TC, GG, CC, AA, AA, AA), X2389 = c(AA,
CC, TT, CC, CC, CC, CC, TT, CC, AG, GG,
AG, GG, TT, TC, AG, CC, AA, AA, AA), X2408 = c(AA,
TT, TT, CC, TT, CC, CC, TT, CC, GG, GA,
GA, GG, TT, CC, GA, CC, AA, AA, AG), X2463 = c(AA,
TT, TT, CC, TT, CC, CC, TT, CC, GG, GG,
GG, GG, TT, CT, GG, CC, AA, AA, GA), X2420 = c(GA,
TC, TT, CC, TC, CC, CC, TT, CC, GG, AG,
GG, GG, TG, TT, GG, CT, AA, AA, AA), X2563 = c(GA,
CC, TT, CC, TC, CC, CC, TT, CC, GG, GA,
GG, GG, GT, TT, GG, CT, AA, AA, AA), X2462 = c(AA,
TT, TT, CC, TT, CC, CC, TT, CC, GG, AA,
GG, GG, GT, TC, GG, CC, AA, AA, AA), X2292 = c(GA,
TT, TT, CC, TT, CC, CC, TT, CC, GG, GA,
AA, GG, TG, TC, AA, TC, AA, AA, AA), X2405 = c(GA,
TT, TT, CC, TT, CC, CC, TT, CC, GG, GG,
AG, GG, TG, TT, AA, CT, AA, AA, AA), X2543 = c(AA,
TC, TT, CC, TC, CC, CC, TT, CC, GA, GA,
GA, GG, TT, CT, GA, TT, AA, AA, GG), X2557 = c(AG,
CT, TT, CC, CT, CC, CC, TT, CC, GG, AG,
GA, GG, GT, CT, GA, CT, AA, AA, AG), X2583 = c(GA,
CT, TT, CC, CT, CC, CC, TT, CC, GG, GA,
GG, GG, GG, CT, GA, CT, AA, AA, AG), X2322 = c(AG,
TT, TT, CC, TT, CC, CC, TT, CC, GG, GG,
GG, GG, GT, TT, GG, CC, AA, AA, GA), X2535 = c(AA,
TC, TT, CC, TT, CC, CC, TT, CC, GG, GA,
GG, GG, TT, CC, GG, CC, AA, AA, AG), X2536 = c(GA,
TC, TT, CC, TC, CC, CC, TT, CC, GG, GG,
AG, GG, TT, TC, AG, TC, AA, AA, GA), X2581 = c(AG,
CT, TT, CC, CT, CC, CC, TT, CC, GG, GG,
GA, GG, TT, CC, GA, CT, AA, AA, AG), X2570 = c(AA,
CT, TT, CC, CT, CC, CC, TT, CC, GG, GG,
GG, GG, TT, TC, GG, CC, AA, AA, GG), X2476 = c(AA,
TT, TT, CC, TT, CC, CC, TT, CC, GG, GG,
GG, GG, GT, TC, AG, CC, AA, AA, AG), X2534 = c(GA,
TC, TT, CC, TC, CC, CC, TT, CC, GG, GA,
AG, GG, TG, CC, AG, TC, AA, AA, AA), X2280 = c(AA,
TC, TT, CC, TC, CC, CC, TT, CC, GG, AG,
AG, GG, TT, CC, GG, CC, AA, AA, AG), X2316 = c(AA,
CC, TT, CC, CC, CC, CC, TT, CC, AG, AA,
AA, AG, TT, TC, GG, CT, AA, GG,
Re: [R] how to count A, C, T, G in each row in a big data.frame?
HI Yao,
You could use this: (Jorge's solution may be faster, I didn't check)
idx-sapply(strsplit(names(res),split=),anyDuplicated) #res from the previous
solution:
res1-do.call(cbind,lapply(LETTERS[c(1,3,7,20)],function(i){rowSums(data.frame(rowSums(res[idx==0][grep(i,names(res)[idx==0])]),2*res[idx!=0][grep(i,names(res)[idx!=0])]))}))
colnames(res1)-LETTERS[c(1,3,7,20)]
res2-data.frame(res,res1)
head(res2)
# id AA AG CC CT GA GG GT TC TG TT A C G T
#1 27412 29 10 0 0 13 1 0 0 0 0 81 0 25 0
#2 27413 0 0 4 9 0 0 0 12 0 28 0 29 0 77
#3 27414 0 0 0 0 0 0 0 0 0 53 0 0 0 106
#4 27415 0 0 53 0 0 0 0 0 0 0 0 106 0 0
#5 27416 0 0 3 9 0 0 0 12 0 29 0 27 0 79
#6 27417 0 0 53 0 0 0 0 0 0 0 0 106 0 0
A.K.
- Original Message -
From: Yao He yao.h.1...@gmail.com
To: arun smartpink...@yahoo.com
Cc: William Dunlap wdun...@tibco.com; R help r-help@r-project.org
Sent: Wednesday, January 9, 2013 11:46 PM
Subject: Re: [R] how to count A,C,T,G in each row in a big data.frame?
Hi arun
Then how could spilt them and get a table of letters count such as:
id AA AG CC CT GA GG GT TC TG TT
id A T C G
#1 27412 81 0 0 25
#2 27413 0 77 29 0
Thanks
2013/1/10 arun smartpink...@yahoo.com:
Hi Yao,
You could also use:
library(reshape2)
dd-dat1[,-(1:4)]
res-dcast(melt(within(dd,{id=row.names(dd)}),id.var=id),id~value,length)
head(res)
# id AA AG CC CT GA GG GT TC TG TT
#1 27412 29 10 0 0 13 1 0 0 0 0
#2 27413 0 0 4 9 0 0 0 12 0 28
#3 27414 0 0 0 0 0 0 0 0 0 53
#4 27415 0 0 53 0 0 0 0 0 0 0
#5 27416 0 0 3 9 0 0 0 12 0 29
#6 27417 0 0 53 0 0 0 0 0 0 0
#Just for comparison:
dat2- dat1[rep(row.names(dat1),2000),]
nrow(dat2)
#[1] 4
row.names(dat2)-1:4
dd - dat2[,-(1:4)]
system.time(res1- table(rownames(dd)[row(dd)], unlist(dd)))
# user system elapsed
# 5.840 0.104 5.954
system.time(res2 -
dcast(melt(within(dd,{id=row.names(dd)}),id.var=id),id~value,length))
# user system elapsed
# 3.100 0.064 3.167
head(res1,3)
# AA AG CC CT GA GG GT TC TG TT
# 1 29 10 0 0 13 1 0 0 0 0
# 10 0 4 0 0 6 43 0 0 0 0
# 100 19 15 0 0 15 4 0 0 0 0
head(res2,3)
# id AA AG CC CT GA GG GT TC TG TT
#1 1 29 10 0 0 13 1 0 0 0 0
#2 10 0 4 0 0 6 43 0 0 0 0
#3 100 19 15 0 0 15 4 0 0 0 0
A.K.
- Original Message -
From: Yao He yao.h.1...@gmail.com
To: R help r-help@r-project.org
Cc:
Sent: Wednesday, January 9, 2013 9:23 AM
Subject: [R] how to count A,C,T,G in each row in a big data.frame?
Dear All
I have a data.frame like that:
structure(list(name = c(Gga_rs10722041, Gga_rs10722249, Gga_rs10722565,
Gga_rs10723082, Gga_rs10723993, Gga_rs10724555, Gga_rs10726238,
Gga_rs10726461, Gga_rs10726774, Gga_rs10726967, Gga_rs10727581,
Gga_rs10728004, Gga_rs10728156, Gga_rs10728177, Gga_rs10728373,
Gga_rs10728585, Gga_rs10729598, Gga_rs10729643, Gga_rs10729685,
Gga_rs10729827), chr = c(7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L), pos = c(11248993L,
20038370L, 16164457L, 38050527L, 20307106L, 13707090L, 12230458L,
36732967L, 2790856L, 1305785L, 29631963L, 13606593L, 13656397L,
2261611L, 32096703L, 13733153L, 16524147L, 558735L, 12514023L,
3619538L), strand = c(+, +, +, +, +, +, +, +,
+, +, +, +, +, +, +, +, +, +, +, +),
X2353 = c(AA, TT, TT, CC, TT, CC, CC, TT,
CC, GG, AG, AG, AG, TT, CC, AG, CC, AA,
GG, GG), X2409 = c(AA, CT, TT, CC, CT, CC,
CC, TT, CC, GG, GG, AG, AG, TT, CC, AG,
CC, AA, AG, GA), X2500 = c(GA, TT, TT, CC,
TT, CC, CC, TT, CC, GG, GG, GG, GG, GT,
CT, GG, CC, AA, AA, AA), X2598 = c(AA, TT,
TT, CC, TT, CC, CC, TT, CC, GG, AA, AG,
GG, TT, CC, AG, TC, AA, AA, AG), X2610 = c(AA,
TT, TT, CC, TT, CC, CC, TT, CC, GG, GA,
GA, GG, TT, CC, GA, CC, AA, AA, GA), X2300 = c(GA,
TT, TT, CC, TT, CC, CC, TT, CC, GG, GA,
AA, AG, TT, TC, AA, TC, AA, AG, AA), X2507 = c(AG,
TT, TT, CC, TT, CC, CC, TT, CC, GG, GG,
GA, GG, TT, TC, GG, CC, AA, GA, AG), X2530 = c(AG,
TC, TT, CC, TC, CC, CC, TT, CC, GG, AA,
GG, GG, TT, CC, GG, CC, AA, AA, AA), X2327 = c(AA,
TT, TT, CC, TT, CC, CC, TT, CC, GG, GA,
GG, GG, TT, TC, GG, CC, AA, AA, AA), X2389 = c(AA,
CC, TT, CC, CC, CC, CC, TT, CC, AG, GG,
AG, GG, TT, TC, AG, CC, AA, AA, AA), X2408 = c(AA,
TT, TT, CC, TT, CC, CC, TT, CC, GG, GA,
GA, GG, TT, CC, GA, CC, AA, AA, AG), X2463 = c(AA,
TT, TT, CC, TT, CC, CC, TT, CC, GG, GG,
GG, GG, TT, CT, GG, CC, AA, AA, GA), X2420 = c(GA,
TC, TT, CC, TC, CC, CC, TT, CC, GG, AG,
GG, GG, TG, TT, GG, CT, AA, AA, AA), X2563 = c(GA,
CC, TT, CC, TC, CC, CC, TT, CC, GG, GA,
GG, GG, GT, TT, GG, CT, AA, AA, AA), X2462 = c(AA,
TT, TT, CC, TT, CC, CC, TT, CC, GG, AA,
GG, GG, GT, TC, GG, CC, AA, AA, AA), X2292
Re: [R] Using objects within functions in formulas
Thanks everyone, very helpful.
On 9 January 2013 18:33, David Winsemius dwinsem...@comcast.net wrote:
On Jan 9, 2013, at 8:53 AM, Aidan MacNamara wrote:
Dear all,
I'm looking to create a formula within a function to pass to glmer()
and I'm having a problem that the following example will illustrate:
library(lme4)
y1 = rnorm(10)
x1 = data.frame(x11=rnorm(10), x12=rnorm(10), x13=rnorm(10))
x1 = data.matrix(x1)
w1 = data.frame(w11=sample(1:3,10, replace=TRUE), w12=sample(1:3,10,
replace=TRUE), w13=sample(1:3,10, replace=TRUE))
test1 - function(x2, y2, w2) {
print(str(w2))
form = as.formula(paste(y2 ~ x2 + ,paste((1|w2$, names(w2),
),
collapse= + , sep=)))
m1 = glmer(form)
return(m1)
}
model1 = test1(x2=x1, y2=y1, w2=w1)
As can be seen from the print statement within the function, the
object w2 is present and is a data frame. However, the following
error occurs:
Error in is.factor(x) : object 'w2' not found
Generally regression functions in R will be expecting to get one 'data'
argument and build formulas using column names from that object.
test1 - function(x2, y2, w2) {
w3 - cbind(w2, x2, x2)
print(str(w3))
form = as.formula(paste(y2 ~ x2 + ,paste((1|, names(w2), ),
collapse= + , sep=)))
m1 = glmer(form, data=w3); print(summary(m1))
return(m1)
}
model1 = test1(x2=x1, y2=y1, w2=w1)
This can be rectified by making 'w2' global - defining it outside the
function. I know there are issues with defining formulas and
environment but I'm not sure why this problem is specific to 'w2' and
not the other objects passed to the function.
Any help would be appreciated.
Aidan MacNamara
EMBL-EBI
David Winsemius, MD
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] how to count A, C, T, G in each row in a big data.frame?
On 10 January 2013 01:04, Yao He yao.h.1...@gmail.com wrote: In fact I want to calculate the gene frequency of each SNP. Why don't you use bioconductor for your analysis instead of trying to develop by your own? For example: http://www.bioconductor.org/help/course-materials/2008/MGED08/BiostringsMGED2008.pdf For example alphabetFrequency function in the ShortRead package. I am sure you can handle I/O somehow to interface with available tools in bioconductor. __ 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] SRS, Stratified, and Cluster sampling
there isn't much, but Dr. Lumley's book is probably your best bet.. up till now, i think most people have learned survey methodology in the abstract and then applied that theory to their language of choice :) http://www.amazon.com/Complex-Surveys-Analysis-Survey-Methodology/dp/0470284307 On Wed, Jan 9, 2013 at 8:26 PM, David Arnold dwarnol...@suddenlink.netwrote: Hi, Has anyone done (or know of) any nice R activities that help introductory students ( and teachers :) ) better understand the concepts of simple vs stratified vs cluster sampling? Any links? David -- View this message in context: http://r.789695.n4.nabble.com/SRS-Stratified-and-Cluster-sampling-tp4655099.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. [[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] multiple versions of function
On Jan 9, 2013, at 9:00 PM, ivo welch ivo.we...@gmail.com wrote:
mea culpa.
f - function(...) {
## parse out the arguments and then do something with them
}
## all of these should result in the same actions
f(2,3) ## interprets a to be first and b to be second
f(a=2,b=3)
f(b=3,a=2)
These will all be the same automatically for functions with the signature f(a,
b, ...) [grammatical not variadic ellipsis there] -- basic call matching,
nothing fancy.
f(data.frame(a=2,b=3))
f(data.frame(b=3,a=1))
Perhaps test if your first arg is a df and, if so, loop over it row-wise
building the function calls with do.call() -- something like:
# Untested
f - function(a, b){
if(is.data.frame(a)) return(lapply(seq_len(NROW(a)), function(n) do.call(f,
a[n,]))
## regular function code here
}
You should probably also fiddle with Recall() to make the recursive structure a
little more robust.
Though it seems that generic functions make more sense here.
Michael Weylandt
On Tue, Jan 8, 2013 at 8:00 AM, David Winsemius dwinsem...@comcast.net
wrote:
On Jan 7, 2013, at 6:58 PM, ivo welch wrote:
hi david---can you give just a little more of an example? the
function should work with call by order, call by name, and data frame
whose columns are the names. /iaw
It is I who should be expecting you to provide an example.
-- David.
__
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] Levels in new data fed to SVM
On 08.01.2013 21:14, Claus O'Rourke wrote: Hi all, I've encountered an issue using svm (e1071) in the specific case of supplying new data which may not have the full range of levels that were present in the training data. I've constructed this really primitive example to illustrate the point: library(e1071) training.data - data.frame(x = c(yellow,red,yellow,red), a = c(alpha,alpha,beta,beta), b = c(a, b, a, c)) my.model - svm(x ~ .,data=training.data) test.data - data.frame(x = c(yellow,red), a = c(alpha,beta), b = c(a, b)) predict(my.model,test.data) Error in predict.svm(my.model, test.data) : test data does not match model ! levels(test.data$b) - levels(training.data$b) predict(my.model,test.data) 1 2 yellowred Levels: red yellow In the first case test.data$b does not have the level c and this results in the input data being rejected. I've debugged this down to the point of model matrix creation in the SVM R code. Once I fill up the levels in the test data with the levels from the original data, then there is no problem at all. Assuming my test data has to come from another source where the number of category levels seen might not always be as large as those for the original training data, is there a better way I should be handling this? You have to tell the factor about the possible levels, it does not necessarily contain examples. That means: levels(test.data$b) - C(a, b, c) predict(my.model,test.data) will help. Best, Uwe Ligges Thanks __ 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] ./R: error while loading shared libraries
I'd move this to the R-SIG-Fedora list and, in doing so, give more info about your install process: built yourself, package manager, etc. MW On Wed, Jan 9, 2013 at 7:31 PM, Adam Dahman adamo...@yahoo.ca wrote: Hi, I have installed R on linux using a non root account. I am getting this error when trying to use it : ./R: error while loading shared libraries: libRblas.so: cannot open shared object file: No such file or directory Linux version I am using : Linux version 2.6.32-131.17.1.el6.x86_64 (mockbu...@x86-007.build.bos.redhat.com) (gcc version 4.4.5 20110214 (Red Hat 4.4.5-6) (GCC) ) #1 SMP Thu Sep 29 10:24:25 EDT 2011 Can someone help ? Regards, Adam [[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] Levels in new data fed to SVM
Thanks for clarifying! On Thu, Jan 10, 2013 at 12:47 PM, Uwe Ligges lig...@statistik.tu-dortmund.de wrote: On 08.01.2013 21:14, Claus O'Rourke wrote: Hi all, I've encountered an issue using svm (e1071) in the specific case of supplying new data which may not have the full range of levels that were present in the training data. I've constructed this really primitive example to illustrate the point: library(e1071) training.data - data.frame(x = c(yellow,red,yellow,red), a = c(alpha,alpha,beta,beta), b = c(a, b, a, c)) my.model - svm(x ~ .,data=training.data) test.data - data.frame(x = c(yellow,red), a = c(alpha,beta), b = c(a, b)) predict(my.model,test.data) Error in predict.svm(my.model, test.data) : test data does not match model ! levels(test.data$b) - levels(training.data$b) predict(my.model,test.data) 1 2 yellowred Levels: red yellow In the first case test.data$b does not have the level c and this results in the input data being rejected. I've debugged this down to the point of model matrix creation in the SVM R code. Once I fill up the levels in the test data with the levels from the original data, then there is no problem at all. Assuming my test data has to come from another source where the number of category levels seen might not always be as large as those for the original training data, is there a better way I should be handling this? You have to tell the factor about the possible levels, it does not necessarily contain examples. That means: levels(test.data$b) - C(a, b, c) predict(my.model,test.data) will help. Best, Uwe Ligges Thanks __ 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] Precision of values 53 bits
Hi, I am working with large numbers and identified that R looses precision for such high numbers. The precision is lost exactly when the number is equal or larger than 53 bits. See the following output which shows that the numbers below 53 bit have proper precision: 2^53 [1] 9007199254740992 2^53-1 [1] 9007199254740991 2^53-2 [1] 9007199254740990 Now, see the numbers above 53 bit: 2^53 [1] 9007199254740992 2^53+1 [1] 9007199254740992 2^53+2 [1] 9007199254740994 2^53+3 [1] 9007199254740996 2^53+4 [1] 9007199254740996 Is there a solution to the problem? Thanks a lot Stephan __ 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] Cartesian co-ordinate generation
Hi. This is my first post to the list so apologies if it is not formatted correctly. I have a dataset from a forest survey with trees marked in an arbitrary co-ordinate system of the form: X from 1000 to 1100, and Y from 1000 to 1100. I have recently resurveyed the forest and found trees that I need to add. For the new trees I have recorded angle and distance to existing trees (3 per new tree). So my question is two-fold: 1. Is there a package that anyone can recommend to help me translate these bearings and distances into the existing cartesian co-ordinate system? 2. Is there a way that I can quantify the locational error, assuming that there is a triangular intersect? Any suggestions would be greatly appreciated. Thanks Austin [[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] Precision of values 53 bits
Hi, I am working with large numbers and identified that R looses precision for such high numbers. The precision is lost exactly when the number is equal or larger than 53 bits. See the following output which shows that the numbers below 53 bit have proper precision: 2^53 [1] 9007199254740992 2^53-1 [1] 9007199254740991 2^53-2 [1] 9007199254740990 Now, see the numbers above 53 bit: 2^53 [1] 9007199254740992 2^53+1 [1] 9007199254740992 2^53+2 [1] 9007199254740994 2^53+3 [1] 9007199254740996 2^53+4 [1] 9007199254740996 Is there a solution to the problem? Thanks a lot Stephan __ 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] Semi Parametric Bootstrap
Greetings to you all,
I am performing a semi parametric bootstrap in R on a Gamma Distributed
data and a Binomial distributed data. The main challenge am facing is
the fact that the residual variance depends on the mean (if I am correct).
I strongly feel that the script below may be wrong due to mean-variance
relationship
#R code###
fit1s -glm(mydata$vzv~mydata$age.c+mydata$age2+mydata$sex1,
family=Gamma(link=log))
x.betahat1-fit1s$fitted.values
res1-fit1s$residuals
b-1000
for (i in 1:b){
b.i - sample(index, size=n, replace=T)
res.star1=res1[b.i]
bst1=x.betahat1+res.star1
mydata1 -data.frame(age,age2,sex,bst1)
Modeling
fit11 -glm(bst1~age+age2+sex, family=Gamma(link=log),data=mydata1)
}
Can someone help me correct this code? Kindly advice on Binomial data
as well
Happy New year2013!
-- ___
Paul K. Musingila
__
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] r-help
-- ___ Paul K. Musingila University of Hasselt, I-Biostat, Diepenbeek, Belgium Mobile: +254-724-423532, +32-48-637-4558 E-mail: paul.musing...@student.uhasselt.be, pmusing...@gmail.com Skype: pmusingila When darkness overtakes the godly, light will come bursting in Psalm 112:4 __ 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] Semi Parametric Bootstrap
Greetings to you all,
I am performing a semi parametric bootstrap in R on a Gamma Distributed
data and a Binomial distributed data. The main challenge am facing is
the fact that the residual variance depends on the mean (if I am correct).
I strongly feel that the script below may be wrong due to mean-variance
relationship
#R code###
fit1s -glm(mydata$vzv~mydata$age.c+mydata$age2+mydata$sex1,
family=Gamma(link=log))
x.betahat1-fit1s$fitted.values
res1-fit1s$residuals
b-1000
for (i in 1:b){
b.i - sample(index, size=n, replace=T)
res.star1=res1[b.i]
bst1=x.betahat1+res.star1
mydata1 -data.frame(age,age2,sex,bst1)
Modeling
fit11 -glm(bst1~age+age2+sex, family=Gamma(link=log),data=mydata1)
}
Can someone help me correct this code? Kindly advice on Binomial data
as well
Happy New year2013!
-- ___
Paul K. Musingila
__
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] Find the functional relationship between two variables in R?
Hi, I have two variables x and y and the functional relationship between x and y is like: y=x^2+log(x). My question is that is it possible to apply some method to reconstruct the functional form based on the training data that is generated from it? I understand that there are many methods for regresstion but none of them can predict the functions that consist of operators. For example in R: x=seq(1:100) y=x^2+log(x) result=somemodel(x,y) summary(result) The output of the method is ideally a mathematical equation that fits to the training data, e.g. y=x^2+log(x) Any ideas will be appreciated! Best, Jing __ 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] how to count A, C, T, G in each row in a big data.frame?
Hi Yao,
Did comparison of the different methods:
dat2- dat1[rep(row.names(dat1),2000),]
nrow(dat2)
#[1] 4
row.names(dat2)-1:4
dd - dat2[,-(1:4)]
foo - function(x) table(factor(unlist(strsplit(as.character(x), )), levels
= c('A','C','G','T')))
system.time(res3-t(apply(dat2[, -c(1:4)], 1, foo)) )
# user system elapsed
# 10.212 0.032 10.265
library(reshape2)
system.time(res2 -
dcast(melt(within(dd,{id=row.names(dd)}),id.var=id),id~value,length))
#user system elapsed
# 3.232 0.096 3.358
system.time({idx-sapply(strsplit(names(res2),split=),anyDuplicated)
res4-do.call(cbind,lapply(LETTERS[c(1,3,7,20)],function(i){rowSums(data.frame(rowSums(res2[idx==0][grep(i,names(res2)[idx==0])]),2*res2[idx!=0][grep(i,names(res2)[idx!=0])]))}))
colnames(res4)-LETTERS[c(1,3,7,20)]})
# user system elapsed
# 0.168 0.000 0.171
#Combined result
system.time({res2 -
dcast(melt(within(dd,{id=row.names(dd)}),id.var=id),id~value,length)
idx-sapply(strsplit(names(res2),split=),anyDuplicated)
res4-do.call(cbind,lapply(LETTERS[c(1,3,7,20)],function(i){rowSums(data.frame(rowSums(res2[idx==0][grep(i,names(res2)[idx==0])]),2*res2[idx!=0][grep(i,names(res2)[idx!=0])]))}))
colnames(res4)-LETTERS[c(1,3,7,20)]
res5-data.frame(res2,res4)
})
# user system elapsed
# 3.025 0.040 3.072
head(res5,5)
# id AA AG CC CT GA GG GT TC TG TT A C G T
#1 1 29 10 0 0 13 1 0 0 0 0 81 0 25 0
#2 10 0 4 0 0 6 43 0 0 0 0 10 0 96 0
#3 100 19 15 0 0 15 4 0 0 0 0 68 0 38 0
#4 1000 19 15 0 0 15 4 0 0 0 0 68 0 38 0
#5 1 19 15 0 0 15 4 0 0 0 0 68 0 38 0
system.time(res6- table(rownames(dd)[row(dd)], unlist(dd)))
# user system elapsed
# 4.840 0.016 4.899
head(res6,5)
#
AA AG CC CT GA GG GT TC TG TT
# 1 29 10 0 0 13 1 0 0 0 0
# 10 0 4 0 0 6 43 0 0 0 0
# 100 19 15 0 0 15 4 0 0 0 0
# 1000 19 15 0 0 15 4 0 0 0 0
# 1 19 15 0 0 15 4 0 0 0 0
system.time({#rows would be:
test2-apply(dat2[,-c(1:4)],1,function(x){table(t(x))})
#find single values in a row
res7- sapply(test2,function(row){
allVars-paste(names(row),collapse=)
u - unique(strsplit(allVars,)[[1]])
parts-sapply(names(row),function(x){u%in%strsplit(x,)[[1]]})
mat-parts%*%row
rownames(mat)-u
mat
})})
#user system elapsed
#21.553 0.000 21.591
A.K.
- Original Message -
From: Yao He yao.h.1...@gmail.com
To: arun smartpink...@yahoo.com
Cc: William Dunlap wdun...@tibco.com; R help r-help@r-project.org
Sent: Wednesday, January 9, 2013 11:46 PM
Subject: Re: [R] how to count A,C,T,G in each row in a big data.frame?
Hi arun
Then how could spilt them and get a table of letters count such as:
id AA AG CC CT GA GG GT TC TG TT
id A T C G
#1 27412 81 0 0 25
#2 27413 0 77 29 0
Thanks
2013/1/10 arun smartpink...@yahoo.com:
Hi Yao,
You could also use:
library(reshape2)
dd-dat1[,-(1:4)]
res-dcast(melt(within(dd,{id=row.names(dd)}),id.var=id),id~value,length)
head(res)
# id AA AG CC CT GA GG GT TC TG TT
#1 27412 29 10 0 0 13 1 0 0 0 0
#2 27413 0 0 4 9 0 0 0 12 0 28
#3 27414 0 0 0 0 0 0 0 0 0 53
#4 27415 0 0 53 0 0 0 0 0 0 0
#5 27416 0 0 3 9 0 0 0 12 0 29
#6 27417 0 0 53 0 0 0 0 0 0 0
#Just for comparison:
dat2- dat1[rep(row.names(dat1),2000),]
nrow(dat2)
#[1] 4
row.names(dat2)-1:4
dd - dat2[,-(1:4)]
system.time(res1- table(rownames(dd)[row(dd)], unlist(dd)))
# user system elapsed
# 5.840 0.104 5.954
system.time(res2 -
dcast(melt(within(dd,{id=row.names(dd)}),id.var=id),id~value,length))
# user system elapsed
# 3.100 0.064 3.167
head(res1,3)
# AA AG CC CT GA GG GT TC TG TT
# 1 29 10 0 0 13 1 0 0 0 0
# 10 0 4 0 0 6 43 0 0 0 0
# 100 19 15 0 0 15 4 0 0 0 0
head(res2,3)
# id AA AG CC CT GA GG GT TC TG TT
#1 1 29 10 0 0 13 1 0 0 0 0
#2 10 0 4 0 0 6 43 0 0 0 0
#3 100 19 15 0 0 15 4 0 0 0 0
A.K.
- Original Message -
From: Yao He yao.h.1...@gmail.com
To: R help r-help@r-project.org
Cc:
Sent: Wednesday, January 9, 2013 9:23 AM
Subject: [R] how to count A,C,T,G in each row in a big data.frame?
Dear All
I have a data.frame like that:
structure(list(name = c(Gga_rs10722041, Gga_rs10722249, Gga_rs10722565,
Gga_rs10723082, Gga_rs10723993, Gga_rs10724555, Gga_rs10726238,
Gga_rs10726461, Gga_rs10726774, Gga_rs10726967, Gga_rs10727581,
Gga_rs10728004, Gga_rs10728156, Gga_rs10728177, Gga_rs10728373,
Gga_rs10728585, Gga_rs10729598, Gga_rs10729643, Gga_rs10729685,
Gga_rs10729827), chr = c(7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L), pos = c(11248993L,
20038370L, 16164457L, 38050527L, 20307106L, 13707090L, 12230458L,
36732967L, 2790856L, 1305785L, 29631963L, 13606593L, 13656397L,
2261611L, 32096703L, 13733153L, 16524147L, 558735L,
Re: [R] Precision of values 53 bits
Perhaps here?: https://r-forge.r-project.org/projects/rmpfr/ M On Thu, Jan 10, 2013 at 10:58 AM, Stephan Mueller stephan.muel...@atsec.com wrote: I am working with large numbers and identified that R looses precision for such high numbers. The precision is lost exactly when the number is equal or larger than 53 bits. See the following output which shows that the numbers below 53 bit have proper precision: __ 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] mgcv: Plotting probabilities for binomial GAM with crossed random intercepts and factor by variable
mgcv: Constructing probabilities for binomial GAM with crossed random intercepts and factor by variable Hello, (I'm sorry if this has been discussed elsewhere; I may not have been looking in the right places.) I ran a binomial GAM in which Correct is modelled in terms of the participant's age and the modality in which the stimulus is presented (written vs spoken). Participants (Subject) and stimuli are also included as crossed random intercepts. age.gam - bam(Correct ~ Modality + s(Age, by=Modality) + s(Subject, bs=re) + s(Stimulus, bs=re), data = dat, family=binomial) Parametric coefficients: Estimate Std. Error z value Pr(|z|) (Intercept) -2.242 1.121 -2.000 0.0455 * ModalityWritten -1.043 1.263 -0.826 0.4087 Approximate significance of smooth terms: edf Ref.df Chi.sq p-value s(Age):ModalitySpoken 4.883 5.477 98.45 2e-16 *** s(Age):ModalityWritten 5.675 6.363 126.81 2e-16 *** s(Subject) 133.296 161.000 1472.35 2e-16 *** s(Stimulus) 83.376 85.000 5100.59 2e-16 *** I would now like to plot the predicted probabilities for Correct == yes, say for written stimuli: plot(age.gam, select=2, shift=(-2.242-1.043), trans=plogis) Though the overall shape of the curve is about what I expected based on a earlier visual exploration, the plotted probabilities are much too low (about 5% for values of Age where I'd expected 40%). Here are the raw counts: xtabs(~ Modality + Correct, dat) Correct Modality 01 Spoken 4146 2693 Written 4323 3011 Is it possible that I need to substract the values of s(Subject).1 and s(Stimulus).1 from coef(age.gam) to the trans argument, too? Those are 0.55 and -2.86, respectively, which would provide a much better match between the plotted and the expected probabilities. (But then what does the (Intercept) represent in this model?) Thanks, Jan __ 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] Precision of values 53 bits
I am working with large numbers and identified that R looses
precision for such high numbers.
Yes. R uses standard 32-bit double precision. See ?double in your R help
system. And welcome to finite precision arithmetic, which is a very widely
known issue in digital comuting ever since it was invented.
Is there a solution to the problem?
Yes, lots, said Bilbo, before he remembered not to give his friends away.
No, none at all, not one, he said immediately afterwards. [1]
R cannot easily be recompiled to use higher precision, so in that sense, 'none
at all'.
However, you could use something like the Rmpfr package for arbitrary precision
arithmetic.
On linux/unix you can use bc (see
http://r.789695.n4.nabble.com/Arbitrary-Precision-Numbers-td855931.html)
Or you could do basic things that address the issue: for example, scale,
mean-centre or transform the numbers or change your parameterisation so that
you do not need high numerical precision on large numbers.
Steve Ellison
[1] JRR Tolkien, ' The Hobbit', Chapter 3 (1937)
Thanks a lot
Stephan
__
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.
***
This email and any attachments are confidential. Any use...{{dropped:8}}
__
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] Titles - main and subtitle won't plot with errbar
Hi, I'm struggling with errbar graphics. I'm trying to plot an x-y graph with correct labelling, however can't seem to get main and sub to show on my graph. They do work when I use title(main=, etc, but this will make it look at lot messier,I'll have to blank out ylab= , and I need to try and get the titles to update automatically according to my excel column headings and paste function. Example code require(Hmisc) data1-array(sample(1:100,35),dim=c(5,7)) data1[,1]-1:5 sd=apply(data1[,2:7],1,sd) mean=rowMeans(data1[,2:7]) #the base plot works errbar(x=data1[,1],y=mean,yplus=mean+sd,yminus=mean-sd) #titles are shown correctly using plot plot(x=data1[,1],y=mean,ylab=value,xlab=time,main=Title,sub=subtitle) #the ylab and xlab update correctly, however main and sub don't? errbar(x=data1[,1],y=mean,yplus=mean+sd,yminus=mean-sd,ylab=value,xlab=time,main=Title,sub=subtitle) (my original is a lot more complicated as I'm reading from excel, and have managed to plot SD and mean for my dataset, but for some reason the main and sub commands aren't working in errbar!) Thanks, Laura -- View this message in context: http://r.789695.n4.nabble.com/Titles-main-and-subtitle-won-t-plot-with-errbar-tp4655149.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.
[R] Lambert W question
Dear All, I am using the following model equation: k*(lambertW_base(b=0,((a)/k)*exp(((a)-z*(t-t0))/k))) I would like to run this through OpenBUGS, but it does not recognize the lambert function. Would you have any thoughts on how to re-vrite this equation matemathically so that it could be run on OpenBUGS? apreciate the help, Sincerely, Andrs [[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] Precision of values 53 bits
On 13-01-10 6:01 AM, Stephan Mueller wrote: Hi, I am working with large numbers and identified that R looses precision for such high numbers. The precision is lost exactly when the number is equal or larger than 53 bits. See the following output which shows that the numbers below 53 bit have proper precision: 2^53 [1] 9007199254740992 2^53-1 [1] 9007199254740991 2^53-2 [1] 9007199254740990 Now, see the numbers above 53 bit: 2^53 [1] 9007199254740992 2^53+1 [1] 9007199254740992 2^53+2 [1] 9007199254740994 2^53+3 [1] 9007199254740996 2^53+4 [1] 9007199254740996 Is there a solution to the problem? R has no native numeric type with more than 53 bit precision, but other packages have added it. For example, the gmp package provides an interface to the gmp library which supports arbitrarily large integers and arbitary precision rationals. 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] Precision of values 53 bits
FAQ 7.31 On Thursday, January 10, 2013, Stephan Mueller smuel...@atsec.com wrote: Hi, I am working with large numbers and identified that R looses precision for such high numbers. The precision is lost exactly when the number is equal or larger than 53 bits. See the following output which shows that the numbers below 53 bit have proper precision: 2^53 [1] 9007199254740992 2^53-1 [1] 9007199254740991 2^53-2 [1] 9007199254740990 Now, see the numbers above 53 bit: 2^53 [1] 9007199254740992 2^53+1 [1] 9007199254740992 2^53+2 [1] 9007199254740994 2^53+3 [1] 9007199254740996 2^53+4 [1] 9007199254740996 Is there a solution to the problem? Thanks a lot Stephan __ 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. -- Jim Holtman Data Munger Guru What is the problem that you are trying to solve? Tell me what you want to do, not how you want to do it. [[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] multiple versions of function
You could make your 'f' a generic function and define methods
for various types. E.g., using S3 generics, define
f - function(a, b) UseMethod(f)
f.default - function(a, b) 10 * a + b
f.data.frame - function(df) f(df$a, df$b)
and use them as
f(b=5:7, a=1:3)
[1] 15 26 37
f(1:3, 5:7)
[1] 15 26 37
d - data.frame(b=5:7, a=1:3)
f(d)
[1] 15 26 37
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 ivo welch
Sent: Wednesday, January 09, 2013 1:00 PM
To: David Winsemius
Cc: r-help
Subject: Re: [R] multiple versions of function
mea culpa.
f - function(...) {
## parse out the arguments and then do something with them
}
## all of these should result in the same actions
f(2,3) ## interprets a to be first and b to be second
f(a=2,b=3)
f(b=3,a=2)
f(data.frame(a=2,b=3))
f(data.frame(b=3,a=1))
On Tue, Jan 8, 2013 at 8:00 AM, David Winsemius dwinsem...@comcast.net
wrote:
On Jan 7, 2013, at 6:58 PM, ivo welch wrote:
hi david---can you give just a little more of an example? the
function should work with call by order, call by name, and data frame
whose columns are the names. /iaw
It is I who should be expecting you to provide an example.
-- David.
__
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] Sweave, Texshop, and sync with included Rnw file
On 13-01-09 9:09 PM, Duncan Murdoch wrote:
On 13-01-09 3:25 PM, michele caseposta wrote:
Hello everyone.
I am in the process of writing a book in Latex with Texshop, on Mac.
This book contains a lot of R code, hence the need to use Sweave.
I was able to compile Rnw files, and to sync back and forth from the pdf to the
source Rnw.
My problem now is that the book is divided in Chapters, and every chapter is in
its own Rnw file.
I can compile them from the main one (book.Rnw) using the directive
\SweaveInput{chapter1.Rnw}
The problem stands in the fact that like this I am missing synchronization
between the pdf and the source Rnw. If part of text is in book.Rnw I can
synchronize, but if the text is in one of the included files, it just doesn't
work.
I am using the sweave engine found in the following webpage:
http://cameron.bracken.bz/synctex-with-sweavepgfsweave-in-texshoptexworks
Has anybody succeeded in synchronizing with included Rnw files?
This is a problem addressed by my patchDVI package, available on
R-forge. You have a main file (which can be .tex or .Rnw), and put code
at the start of each .Rnw file to indicate where to find it. Then you
just run Sweave on one of the chapters, and it automatically produces
the full document.
The sample document here:
http://www.umanitoba.ca/statistics/seminars/2011/3/4/duncan-murdoch-using-sweave-R/
includes an appendix describing how to set this up with TeXShop.
I just committed an update to the vignette in patchDVI giving a quick
version of the instructions for basic use. Version 1.8.1585 has the new
vignette.
I should get around to pushing it to CRAN one of these days...
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] Find the functional relationship between two variables in R?
On 10 January 2013 15:04, jt...@mappi.helsinki.fi wrote: Hi, I have two variables x and y and the functional relationship between x and y is like: y=x^2+log(x). My question is that is it possible to apply some method to reconstruct the functional form based on the training data that is generated from it? I understand that there are many methods for regresstion The answer is 42. __ 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] help with knit_hooks
This is the only public reference at the moment:
http://yihui.name/knitr/hooks (which has explained why your hook does
not work)
Or learn by examples: https://github.com/yihui/knitr-examples (e.g. example 045)
Or in the spirit of Luke, use the source!, see https://github.com/yihui/knitr
Regards,
Yihui
--
Yihui Xie xieyi...@gmail.com
Phone: 515-294-2465 Web: http://yihui.name
Department of Statistics, Iowa State University
2215 Snedecor Hall, Ames, IA
On Thu, Jan 10, 2013 at 6:51 AM, Francesco Sarracino
f.sarrac...@gmail.com wrote:
Dear R-listers,
does anybody can suggest some manual where I can learn more about how the
hooks in knitr work?
I am trying to enclose the output of an R command in the Latex verbatim
environment.
I defined a hook as follows:
knit_hooks$set(fsverb = function(x, options) {
paste(\\begin(verbatim)\n, x, \\end(verbatim)\n, sep = )
}
then I set a chunk as follows:
expl-emp, echo = TRUE, results = 'asis', fsverb = TRUE=
names(data)
@
but when I compile, I get a message of error saying that the compiler (in R
studio) is quitting from lines ...
I know that this is my fault because I don't know how to properly use the
hooks, therefore I am asking if anybody can advice a reference where I can
understand how the hooks and their options work and how to set them.
Thanks in advance for your kind help,
f.
--
Francesco Sarracino, Ph.D.
https://sites.google.com/site/fsarracino/
__
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] two phases sampling
Dear all, I have a question about estimation in two phases sampling. I' have a first sample of household from a complex sampling S1, a second sample is drawned from S2. from S2, I compute an estimator y2, for households of S1 not in S2 I set y2=0. I have an estimator y1 on S1 My indicator is y=y1+y2. So the variance of y is v(y)=v(y1)+v(y2)+cov(y1,y2). Sampling theory indicate how to compute v(y1) and v(y2) but how can I compute cov(y1,y2) ? Can the survey package help me for that ? Justin BEM BP 1917 Yaoundé Tél (237) 76043774 [[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] asCairoDevice issue
Hi All, I found this issue when using asCairoDevice to transforming splom scatter plot to my RGtk2 GUI: If I put the code in R GUI or using CairoPNG or Cairo_pdf() to draw the scatter plot, I can get it correctly: The codes are: (you can copy and paste to your R GUI) super.sym - trellis.par.get(superpose.symbol) plot.call-splom(~iris[1:4], groups = Species, data = iris, panel = panel.superpose, key = list(title = Three Varieties of Iris, columns = 3, points = list(pch = super.sym$pch[1:3], col = super.sym$col[1:3]), text = list(c(Setosa, Versicolor, Virginica plot.call However if I want to draw the same plot in the drawing area by asCairoDevice I lost all the colorful dots in the upper left and lower right, having only the diagonal charts: The codes are: (you can copy and paste to your R GUI) win- gtkWindowNew(show= FALSE) DA- gtkDrawingArea() asCairoDevice(DA) win$add(DA) win$show() plot.call Did I miss anything here? Thanks a lot!!! Regards, Yan [[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] sort matrix based on a specific order
Hi I have a character matrix with 2 columns A and B, If I want to sort the
matrix based on the column B, but based on a specific order of characters:
mat-cbind(c('w','x','y','z'),c('a','b','c','d'))
ind-c('c','b','d','a')
I want mat to be sorted by the sequence in ind:
[,1] [,2]
[1,] y c
[2,] x b
[3,] z d
[4,] w a
Is there any simple function that can do this?
Thanks
John
[[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] sort matrix based on a specific order
Define them as factors with a specified order for your sorting.
e.g.
x - factor(your_data, levels = c('c', 'b','d', 'a'))
On Thu, Jan 10, 2013 at 1:21 PM, array chip arrayprof...@yahoo.com wrote:
Hi I have a character matrix with 2 columns A and B, If I want to sort the
matrix based on the column B, but based on a specific order of characters:
mat-cbind(c('w','x','y','z'),c('a','b','c','d'))
ind-c('c','b','d','a')
I want mat to be sorted by the sequence in ind:
[,1] [,2]
[1,] y c
[2,] x b
[3,] z d
[4,] w a
Is there any simple function that can do this?
Thanks
John
[[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.
--
Jim Holtman
Data Munger Guru
What is the problem that you are trying to solve?
Tell me what you want to do, not how you want to do it.
__
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] how to generate a matrix by an my data.frame
Hello, Here are two ways. dat - read.table(text = id1id2 value 2353 2353 0.096313 2353 2409 0.301773 [...etc...] 2356 2356 0 2356 2611 0 2611 2611 0 , header = TRUE) mat1 - matrix(nrow = 53, ncol = 53) # initialize with NA's mat1[upper.tri(mat1, diag = TRUE)] - dat$value mat2 - matrix(0, nrow = 53, ncol = 53) # initialize with zeros mat2[upper.tri(mat2, diag = TRUE)] - dat$value Hope this helps, Rui Barradas Em 10-01-2013 15:21, Yao He escreveu: Dear All It is a little hard to give a good small example of my question,so I will show the full data on the bottom and the attachment.Maybe some one could tell me an appropriate way to show it.I'm sorry for the inconvenience. Q:How to generate a 53*53 diagonal matrix by my data Some problems confused me are that: 1.Since it is a diagonal matrix,I have tried to transform col1 and col2 to rowindex and colindex ,but I don't know how to generate matrix by its value's index 2. As you see, the number of 2353 corresponding to other ids in col2 is 53,however,the number of 2409 corresponding to other ids in col2 is 52 and 2500 corresponding to 51 values and so on,so it is hard to use matrix() to generate it id1id2 value 2353 23530.096313 2353 24090.301773 2353 25000.169518 2353 25980.11274 2353 26100.107414 2353 23000.034492 2353 25070.037521 2353 25300.064125 2353 23270.029259 2353 23890.036423 2353 24080.029259 2353 24630.036423 2353 24200.04409 2353 25630.055038 2353 24620.046478 2353 22920.036369 2353 24050.036369 2353 25430.053413 2353 25570.058151 2353 25830.081512 2353 23220.044373 2353 25350.04847 2353 25360.035538 2353 25810.035538 2353 25700.07711 2353 24760.047081 2353 25340.047081 2353 22800.088264 2353 23160.073608 2353 23390.067307 2353 23310.061172 2353 23430.060425 2353 23520.041153 2353 22930.040764 2353 23380.045128 2353 24490.040764 2353 22960.061333 2353 24530.046074 2353 24600.060387 2353 24740.060387 2353 26030.060387 2353 22820.048065 2353 23130.05584 2353 25380.050873 2353 25220.065727 2353 24890.041023 2353 25640.039696 2353 25940.056946 2353 22740.060875 2353 24510.037468 2353 23210 2353 23560 2353 26110 2409 24090.096313 2409 25000.169518 2409 25980.11274 2409 26100.107414 2409 23000.034492 2409 25070.037521 2409 25300.064125 2409 23270.029259 2409 23890.036423 2409 24080.029259 2409 24630.036423 2409 24200.04409 2409 25630.055038 2409 24620.046478 2409 22920.036369 2409 24050.036369 2409 25430.053413 2409 25570.058151 2409 25830.081512 2409 23220.044373 2409 25350.04847 2409 25360.035538 2409 25810.035538 2409 25700.07711 2409 24760.047081 2409 25340.047081 2409 22800.088264 2409 23160.073608 2409 23390.067307 2409 23310.061172 2409 23430.060425 2409 23520.041153 2409 22930.040764 2409 23380.045128 2409 24490.040764 2409 22960.061333 2409 24530.046074 2409 24600.060387 2409 24740.060387 2409 26030.060387 2409 22820.048065 2409 23130.05584 2409 25380.050873 2409 25220.065727 2409 24890.041023 2409 25640.039696 2409 25940.056946 2409 22740.060875 2409 24510.037468 2409 23210 2409 23560 2409 26110 2500 25000.048615 2500 25980.051979 2500 26100.041031 2500 23000.032974 2500 25070.052788 2500 25300.041435 2500 23270.038071 2500 23890.051659 2500 24080.038071 2500 24630.051659 2500 24200.052635 2500 25630.07872 2500 24620.048615 2500 22920.044365 2500 24050.044365 2500 25430.04277 2500 25570.051109 2500 25830.047409 2500 23220.054512 2500 25350.039368 2500 25360.041763 2500 25810.041763 2500 25700.063148 2500 24760.043755 2500 25340.043755 2500 22800.063164 2500 23160.083961 2500 23390.074127 2500 23310.051094 2500 23430.060066 2500 23520.038208 2500 22930.048698 2500 23380.048218 2500 24490.048698 2500 22960.073212 2500 24530.070595 2500 24600.073677 2500 24740.073677 2500 26030.073677 2500 22820.073677 2500 23130.068443 2500 25380.079865 2500 25220.059723 2500 24890.049873 2500 25640.087639 2500 25940.05175 2500 22740.043396 2500 24510.046532 2500 23210 2500 23560 2500 26110 2598 25980.040619 2598 26100.156075 2598 23000.049416
Re: [R] sort matrix based on a specific order
You can use factor() or match() to specify a particular order. E.g.,
mat-cbind(c('w','x','y','z'),c('a','b','c','d'))
ind-c('c','b','d','a')
mat[ order(match(mat[,2], ind)), ]
[,1] [,2]
[1,] y c
[2,] x b
[3,] z d
[4,] w a
mat[ order( factor(mat[,2], levels=ind) ), ]
[,1] [,2]
[1,] y c
[2,] x b
[3,] z d
[4,] w a
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 array chip
Sent: Thursday, January 10, 2013 10:22 AM
To: r-help@r-project.org
Subject: [R] sort matrix based on a specific order
Hi I have a character matrix with 2 columns A and B, If I want to sort the
matrix based on
the column B, but based on a specific order of characters:
mat-cbind(c('w','x','y','z'),c('a','b','c','d'))
ind-c('c','b','d','a')
I want mat to be sorted by the sequence in ind:
[,1] [,2]
[1,] y c
[2,] x b
[3,] z d
[4,] w a
Is there any simple function that can do this?
Thanks
John
[[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] sort matrix based on a specific order
more complete example
mat-cbind(c('w','x','y','z'),c('a','b','c','d'))
matOrd - mat[order(factor(mat[,2], levels = c('c', 'b', 'd','a'))), ]
matOrd
[,1] [,2]
[1,] y c
[2,] x b
[3,] z d
[4,] w a
On Thu, Jan 10, 2013 at 1:21 PM, array chip arrayprof...@yahoo.com wrote:
Hi I have a character matrix with 2 columns A and B, If I want to sort the
matrix based on the column B, but based on a specific order of characters:
mat-cbind(c('w','x','y','z'),c('a','b','c','d'))
ind-c('c','b','d','a')
I want mat to be sorted by the sequence in ind:
[,1] [,2]
[1,] y c
[2,] x b
[3,] z d
[4,] w a
Is there any simple function that can do this?
Thanks
John
[[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.
--
Jim Holtman
Data Munger Guru
What is the problem that you are trying to solve?
Tell me what you want to do, not how you want to do it.
__
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] sort matrix based on a specific order
Hello,
Try the following. order() gives you a permutation of the vector 'ind'
and to order that permutation gives its inverse.
mat - cbind(c('w','x','y','z'),c('a','b','c','d'))
ind - c('c','b','d','a')
ord - order(ind)
mat[order(ord), ]
Hope this helps,
Rui Barradas
Em 10-01-2013 18:21, array chip escreveu:
Hi I have a character matrix with 2 columns A and B, If I want to sort the
matrix based on the column B, but based on a specific order of characters:
mat-cbind(c('w','x','y','z'),c('a','b','c','d'))
ind-c('c','b','d','a')
I want mat to be sorted by the sequence in ind:
[,1] [,2]
[1,] y c
[2,] x b
[3,] z d
[4,] w a
Is there any simple function that can do this?
Thanks
John
[[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.
Re: [R] sort matrix based on a specific order
Thank you all for the suggestions, fantastic!
From: Rui Barradas ruipbarra...@sapo.pt
Cc: r-help@r-project.org r-help@r-project.org
Sent: Thursday, January 10, 2013 10:32 AM
Subject: Re: [R] sort matrix based on a specific order
Hello,
Try the following. order() gives you a permutation of the vector
'ind' and to order that permutation gives its inverse.
mat - cbind(c('w','x','y','z'),c('a','b','c','d'))
ind - c('c','b','d','a')
ord - order(ind)
mat[order(ord), ]
Hope this helps,
Rui Barradas
Em 10-01-2013 18:21, array chip escreveu:
Hi I have a character matrix with 2 columns A and B, If I want to sort the
matrix based on the column B, but based on a specific order of characters:
mat-cbind(c('w','x','y','z'),c('a','b','c','d'))
ind-c('c','b','d','a') I want mat to be sorted by the sequence in ind:
[,1] [,2]
[1,] y c
[2,] x b
[3,] z d
[4,] w a Is there any simple function that can do this? Thanks John
[[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.
Re: [R] Titles - main and subtitle won't plot with errbar
You should contact the maintainer of package Hmisc: Maintainer: Charles Dupont charles.dupont at vanderbilt.edu As you note, it would not be difficult to use the titles() function to get what you want or a plot command to set up but not plot data followed by the errbar() with add=TRUE. -- David L Carlson Associate Professor of Anthropology Texas AM University College Station, TX 77843-4352 -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-bounces@r- project.org] On Behalf Of masepot Sent: Thursday, January 10, 2013 8:54 AM To: r-help@r-project.org Subject: [R] Titles - main and subtitle won't plot with errbar Hi, I'm struggling with errbar graphics. I'm trying to plot an x-y graph with correct labelling, however can't seem to get main and sub to show on my graph. They do work when I use title(main=, etc, but this will make it look at lot messier,I'll have to blank out ylab= , and I need to try and get the titles to update automatically according to my excel column headings and paste function. Example code require(Hmisc) data1-array(sample(1:100,35),dim=c(5,7)) data1[,1]-1:5 sd=apply(data1[,2:7],1,sd) mean=rowMeans(data1[,2:7]) #the base plot works errbar(x=data1[,1],y=mean,yplus=mean+sd,yminus=mean-sd) #titles are shown correctly using plot plot(x=data1[,1],y=mean,ylab=value,xlab=time,main=Title,sub=subt itle) #the ylab and xlab update correctly, however main and sub don't? errbar(x=data1[,1],y=mean,yplus=mean+sd,yminus=mean- sd,ylab=value,xlab=time,main=Title,sub=subtitle) (my original is a lot more complicated as I'm reading from excel, and have managed to plot SD and mean for my dataset, but for some reason the main and sub commands aren't working in errbar!) Thanks, Laura -- View this message in context: http://r.789695.n4.nabble.com/Titles- main-and-subtitle-won-t-plot-with-errbar-tp4655149.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. __ 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] sort matrix based on a specific order
HI,
Try this:
mat[match(ind,mat[,2]),]
# [,1] [,2]
#[1,] y c
#[2,] x b
#[3,] z d
#[4,] w a
A.K.
- Original Message -
From: array chip arrayprof...@yahoo.com
To: r-help@r-project.org r-help@r-project.org
Cc:
Sent: Thursday, January 10, 2013 1:21 PM
Subject: [R] sort matrix based on a specific order
Hi I have a character matrix with 2 columns A and B, If I want to sort the
matrix based on the column B, but based on a specific order of characters:
mat-cbind(c('w','x','y','z'),c('a','b','c','d'))
ind-c('c','b','d','a')
I want mat to be sorted by the sequence in ind:
[,1] [,2]
[1,] y c
[2,] x b
[3,] z d
[4,] w a
Is there any simple function that can do this?
Thanks
John
[[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.
[R] Questions about the glht function for planned comparison
Hi all, I've posted this question before, but did not get any reply. I post it again here and see if anybody can help. Thank you. I have a nested model with the following effects fixed: treatments random: experiment_date I used lme() to model the data mod1 - lme(N_cells ~treatments-1, random=~1|experiment_date, method='ML') Then I want to compare all the other treatments to the control (included in the treatments in mod1). After a fair amount of searching around, I decided to use glht() from the multcomp package (any other suggestions?). lvl.treatments=table(treatments) K = contrMat(lvl.treatments,type='Dunnett',base=1) mc-glht(mod1, linfct=mcp(treatments=K),alternative='greater') But I got the following error: Error in contr.treatment(n = 0L) : not enough degrees of freedom to define contrasts I tried to extract the df parameter using modelparm(), but the function couldn't be applied to lme: Error in UseMethod(modelparm) : no applicable method for 'modelparm' applied to an object of class lme The degree of freedom of the fixed effect was 194. I tried to specify the number in glht, but got the same error as not enough degrees of freedom to define contrasts. Does anyone know what's happening and how I could possibly solve the problem? Thank you so much. Best, Wendy [[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] same model, different coefficients
Hello R-help subscribers, I am analyzing a data set using a mixed logit model, and I have recently discovered some curious behavior. I am hoping you all can help. I first ran the following model in December 2012. lmer(Response.binary ~ ItemType.c * Block + (1 | Subject) + (1 | Word), data=lexdec, family=binomial) I then took a break from the data for the holidays. I returned to the data yesterday and discovered that running the exact same model on the exact same data set yields different output. The overall patterns are the same, but the coefficients, variance estimates, and model fits (AIC, BIC) differ. The model outputs from the old and current attempt are appended below. I have triple checked the code and the data set to ensure that what I'm working with now is the same as in December. Having found no differences, I can only suspect that some function has changed. During my hiatus from these data, I updated plyr and its dependencies (and maybe some other packages). But to my understanding, these updates mostly concerned documentation, not algorithms. Any ideas then about why the model outputs differ? Thank you for your help! Kodi * OLD MODEL OUTPUT*: Generalized linear mixed model fit by the Laplace approximation Formula: Response.binary ~ ItemType.c * Block + (1 + ItemType.c + Block | Subject) + (1 | Word) Data: lexdec AIC BIC logLik deviance 4788 4957 -2370 4740 Random effects: Groups Name Variance Std.Dev. Corr Word(Intercept) 1.66447 1.29014 Subject (Intercept) 0.50865 0.71320 ItemType.cFV-L 0.89270 0.94483 0.261 ItemType.cFV-R 1.26385 1.12421 0.210 0.978 ItemType.cFV-B 1.33556 1.15566 0.143 0.916 0.979 Blockpost0.93878 0.96891 -0.349 -0.093 0.037 0.163 Number of obs: 8500, groups: Word, 298; Subject, 17 Fixed effects: Estimate Std. Error z value Pr(|z|) (Intercept) 3.8233 0.2688 14.225 2e-16 *** ItemType.cFV-L -6.0547 0.3749 -16.149 2e-16 *** ItemType.cFV-R -6.8649 0.4130 -16.621 2e-16 *** ItemType.cFV-B -7.3542 0.4285 -17.164 2e-16 *** Blockpost 0.9754 0.3238 3.013 0.00259 ** ItemType.cFV-L:Blockpost0.5921 0.2725 2.173 0.02980 * ItemType.cFV-R:Blockpost0.5835 0.2926 1.994 0.04612 * ItemType.cFV-B:Blockpost -0.2718 0.3083 -0.882 0.37793 --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 *NEW MODEL OUTPUT: *Generalized linear mixed model fit by the Laplace approximation Formula: Response.binary ~ ItemType.c * Block + (1 + ItemType.c + Block | Subject) + (1 | Word) Data: lexdec AIC BIC logLik deviance 4791 4961 -2372 4743 Random effects: Groups NameVariance Std.Dev. Corr Word(Intercept) 1.57837 1.25633 Subject (Intercept) 0.58476 0.76470 ItemType.cFV-L 0.92922 0.96396 -0.105 ItemType.cFV-R 1.36398 1.16790 -0.241 0.990 ItemType.cFV-B 1.59667 1.26360 -0.323 0.956 0.978 Blockpost 1.03413 1.01692 -0.511 0.198 0.264 0.406 Number of obs: 8500, groups: Word, 298; Subject, 17 Fixed effects: Estimate Std. Error z value Pr(|z|) (Intercept) 3.3659 0.2704 12.448 2e-16 *** ItemType.cFV-L -5.6366 0.3710 -15.193 2e-16 *** ItemType.cFV-R -6.3466 0.4128 -15.376 2e-16 *** ItemType.cFV-B -6.6767 0.4372 -15.271 2e-16 *** Blockpost 1.1657 0.3288 3.546 0.000391 *** ItemType.cFV-L:Blockpost0.5119 0.2681 1.909 0.056243 . ItemType.cFV-R:Blockpost0.4834 0.2859 1.691 0.090886 . ItemType.cFV-B:Blockpost -0.4394 0.3002 -1.463 0.143336 --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 * * -- Kodi Weatherholtz Ph.D. Student Department of Linguistics The Ohio State University [[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] Subset in, not in
Hello, I need to subset my dataframe into 2 parts; in: mm - subset(agr1, subset=lmpcrd %in% c(11697,149823,7654)) not in: but where do I stick the ! in the above? I've tried every position. Thanks for your help. -- View this message in context: http://r.789695.n4.nabble.com/Subset-in-not-in-tp4655178.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] Problem getting loess tricubic weights
To further the understanding of the loess fit and how the tricube weight work you may want to look at the loess.demo function in the TeachingDemos package. It will create a scatterplot of the data and show the loess fit, then when you click on the plot it will show the weights used for predicting at that point and the fitted line/curve for that point. Click at another place in the graph and it will show the weights and local fit corresponding to that point. On Tue, Jan 8, 2013 at 9:24 PM, Joyce Lin joyceli...@gmail.com wrote: Thank you Mr Gunter! I will look into it. On Wed, Jan 9, 2013 at 11:59 AM, Bert Gunter gunter.ber...@gene.com wrote: As this does not seem to have been answered... I believe you may misunderstand how loess works. The tricube weights are part of the smoothing algorithm and change with each local fit, not fixed weights for observations, which is what the weights argument provides (and initially multiplies the tricube weight, IIRC). I suggest you consult ?predict.loess to get standard deviations of fitted values at existing or new points. -- Bert On Tue, Jan 8, 2013 at 12:57 AM, Joyce Lin joyceli...@gmail.com wrote: Hi I am trying to get the tricube weights from the loess outputs as I need to calculate an error function which requires the weight. So I have used the following example from the R: cars.lo - loess(dist ~ speed, cars, span=0.5, degree=1, family=symmetric) Then i try to get the weights: cars.lo$weights [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 The results are all 1 so i dont think that the tricube weighting are set. May I know what other parameters do i need to tweak to set the weights to tricube weights? Thank you. -- Best regards Joyce Lin [[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. -- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm -- Best regards Joyce Lin [[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. -- Gregory (Greg) L. Snow Ph.D. 538...@gmail.com [[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] Subset in, not in
Did you try: mm - subset(agr1, subset= !(lmpcrd %in% c(11697,149823,7654))) (Actually the parentheses that I added are not necessary, but they make me feel better.) Pat On 10/01/2013 19:54, ramoss wrote: Hello, I need to subset my dataframe into 2 parts; in: mm - subset(agr1, subset=lmpcrd %in% c(11697,149823,7654)) not in: but where do I stick the ! in the above? I've tried every position. Thanks for your help. -- View this message in context: http://r.789695.n4.nabble.com/Subset-in-not-in-tp4655178.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. -- Patrick Burns pbu...@pburns.seanet.com twitter: @portfolioprobe http://www.portfolioprobe.com/blog http://www.burns-stat.com (home of 'Some hints for the R beginner' and 'The R Inferno') __ 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] problem adding curve/abline
I believe the problem could be that xyplot uses grid graphics and plot.new and curve are base graphics functions and the 2 graphics systems (grid and base) don't play nicely together without a little extra work. In general the gridBase package helps them play nicely, but I am not sure that it will be the best approach in your case. I would look at rewriting the effects of curve in a panel function for xyplot using the llines function (notice the 2 'l's at the front) or using only base graphics. On Wed, Jan 9, 2013 at 9:57 AM, Elisabeth Van Beveren eli...@hotmail.comwrote: Hey, I'm stuck on something I already did before (just a different kind of database), and whatever I try, it doesn't work anymore. So thanks for your help. Here's how my data approximately looks like: year season replicate sizefreq weight 2000 summer ch1 6 1 45 2000 summer ch1 6.5 12 46 2000 summer ch1 7 33 470 I have 2 years (2000 and 2001) and 2 seizons (winter and summer). I wanted to plot weight~size, with 2 groups (year and seizon), so here's my shortened script for that: database$groups=paste(database$seizon,database2$year,sep= ) xyplot(database$weight~database2$size, groups=database$groups, par.settings=list(superpose.symbol=list(col=col.list,pch=c(21,16,21,16))), auto.key=list(corner=c(0.1,0.9),lines=F,points=T)) Which works fine, the problem comes when I try to add 2 exponential curves to the data (the 2 seizons). I tried this: summ=subset(database,seizon==summer) modsumm=nls(summ$weight~exp(a+b*summ$size), data=summ, start=list(a=0,b=0)) exposumm=curve(exp(0.05354+0.19872*x), from=0, to=22, add=T, lwd=1, col=blue,lty=1) After having to add plot.new() in the front, the line does or not show up, or shows up but wrongly placed. I thought this might be because of the subset, so I wanted to do something like this: modsumm=nls(weight~exp(a+b* size), data=engsAGG2[seizon==summer], start=list(a=0,b=0)) which returns: undefined columns selected Thanks in advance for the reply. [[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. -- Gregory (Greg) L. Snow Ph.D. 538...@gmail.com [[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] Subset in, not in
On Jan 10, 2013, at 12:04 PM, Patrick Burns wrote: Did you try: mm - subset(agr1, subset= !(lmpcrd %in% c(11697,149823,7654))) And it might be noted that this is part of the last example on the help page: ?'%in%' -- David. (Actually the parentheses that I added are not necessary, but they make me feel better.) Pat On 10/01/2013 19:54, ramoss wrote: Hello, I need to subset my dataframe into 2 parts; in: mm - subset(agr1, subset=lmpcrd %in% c(11697,149823,7654)) not in: but where do I stick the ! in the above? I've tried every position. Thanks for your help. -- David Winsemius, MD 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] merging command
Dear Arun,thankyou very much...
Date: Thu, 10 Jan 2013 12:02:31 -0800
From: smartpink...@yahoo.com
Subject: Re: merging command
To: eliza_bo...@hotmail.com
CC: r-help@r-project.org
HI Eliza,
You could do this:
set.seed(15)
mat1-matrix(sample(1:800,124*12,replace=TRUE),nrow=12) # smaller dataset
#Your codes
list1-list()
for(i in 1:ncol(mat1)){
list1[[i]]-t(apply(mat1,1,function(x) x[i]-x))
list1}
x-list1
x-matrix(unlist(x),nrow=12)
x-abs(x)
y-colSums(x, na.rm=FALSE)
z-matrix(y,ncol=10)
z-as.dist(z)
z
#1 2 3 4 5 6 7 8 9
#2 319
#3 459 516
#4 385 504 260
#5 365 282 506 520
#6 318 363 373 305 383
#7 382 277 459 457 363 370
#8 526 521 431 443 523 472 608
#9 329 534 358 374 382 393 467 429
#10 364 377 393 365 419 420 346 472 489
#Modified code
z1-as.dist(do.call(cbind,lapply(seq_len(ncol(mat1)),function(i)
colSums(abs(t(apply(mat1,1, function(x) x[i]-x))),na.rm=FALSE
z1
# 1 2 3 4 5 6 7 8 9
#2 319
#3 459 516
#4 385 504 260
#5 365 282 506 520
#6 318 363 373 305 383
#7 382 277 459 457 363 370
#8 526 521 431 443 523 472 608
#9 329 534 358 374 382 393 467 429
#10 364 377 393 365 419 420 346 472 489
A.K.
From: eliza botto eliza_bo...@hotmail.com
To: smartpink...@yahoo.com smartpink...@yahoo.com
Sent: Thursday, January 10, 2013 9:13 AM
Subject: merging command
Dear Arun,
i need you expertise to merge the following commands in to one step command.
mat1-m
list1-list()
for(i in 1:ncol(mat1)){
list1[[i]]-t(apply(mat1,1,function(x) x[i]-x))
list1}
x-list1
x-matrix(unlist(x),nrow=12)
x-abs(x)
y-colSums(x,
na.rm=FALSE)
z-matrix(y,
ncol=124)
z-as.dist(z)
i needed that distance to be executed in one command by merging all these
commands.
is it possible??
thanks in advance
elisa
[[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] Titles - main and subtitle won't plot with errbar
On Jan 10, 2013, at 6:54 AM, masepot wrote: Hi, I'm struggling with errbar graphics. I'm trying to plot an x-y graph with correct labelling, however can't seem to get main and sub to show on my graph. If you are like me you are perhaps being surprised by the fact that the subtitle s not appearing just below the title but rather below the xlab. If you are also like me it will take you some time to learn that the 'xlab' and the 'ylab' have nothing to do with labels. -- David. They do work when I use title(main=, etc, but this will make it look at lot messier,I'll have to blank out ylab= , and I need to try and get the titles to update automatically according to my excel column headings and paste function. Example code require(Hmisc) data1-array(sample(1:100,35),dim=c(5,7)) data1[,1]-1:5 sd=apply(data1[,2:7],1,sd) mean=rowMeans(data1[,2:7]) #the base plot works errbar(x=data1[,1],y=mean,yplus=mean+sd,yminus=mean-sd) #titles are shown correctly using plot plot(x=data1[,1],y=mean,ylab=value,xlab=time,main=Title,sub=subtitle) #the ylab and xlab update correctly, however main and sub don't? errbar(x=data1[,1],y=mean,yplus=mean+sd,yminus=mean-sd,ylab=value,xlab=time,main=Title,sub=subtitle) (my original is a lot more complicated as I'm reading from excel, and have managed to plot SD and mean for my dataset, but for some reason the main and sub commands aren't working in errbar!) Thanks, Laura -- View this message in context: http://r.789695.n4.nabble.com/Titles-main-and-subtitle-won-t-plot-with-errbar-tp4655149.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. 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.
[R] transparency in segments()
Dear all, i would like to plot each value from my datasets as segment with defined transparency However, I didnt find out how to set the transparency. definition by col= in par() or segments() doesnt seem to work any suggestions? Thanks in advance. Kind regards, Robert Pazur example code: xx2 -read.table(http://www.scandinavia.sk/data/R/0_05.csv;, sep=;, header=T) plot(xx2$MEAN_PERI,ylim=c(50,350),xlim=c(1,3),log= x, yaxt='n',ylab=,xlab=mean, type = n) segments(xx2$MEAN_PERI, 60,xx2$MEAN_PERI, 100, tcl=-.2) xx3 -read.table(http://www.scandinavia.sk/data/R/0_10.csv;, sep=;, header=T) segments(xx3$MEAN_PERI, 120,xx3$MEAN_PERI, 160, tcl=-.2) xx4 -read.table(http://www.scandinavia.sk/data/R/0_15.csv;, sep=;, header=T) segments(xx4$MEAN_PERI, 180,xx4$MEAN_PERI, 220, tcl=-.2) xx5 -read.table(http://www.scandinavia.sk/data/R/0_20.csv;, sep=;, header=T) segments(xx5$MEAN_PERI, 240,xx5$MEAN_PERI, 280, tcl=-.2) xx6 -read.table(http://www.scandinavia.sk/data/R/0_25.csv;, sep=;, header=T) segments(xx6$MEAN_PERI, 300,xx6$MEAN_PERI, 340, tcl=-.2) [[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] sort matrix based on a specific order
mat[match(ind, mat[, 2]), ]
[,1] [,2]
[1,] y c
[2,] x b
[3,] z d
[4,] w a
though you need to take care if you have cases where ind will contains letters
that are not in mat[, 2] and so on (check ?match).
Best,
I
On 10 Jan 2013, at 18:21, array chip arrayprof...@yahoo.com wrote:
Hi I have a character matrix with 2 columns A and B, If I want to sort the
matrix based on the column B, but based on a specific order of characters:
mat-cbind(c('w','x','y','z'),c('a','b','c','d'))
ind-c('c','b','d','a')
I want mat to be sorted by the sequence in ind:
[,1] [,2]
[1,] y c
[2,] x b
[3,] z d
[4,] w a
Is there any simple function that can do this?
Thanks
John
[[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.
--
Dr Ioannis Kosmidis
Department of Statistical Science,
University College,
London, WC1E 6BT, UK
Webpage: http://www.ucl.ac.uk/~ucakiko
__
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] Sweave, Texshop, and sync with included Rnw file
Hi everybody,
thanks for the replies.
I might have not explained the problem completely.
Duncan Mackay:
Yes, I am already having a master file and separate Rnw files.
Duncan Murdock:
I am using patchDVI in the TexShop Sweave engine.
Sync works flawlessly between the master file and the pdf produced by pdflatex.
My problem is that I don't seem to be able to obtain sync between the
*included* Rnws and the pdf, either way.
The sweave engine is as follows:
#!/bin/bash
R CMD Sweave $1
latexmk -pdf -silent -pdflatex=âpdflatex âshell-escape
âsynctex=1â²${1%.*}
Rscript -e library(âpatchDVIâ);patchSynctex(â${1%.*}.synctex.gzâ)
Funny thing is that the sync works in texworks, using the following Rscript line
patchDVI::SweavePDF('$fullname',stylepath=FALSE)
I tried to mix and match configurations between texshop and texworks but I had
no luck
On Jan 10, 2013, at 11:23 AM, Duncan Murdoch wrote:
On 13-01-09 9:09 PM, Duncan Murdoch wrote:
On 13-01-09 3:25 PM, michele caseposta wrote:
Hello everyone.
I am in the process of writing a book in Latex with Texshop, on Mac.
This book contains a lot of R code, hence the need to use Sweave.
I was able to compile Rnw files, and to sync back and forth from the pdf to
the source Rnw.
My problem now is that the book is divided in Chapters, and every chapter
is in its own Rnw file.
I can compile them from the main one (book.Rnw) using the directive
\SweaveInput{chapter1.Rnw}
The problem stands in the fact that like this I am missing synchronization
between the pdf and the source Rnw. If part of text is in book.Rnw I can
synchronize, but if the text is in one of the included files, it just
doesn't work.
I am using the sweave engine found in the following webpage:
http://cameron.bracken.bz/synctex-with-sweavepgfsweave-in-texshoptexworks
Has anybody succeeded in synchronizing with included Rnw files?
This is a problem addressed by my patchDVI package, available on
R-forge. You have a main file (which can be .tex or .Rnw), and put code
at the start of each .Rnw file to indicate where to find it. Then you
just run Sweave on one of the chapters, and it automatically produces
the full document.
The sample document here:
http://www.umanitoba.ca/statistics/seminars/2011/3/4/duncan-murdoch-using-sweave-R/
includes an appendix describing how to set this up with TeXShop.
I just committed an update to the vignette in patchDVI giving a quick version
of the instructions for basic use. Version 1.8.1585 has the new vignette.
I should get around to pushing it to CRAN one of these days...
Duncan Murdoch
[[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] RPART: Including the expense of predictor variables
I know that rpart has a complexity parameter that adjusts the number of nodes in a model. I also know that a loss function allows one to weight misclassifications of different types. However, some of my predictor variables are much more expensive dollar-wise to use than others. Is there a way to weight the predictor variables such that rpart will not use an expensive variable (or will only send limited fractions of the population to the node) if there is not a comparatively large decrease in misclassifications after splitting by that variable? Thanks, Lee __ 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] Lambert W question
RSiteSearch(lambert) I don't know anything about OpenBUGS, but implementations of the Lambert W function exist, or you could roll your own. --- Jeff NewmillerThe . . Go Live... DCN:jdnew...@dcn.davis.ca.usBasics: ##.#. ##.#. Live Go... Live: OO#.. Dead: OO#.. Playing Research Engineer (Solar/BatteriesO.O#. #.O#. with /Software/Embedded Controllers) .OO#. .OO#. rocks...1k --- Sent from my phone. Please excuse my brevity. Andras Farkas motyoc...@yahoo.com wrote: Dear All, � I am using the following model equation: � k*(lambertW_base(b=0,((a)/k)*exp(((a)-z*(t-t0))/k))) � I would like to run this through OpenBUGS, but it does not recognize the lambert function. Would you have any thoughts on how to re-vrite this equation matemathically so that it could be run on OpenBUGS? � apreciate the help, � Sincerely, � Andrs [[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] merging command
HI Eliza,
You could do this:
set.seed(15)
mat1-matrix(sample(1:800,124*12,replace=TRUE),nrow=12) # smaller dataset
#Your codes
list1-list()
for(i in 1:ncol(mat1)){
list1[[i]]-t(apply(mat1,1,function(x) x[i]-x))
list1}
x-list1
x-matrix(unlist(x),nrow=12)
x-abs(x)
y-colSums(x, na.rm=FALSE)
z-matrix(y,ncol=10)
z-as.dist(z)
z
# 1 2 3 4 5 6 7 8 9
#2 319
#3 459 516
#4 385 504 260
#5 365 282 506 520
#6 318 363 373 305 383
#7 382 277 459 457 363 370
#8 526 521 431 443 523 472 608
#9 329 534 358 374 382 393 467 429
#10 364 377 393 365 419 420 346 472 489
#Modified code
z1-as.dist(do.call(cbind,lapply(seq_len(ncol(mat1)),function(i)
colSums(abs(t(apply(mat1,1, function(x) x[i]-x))),na.rm=FALSE
z1
# 1 2 3 4 5 6 7 8 9
#2 319
#3 459 516
#4 385 504 260
#5 365 282 506 520
#6 318 363 373 305 383
#7 382 277 459 457 363 370
#8 526 521 431 443 523 472 608
#9 329 534 358 374 382 393 467 429
#10 364 377 393 365 419 420 346 472 489
A.K.
From: eliza botto eliza_bo...@hotmail.com
To: smartpink...@yahoo.com smartpink...@yahoo.com
Sent: Thursday, January 10, 2013 9:13 AM
Subject: merging command
Dear Arun,
i need you expertise to merge the following commands in to one step command.
mat1-m
list1-list()
for(i in 1:ncol(mat1)){
list1[[i]]-t(apply(mat1,1,function(x) x[i]-x))
list1}
x-list1
x-matrix(unlist(x),nrow=12)
x-abs(x)
y-colSums(x,
na.rm=FALSE)
z-matrix(y,
ncol=124)
z-as.dist(z)
i needed that distance to be executed in one command by merging all these
commands.
is it possible??
thanks in advance
elisa
__
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] Wald test for comparing coefficients across groups
Dear R users, my question concerns my interest in comparing the beta coefficients between two identical regression models in two (actually 3) groups. Disclaimer: I am quite new to R, so I might be missing some terminology that I have not come across. I am trying to find out if I can easily implement a Wald test in R for this, but the only relevant thing that I came across is this link (http://hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/aod/html/wald.test.html), which seems to be referring to comparing the coefficients of separate IVs within the same regression. What I am interested is to compare the coefficients of the same IVs across three regressions and across three groups. Some basic information about my model/design: My model includes 2 fixed factors (task, group) and 5 more IVs of interest. Choice ~ Task * Group + IV1 + IV2 + IV3 + IV4 + IV5 I established that the Task * Group interaction is significant (by running two models of a regression, one comparing a model with and one without the interaction term) and I am interested to explore group differences within each task in a separate logistic regression. Task has 3 levels, therefore I do 3 of the following regressions: Choice ~ Group + IV1 + IV2 + IV3 + IV4 + IV5, one for each level of Task. At this point I want to compare betas for the 5 IVs across my 3 groups (actually I don't even need to compare all 5, based on the theory behind it, if that is a concern for some in terms of groupwise errors). PS: I did a z-test across my critical IVs between groups, which is quite similar to Wald. I just want to do the Wald as well to doublecheck my results. PS2: I also did another thing that some people might suggest, which is to recode Group in 2 dummy variables and include those and the interaction terms within each of my 3 regressions. I could go one step back and do it in the original regression (Choice ~ Task * Group + IV1 + IV2 + IV3) but I want to avoid having to include the complicated 3-way interactions. Also, I don't think it is much of a problem to analyze different tasks individually, especially after I have already established that the Task * Group interaction is significant. So, can someone provide the code or point me to a place that a similar issue might already have been discussed, for the case of doing a Wald test to compare the beta coefficients of specific IVs of interest between different regressions/groups (same model used for all thesegroups)? Thanks in advance, Kostas [[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] piece-wise linear regression nls function
Instead of reinventing the wheel, why not use the segmented package? Having stored your data in a data frame X I did: require(segmented) m1 - lm(FM ~ BMIJS,data=X) m2 - segmented(m1,seg.Z=~BMIJS,psi=list(BMIJS=35)) which worked instantaneously. The break point is estimated as 41.580, with a standard error of 2.094 I then did: with(X, plot(BMIJS,FM)) plot(m2,add=TRUE) which seems to look as good as one can expect. I must say however that the plot of your data does not look to me as though a broken-stick model is appropriate. Why not just a straight line? cheers, Rolf Turner On 01/10/2013 02:33 PM, John Sorkin wrote: windows 7, R 2.12 I am trying to run a piecewise linear regression with a single knot, i.e. a regression composed of two straight lines where the two lines intersect at an x value given by the variable knot. I wish to estimate the slope of both lines, the value of knot, the x value where the two lines intersect, and an intercept. I am using the nls code below, and get the following error message: Error in nls(FM ~ blow * BMIJS + bhi * sapply(BMIJS - knot, max, 0), start = list(knot = 25, : singular gradient nls code: test - nls(FM~blow*BMIJS+bhi*sapply(BMIJS-knot,max,0),start=list(knot=25,blow=1,bhi=1),data=FeMaleData) summary(test) greatly shortened version of my data (the full data set has 450 records) FMBMIJS 2 55.878 40.57273 4 34.270 27.76939 5 20.123 21.73818 6 19.320 19.71203 9 49.701 43.55356 10 51.188 37.84742 11 46.753 37.71003 13 65.079 37.23438 14 37.097 36.81806 15 30.625 29.92783 17 50.617 42.42754 18 63.954 48.78709 20 29.790 26.97648 21 36.558 34.79373 22 41.275 33.03063 24 27.682 27.24508 26 37.968 35.41399 28 24.878 27.20250 30 47.513 35.77961 31 51.315 37.46032 33 41.944 36.40212 34 38.150 32.83818 35 60.719 42.48594 36 42.643 34.29355 38 40.728 32.42817 42 34.814 30.57573 43 32.896 29.32912 44 30.430 25.44183 46 48.986 37.90910 49 47.485 36.34642 52 46.312 38.64647 54 45.228 33.08783 55 45.391 35.86965 59 37.256 32.66507 60 27.367 28.49880 63 38.663 34.34131 64 34.527 29.57858 67 58.368 38.97266 68 13.473 17.35397 69 22.456 20.80958 71 28.829 25.50056 73 15.487 20.22202 76 18.313 21.38991 77 41.535 36.85707 78 56.124 40.51978 80 52.587 40.77256 81 24.991 25.48543 83 56.327 39.97214 84 70.836 36.52915 85 62.294 42.45244 86 39.689 35.18527 87 35.006 35.15136 88 47.378 37.54779 89 18.149 23.99236 90 33.041 28.10476 91 28.884 26.74443 92 37.670 32.25230 94 55.410 43.72364 99 34.461 35.05930 101 59.727 42.83035 102 41.913 35.64677 104 66.644 41.01642 105 55.250 43.86426 107 45.196 31.78370 108 36.476 33.45537 109 34.386 29.08402 110 39.277 36.98500 111 53.789 45.54654 112 33.077 29.09559 116 57.246 39.98031 120 52.546 40.12191 122 34.409 29.70977 123 31.188 28.75295 126 54.567 38.15226 129 19.193 22.71878 133 39.322 33.45712 134 41.415 31.28980 136 57.616 36.94016 140 28.162 24.40219 142 37.524 29.92673 143 29.611 29.15452 144 26.780 26.53462 146 47.219 35.14919 147 35.341 28.68955 148 44.827 37.68317 149 54.180 41.12226 150 41.636 30.00930 151 33.626 28.00164 156 34.334 29.64970 160 36.317 30.12031 161 46.823 35.64603 163 39.506 34.27740 164 61.619 39.20019 169 48.984 35.77558 171 66.467 41.59008 172 70.144 42.79996 173 37.324 31.56521 174 66.882 46.04938 182 54.239 38.21065 184 48.800 32.01630 __ 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] transparency in segments()
Hello, You can use ?rgb to set the transparency level. As an example, with alpha = 0.5 clr - c(rgb(1, 0, 0, 0.5), rgb(0, 0, 1, 0.5)) plot(0:1, 0:1, col = clr[1], lwd = 10, type = l) lines(0:1, 1:0, col = clr[2], lwd = 10) Hope this helps, Rui Barradas Em 10-01-2013 21:29, Robert Pazur escreveu: Dear all, i would like to plot each value from my datasets as segment with defined transparency However, I didnt find out how to set the transparency. definition by col= in par() or segments() doesnt seem to work any suggestions? Thanks in advance. Kind regards, Robert Pazur example code: xx2 -read.table(http://www.scandinavia.sk/data/R/0_05.csv;, sep=;, header=T) plot(xx2$MEAN_PERI,ylim=c(50,350),xlim=c(1,3),log= x, yaxt='n',ylab=,xlab=mean, type = n) segments(xx2$MEAN_PERI, 60,xx2$MEAN_PERI, 100, tcl=-.2) xx3 -read.table(http://www.scandinavia.sk/data/R/0_10.csv;, sep=;, header=T) segments(xx3$MEAN_PERI, 120,xx3$MEAN_PERI, 160, tcl=-.2) xx4 -read.table(http://www.scandinavia.sk/data/R/0_15.csv;, sep=;, header=T) segments(xx4$MEAN_PERI, 180,xx4$MEAN_PERI, 220, tcl=-.2) xx5 -read.table(http://www.scandinavia.sk/data/R/0_20.csv;, sep=;, header=T) segments(xx5$MEAN_PERI, 240,xx5$MEAN_PERI, 280, tcl=-.2) xx6 -read.table(http://www.scandinavia.sk/data/R/0_25.csv;, sep=;, header=T) segments(xx6$MEAN_PERI, 300,xx6$MEAN_PERI, 340, tcl=-.2) [[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] problem adding curve/abline
On Jan 10, 2013, at 12:04 PM, Greg Snow wrote: I believe the problem could be that xyplot uses grid graphics and plot.new and curve are base graphics functions and the 2 graphics systems (grid and base) don't play nicely together without a little extra work. In general the gridBase package helps them play nicely, but I am not sure that it will be the best approach in your case. I would look at rewriting the effects of curve in a panel function for xyplot using the llines function (notice the 2 'l's at the front) or using only base graphics. There is a lattice function `panel.curve` that might make this all hang together rather than hanging separately. -- David. On Wed, Jan 9, 2013 at 9:57 AM, Elisabeth Van Beveren eli...@hotmail.comwrote: Hey, I'm stuck on something I already did before (just a different kind of database), and whatever I try, it doesn't work anymore. So thanks for your help. Here's how my data approximately looks like: year season replicate sizefreq weight 2000 summer ch1 6 1 45 2000 summer ch1 6.5 12 46 2000 summer ch1 7 33 470 I have 2 years (2000 and 2001) and 2 seizons (winter and summer). I wanted to plot weight~size, with 2 groups (year and seizon), so here's my shortened script for that: database$groups=paste(database$seizon,database2$year,sep= ) xyplot(database$weight~database2$size, groups=database$groups, par.settings=list(superpose.symbol=list(col=col.list,pch=c(21,16,21,16))), auto.key=list(corner=c(0.1,0.9),lines=F,points=T)) Which works fine, the problem comes when I try to add 2 exponential curves to the data (the 2 seizons). I tried this: summ=subset(database,seizon==summer) modsumm=nls(summ$weight~exp(a+b*summ$size), data=summ, start=list(a=0,b=0)) exposumm=curve(exp(0.05354+0.19872*x), from=0, to=22, add=T, lwd=1, col=blue,lty=1) After having to add plot.new() in the front, the line does or not show up, or shows up but wrongly placed. I thought this might be because of the subset, so I wanted to do something like this: modsumm=nls(weight~exp(a+b* size), data=engsAGG2[seizon==summer], start=list(a=0,b=0)) which returns: undefined columns selected Thanks in advance for the reply. [[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. -- Gregory (Greg) L. Snow Ph.D. 538...@gmail.com [[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. 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] transparency in segments()
On Jan 10, 2013, at 1:29 PM, Robert Pazur wrote: Dear all, i would like to plot each value from my datasets as segment with defined transparency However, I didnt find out how to set the transparency. definition by col= in par() or segments() doesnt seem to work any suggestions? Try this: ?colors example code: xx2 -read.table(http://www.scandinavia.sk/data/R/0_05.csv;, sep=;, header=T) plot(xx2$MEAN_PERI,ylim=c(50,350),xlim=c(1,3),log= x, yaxt='n',ylab=,xlab=mean, type = n) segments(xx2$MEAN_PERI, 60,xx2$MEAN_PERI, 100, tcl=-.2) xx3 -read.table(http://www.scandinavia.sk/data/R/0_10.csv;, sep=;, header=T) segments(xx3$MEAN_PERI, 120,xx3$MEAN_PERI, 160, tcl=-.2) xx4 -read.table(http://www.scandinavia.sk/data/R/0_15.csv;, sep=;, header=T) segments(xx4$MEAN_PERI, 180,xx4$MEAN_PERI, 220, tcl=-.2) xx5 -read.table(http://www.scandinavia.sk/data/R/0_20.csv;, sep=;, header=T) segments(xx5$MEAN_PERI, 240,xx5$MEAN_PERI, 280, tcl=-.2) xx6 -read.table(http://www.scandinavia.sk/data/R/0_25.csv;, sep=;, header=T) segments(xx6$MEAN_PERI, 300,xx6$MEAN_PERI, 340, tcl=-.2) [[alternative HTML version deleted]] _ 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] Problem getting loess tricubic weights
Thank you Mr Snow. I will look into it. Best regards Joyce Lin On 11 Jan, 2013, at 3:55 AM, Greg Snow 538...@gmail.com wrote: To further the understanding of the loess fit and how the tricube weight work you may want to look at the loess.demo function in the TeachingDemos package. It will create a scatterplot of the data and show the loess fit, then when you click on the plot it will show the weights used for predicting at that point and the fitted line/curve for that point. Click at another place in the graph and it will show the weights and local fit corresponding to that point. On Tue, Jan 8, 2013 at 9:24 PM, Joyce Lin joyceli...@gmail.com wrote: Thank you Mr Gunter! I will look into it. On Wed, Jan 9, 2013 at 11:59 AM, Bert Gunter gunter.ber...@gene.com wrote: As this does not seem to have been answered... I believe you may misunderstand how loess works. The tricube weights are part of the smoothing algorithm and change with each local fit, not fixed weights for observations, which is what the weights argument provides (and initially multiplies the tricube weight, IIRC). I suggest you consult ?predict.loess to get standard deviations of fitted values at existing or new points. -- Bert On Tue, Jan 8, 2013 at 12:57 AM, Joyce Lin joyceli...@gmail.com wrote: Hi I am trying to get the tricube weights from the loess outputs as I need to calculate an error function which requires the weight. So I have used the following example from the R: cars.lo - loess(dist ~ speed, cars, span=0.5, degree=1, family=symmetric) Then i try to get the weights: cars.lo$weights [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 The results are all 1 so i dont think that the tricube weighting are set. May I know what other parameters do i need to tweak to set the weights to tricube weights? Thank you. -- Best regards Joyce Lin [[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. -- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm -- Best regards Joyce Lin [[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. -- Gregory (Greg) L. Snow Ph.D. 538...@gmail.com [[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] Problem with inconsolata font (again) --- on Fedora 17 this time.
Some while ago I posted a problem on this list concerning a failure of R CMD check on one of my packages that resulted from LaTeX being unable to find the inconsolata font. This was under the Ubuntu OS. This was solved, thanks to advice I got from this list, basically by doing sudo apt-get install texlive-fonts-extra I am now (for reasons which I won't go into) running Fedora 17, on a different computer. Now LaTeX cannot find the necessary inconsolata.sty file. I googled around a bit, and an item I found led me to check whether I actually had texlive installed on my (new) system. I hadn't!!! So I did sudo yum install texlive and that seemed to work. But then the item I'd found indicated that I should do sudo yum install texlive-inconsolata, texlive-inconsolata-font But I got: No package texlive-inconsolata, available. No package texlive-inconsolata-font available. I also tried sudo yum install texlive-fonts-extra with a similar result. Can anyone give me a simple recipe as to how to get the inconsolata font and the associated inconsolata.sty? Thanks. cheers, Rolf Turner __ 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] another R2HTML question, please
Dear R People: Is there a way to just print the commands without output into R2HTML, please? What I would like to do is to put up some commands for the students and see if they can get results. Thanks, Erin -- Erin Hodgess Associate Professor Department of Computer and Mathematical Sciences University of Houston - Downtown mailto: erinm.hodg...@gmail.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.
[R] Manual two-way demeaning of unbalanced panel data (Wansbeek/Kapteyn transformation)
Dear R users,
I wish to manually demean a panel over time and entities. I tried to code
the Wansbeek and Kapteyn (1989) transformation (from Baltagi's book Ch. 9).
As a benchmark I use both the pmodel.response() and model.matrix() functions
in package plm and the results from using dummy variables. As far as I
understood the transformation (Ch.3), Q%*%y (with y being the dependent
variable) should yield the demeaned series.
However, ...
...I find that the results do not match, if I do so.
...that if I am looking at a balanced panel, I get the correct
results
when multiplying Q with the already demeaned series y_it, Q%*%y_it.
...that if I am looking at an unbalanced panel, I receive
results which
differ (even though being close).
I guess I am missing something. Every comment pointing me to the right
solution would be of great value to me. Also, comments on how to increase
the efficiency of my function would help!
Please find an example based on the Grunfeld data below.
Thank you very much!
Philipp Grueber
##
library(MASS)
getQ-function(object,t.index,i.index){
t.ind-object[,t.index]
i.ind-object[,i.index]
nam-unique(i.ind)
tim-unique(t.ind)
n-nrow(object)
N-length(nam)
T-length(tim)
I-matrix(0,n,n)
I[row(I)==col(I)] - 1
I_N-matrix(0,N,N)
I_N[row(I_N)==col(I_N)] - 1
I_T-matrix(0,T,T)
I_T[row(I_T)==col(I_T)] - 1
D1-data.frame()
for(t in tim){
D1-rbind(D1,I_N)
}
D1-data.matrix(D1[as.vector(table(i.ind,t.ind))==1,])
D2-data.frame()
for(i in nam){
D2-rbind(D2,I_T)
}
D2-data.matrix(D2[as.vector(table(t.ind,i.ind))==1,])
D-data.matrix(cbind(D1,D2))
Q-I-D%*%ginv(crossprod(D))%*%t(D)#fQ(D1)-fQ(D1)%*%D2%*%ginv(t(D2)%*%fQ(D1)%*%D2)%*%t(D2)%*%fQ(D1)
Q
}
##
library(plm)
library(lmtest)
data(Grunfeld)
y_i-Grunfeld$inv-ave(Grunfeld$inv,index=Grunfeld$firm)
y_t-Grunfeld$inv-ave(Grunfeld$inv,index=Grunfeld$year)
y_it-(Grunfeld$inv-ave(Grunfeld$inv,index=Grunfeld$firm)-ave(Grunfeld$inv,index=Grunfeld$year)+rep(mean(Grunfeld$inv),length(Grunfeld$inv)))
x1_it-(Grunfeld$value-ave(Grunfeld$value,index=Grunfeld$firm)-ave(Grunfeld$value,index=Grunfeld$year)+rep(mean(Grunfeld$value),length(Grunfeld$inv)))
dem_y_i-pmodel.response(plm(formula=inv~value,data=Grunfeld,index=c(firm,year),model=within,effect=individual))
dem_y_t-pmodel.response(plm(formula=inv~value,data=Grunfeld,index=c(firm,year),model=within,effect=time))
dem_y_it-pmodel.response(plm(formula=inv~value,data=Grunfeld,index=c(firm,year),model=within,effect=twoways))
dem_X_it-model.matrix(plm(formula=inv~value,data=Grunfeld,index=c(firm,year),model=within,effect=twoways))
sum(y_i!=dem_y_i)
#y_i[1:10]
#dem_y_i[1:10]
sum(y_t!=dem_y_t)
#y_t[1:10]
#dem_y_t[1:10]
sum(y_it!=dem_y_it)
#y_it[1:10]
#dem_y_it[1:10]
sum(x1_it!=dem_X_it)
#x1_it[1:10]
#dem_X_it[1:10,]
(getQ(Grunfeld,t.index=year,i.index=firm)%*%Grunfeld$inv)[1:10]
(getQ(Grunfeld,t.index=year,i.index=firm)%*%y_it)[1:10]
dem_y_it[1:10]
Grunfeld1-Grunfeld[-1,]
sum(ave(Grunfeld1$inv,index=Grunfeld1$firm)!=(tapply(Grunfeld1$inv,Grunfeld1$firm,function(x){mean(x)})[Grunfeld1$firm]))==0
y_i-Grunfeld1$inv-ave(Grunfeld1$inv,index=Grunfeld1$firm)
y_t-Grunfeld1$inv-ave(Grunfeld1$inv,index=Grunfeld1$year)
y_it-(Grunfeld1$inv-ave(Grunfeld1$inv,index=Grunfeld1$firm)-ave(Grunfeld1$inv,index=Grunfeld1$year)+rep(mean(Grunfeld1$inv),length(Grunfeld1$inv)))
x1_it-(Grunfeld1$value-ave(Grunfeld1$value,index=Grunfeld1$firm)-ave(Grunfeld1$value,index=Grunfeld1$year)+rep(mean(Grunfeld1$value),length(Grunfeld1$inv)))
dem_y_i-pmodel.response(plm(formula=inv~value,data=Grunfeld1,index=c(firm,year),model=within,effect=individual))
dem_y_t-pmodel.response(plm(formula=inv~value,data=Grunfeld1,index=c(firm,year),model=within,effect=time))
dem_y_it-pmodel.response(plm(formula=inv~value,data=Grunfeld1,index=c(firm,year),model=within,effect=twoways))
dem_X_it-model.matrix(plm(formula=inv~value,data=Grunfeld1,index=c(firm,year),model=within,effect=twoways))
sum(y_i!=dem_y_i)
#y_i[1:10]
#dem_y_i[1:10]
sum(y_t!=dem_y_t)
#y_t[1:10]
#dem_y_t[1:10]
sum(y_it!=dem_y_it)
#y_it[1:10]
#dem_y_it[1:10]
#y_it[1:10]-dem_y_it[1:10]
sum(x1_it!=dem_X_it)
#x1_it[1:10]
#dem_X_it[1:10,]
#x1_it[1:10]-dem_X_it[1:10,]
(getQ(Grunfeld1,t.index=year,i.index=firm)%*%Grunfeld1$inv)[1:10]
(getQ(Grunfeld1,t.index=year,i.index=firm)%*%y_it)[1:10]
dem_y_it[1:10]
-
EBS Universitaet fuer Wirtschaft und Recht
FARE Department
Wiesbaden/ Germany
http://www.ebs.edu/index.php?id=finaccL=0
--
View this message in context:
http://r.789695.n4.nabble.com/Manual-two-way-demeaning-of-unbalanced-panel-data-Wansbeek-Kapteyn-transformation-tp4655202.html
Sent from the R
Re: [R] Sweave, Texshop, and sync with included Rnw file
On 13-01-10 4:54 PM, michele caseposta wrote:
Hi everybody,
thanks for the replies.
I might have not explained the problem completely.
Duncan Mackay:
Yes, I am already having a master file and separate Rnw files.
Duncan Murdock:
I am using patchDVI in the TexShop Sweave engine.
Sync works flawlessly between the master file and the pdf produced by pdflatex.
My problem is that I don't seem to be able to obtain sync between the
*included* Rnws and the pdf, either way.
I think you said before that you were using \SweaveInput to include the
files. I thought this had been handled, but perhaps not, or perhaps
there's a bug. What I'd recommend is that you don't use \SweaveInput.
Set up your main file main.Rnw like this:
... usual header stuff ...
echo=FALSE=
.SweaveFiles - c(chap1.Rnw, chap2.Rnw)
@
\input{chap1}
\input{chap2}
Set up the chapters like this. They aren't standalone files, so they
don't need the usual LaTeX header lines, but they'll need some Sweave stuff:
% Put a commented usepackage so Sweave doesn't insert one
%\usepackage{Sweave}
% Make sure to set a unique prefix in each chapter
% so that figures and concordances don't clash
\SweaveOpts{concordance=TRUE,prefix=chap1}
% Tell SweaveAll about the other files
echo=FALSE=
.SweaveFiles - main.Rnw
.TexRoot - main.tex
@
Then patchDVI::SweaveAll (or SweavePDF, etc.) on any chapter or on the
main.Rnw file will run Sweave on all the chapters (in a slightly
unpredictable order, so don't count on it). You should get the
concordances working for each file.
I'll take a look at what is happening with SweaveInput, but not tonight.
(Besides working, the setup described above has the advantage of
making things go faster: you won't need to run Sweave on unmodified
chapters. Only modified ones get run each time. You can also use
\include in place of \input, and then LaTeX will run faster with
\includeonly to select particular chapters.
Duncan Murdoch
The sweave engine is as follows:
#!/bin/bash
R CMD Sweave $1
latexmk -pdf -silent -pdflatex=‘pdflatex –shell-escape –synctex=1′${1%.*}
Rscript -e library(‘patchDVI’);patchSynctex(‘${1%.*}.synctex.gz’)
Funny thing is that the sync works in texworks, using the following Rscript line
patchDVI::SweavePDF('$fullname',stylepath=FALSE)
I tried to mix and match configurations between texshop and texworks but I had
no luck
On Jan 10, 2013, at 11:23 AM, Duncan Murdoch wrote:
On 13-01-09 9:09 PM, Duncan Murdoch wrote:
On 13-01-09 3:25 PM, michele caseposta wrote:
Hello everyone.
I am in the process of writing a book in Latex with Texshop, on Mac.
This book contains a lot of R code, hence the need to use Sweave.
I was able to compile Rnw files, and to sync back and forth from the pdf to the
source Rnw.
My problem now is that the book is divided in Chapters, and every chapter is in
its own Rnw file.
I can compile them from the main one (book.Rnw) using the directive
\SweaveInput{chapter1.Rnw}
The problem stands in the fact that like this I am missing synchronization
between the pdf and the source Rnw. If part of text is in book.Rnw I can
synchronize, but if the text is in one of the included files, it just doesn't
work.
I am using the sweave engine found in the following webpage:
http://cameron.bracken.bz/synctex-with-sweavepgfsweave-in-texshoptexworks
Has anybody succeeded in synchronizing with included Rnw files?
This is a problem addressed by my patchDVI package, available on
R-forge. You have a main file (which can be .tex or .Rnw), and put code
at the start of each .Rnw file to indicate where to find it. Then you
just run Sweave on one of the chapters, and it automatically produces
the full document.
The sample document here:
http://www.umanitoba.ca/statistics/seminars/2011/3/4/duncan-murdoch-using-sweave-R/
includes an appendix describing how to set this up with TeXShop.
I just committed an update to the vignette in patchDVI giving a quick version
of the instructions for basic use. Version 1.8.1585 has the new vignette.
I should get around to pushing it to CRAN one of these days...
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] Learning to speak R: simple data processing
So, I am just trying to learn R... Here is a rather contrived example that would be pretty obvious to me in terms of writing code to loop through elements, but the slick, fast, compact way of expressing this in R is not obvious to me. Here's code to generate a simple matrix of data: m - floor(matrix(runif(24, 1, 100), 8, 3)) m [,1] [,2] [,3] [1,] 755 15 [2,] 82 792 [3,] 72 24 55 [4,] 88 383 [5,] 42 82 98 [6,] 38 78 43 [7,] 79 88 60 [8,]6 89 43 I see how I can get the max on each row as: apply(m, 1, max) [1] 75 82 72 88 98 78 88 89 I want to generate a vector with the column position of each row max. The whole vector for the matrix shown above would be: 1 1 1 1 3 2 2 2 Is there some easy, straightforward way to compute such a thing short of looping over the rows and columns of the matrix? Thanks, -ej -- View this message in context: http://r.789695.n4.nabble.com/Learning-to-speak-R-simple-data-processing-tp4655194.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.
[R] polr model, out-of-sample probabilities
Hi, Is there a function to calculate probabilities for new out-of-sample data once we fit a model using the in-sample data? predict(model, newdata=... ) seems to require the new data to be the same size as the original data used to fit the model. In short, I would like to fit a model and then pass out-of-sample data to calculate probabilities. Regards, Alp [[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] Error with looping through a list of strings as variables
Dear R users:
I have been trying to figure out how to include string variables in a for
loop to run multiple random forests with little success. The current code
returns the following error:
Error in trafo(data = data, numeric_trafo = numeric_trafo, factor_trafo =
factor_trafo, :
data class character is not supported
In addition: Warning message:
In storage.mode(RET@predict_trafo) - double : NAs introduced by coercion
The code runs fine with the data before I add the for (h in varlist){ loop.
Loops i, k work without issue as long as I manually enter
the response variable into the code below for h.
Using R 2.15.0 (64bit), with cforest from the party package.
Any thoughts would be of great help.
Cheers
Note : The data used in the script below is not the actual data but
a substitute set which results in the same start-up errors. Once these
start-up errors are resolved there should be a (data, ...) : fraction of
0.00 is too small error will be seen which is simply due to the
small substitute data set and of no concern.
rm(list=ls())
library(party)
library(reshape)
puthere - c(TEST_RESULTS.csv)
hsb2 - read.csv(http://www.ats.ucla.edu/stat/data/hsb2.csv;)
names(hsb2)
set.seed(8296)
ctrl - cforest_unbiased(ntree=500, mtry=2)
varlist - names(hsb2)[3:4]
for (h in varlist){
for (k in c(1,0)){
for (i in c(1,2)){ ## Data subset
filtered - subset(hsb2,
schtyp == i
female == k,
select = c(id:socst))
rank.cf - cforest(h ~ write + math + science + socst,
data = filtered,
control = ctrl)
print(rank.cf)
## Standard importance values __
imp=varimp(rank.cf,
conditional = TRUE)
print(imp)
## predict variables _
predicted=predict(rank.cf,OOB = TRUE)
residual=filtered$h-predicted
mse=mean(residual^2)
rsq=1-mse/var(filtered$h)
##Correlation between fitted values and original values:
correl - paste(cor(filtered$h,predicted))
Correlation -paste(
MSE:,mse,
Rsq:,rsq,
Correlation between fitted values and original values:,correl)
print(Correlation)
## combine results for output ___
TestVar - paste(Dependent =,h, sep= )
namCL - paste(schtyp =,i, sep= )
namSE - paste(female =,k, sep= )
assign(namCL, 1:i)
assign(namSE, 1:k)
results - rbind(TestVar, namCL, namSE, mse, rsq, correl)
## Writing data to csv file _
write.table(results, file = puthere,
append = TRUE,
quote = FALSE,
sep = ,
col.names = TRUE,
row.names = TRUE,)
write.table(imp, file = puthere,
append = TRUE,
quote = FALSE,
sep = ,
eol = \r,
na = N/A,
row.names = TRUE,
col.names = TRUE,
qmethod = double)
}
}
}
[[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] how to generate a matrix by an my data.frame
Thanks a lot it works! 2013/1/11 Rui Barradas ruipbarra...@sapo.pt: Hello, Here are two ways. dat - read.table(text = id1id2 value 2353 2353 0.096313 2353 2409 0.301773 [...etc...] 2356 2356 0 2356 2611 0 2611 2611 0 , header = TRUE) mat1 - matrix(nrow = 53, ncol = 53) # initialize with NA's mat1[upper.tri(mat1, diag = TRUE)] - dat$value mat2 - matrix(0, nrow = 53, ncol = 53) # initialize with zeros mat2[upper.tri(mat2, diag = TRUE)] - dat$value Hope this helps, Rui Barradas Em 10-01-2013 15:21, Yao He escreveu: Dear All It is a little hard to give a good small example of my question,so I will show the full data on the bottom and the attachment.Maybe some one could tell me an appropriate way to show it.I'm sorry for the inconvenience. Q:How to generate a 53*53 diagonal matrix by my data Some problems confused me are that: 1.Since it is a diagonal matrix,I have tried to transform col1 and col2 to rowindex and colindex ,but I don't know how to generate matrix by its value's index 2. As you see, the number of 2353 corresponding to other ids in col2 is 53,however,the number of 2409 corresponding to other ids in col2 is 52 and 2500 corresponding to 51 values and so on,so it is hard to use matrix() to generate it id1id2 value 2353 23530.096313 2353 24090.301773 2353 25000.169518 2353 25980.11274 2353 26100.107414 2353 23000.034492 2353 25070.037521 2353 25300.064125 2353 23270.029259 2353 23890.036423 2353 24080.029259 2353 24630.036423 2353 24200.04409 2353 25630.055038 2353 24620.046478 2353 22920.036369 2353 24050.036369 2353 25430.053413 2353 25570.058151 2353 25830.081512 2353 23220.044373 2353 25350.04847 2353 25360.035538 2353 25810.035538 2353 25700.07711 2353 24760.047081 2353 25340.047081 2353 22800.088264 2353 23160.073608 2353 23390.067307 2353 23310.061172 2353 23430.060425 2353 23520.041153 2353 22930.040764 2353 23380.045128 2353 24490.040764 2353 22960.061333 2353 24530.046074 2353 24600.060387 2353 24740.060387 2353 26030.060387 2353 22820.048065 2353 23130.05584 2353 25380.050873 2353 25220.065727 2353 24890.041023 2353 25640.039696 2353 25940.056946 2353 22740.060875 2353 24510.037468 2353 23210 2353 23560 2353 26110 2409 24090.096313 2409 25000.169518 2409 25980.11274 2409 26100.107414 2409 23000.034492 2409 25070.037521 2409 25300.064125 2409 23270.029259 2409 23890.036423 2409 24080.029259 2409 24630.036423 2409 24200.04409 2409 25630.055038 2409 24620.046478 2409 22920.036369 2409 24050.036369 2409 25430.053413 2409 25570.058151 2409 25830.081512 2409 23220.044373 2409 25350.04847 2409 25360.035538 2409 25810.035538 2409 25700.07711 2409 24760.047081 2409 25340.047081 2409 22800.088264 2409 23160.073608 2409 23390.067307 2409 23310.061172 2409 23430.060425 2409 23520.041153 2409 22930.040764 2409 23380.045128 2409 24490.040764 2409 22960.061333 2409 24530.046074 2409 24600.060387 2409 24740.060387 2409 26030.060387 2409 22820.048065 2409 23130.05584 2409 25380.050873 2409 25220.065727 2409 24890.041023 2409 25640.039696 2409 25940.056946 2409 22740.060875 2409 24510.037468 2409 23210 2409 23560 2409 26110 2500 25000.048615 2500 25980.051979 2500 26100.041031 2500 23000.032974 2500 25070.052788 2500 25300.041435 2500 23270.038071 2500 23890.051659 2500 24080.038071 2500 24630.051659 2500 24200.052635 2500 25630.07872 2500 24620.048615 2500 22920.044365 2500 24050.044365 2500 25430.04277 2500 25570.051109 2500 25830.047409 2500 23220.054512 2500 25350.039368 2500 25360.041763 2500 25810.041763 2500 25700.063148 2500 24760.043755 2500 25340.043755 2500 22800.063164 2500 23160.083961 2500 23390.074127 2500 23310.051094 2500 23430.060066 2500 23520.038208 2500 22930.048698 2500 23380.048218 2500 24490.048698 2500 22960.073212 2500 24530.070595 2500 24600.073677 2500 24740.073677 2500 26030.073677 2500 22820.073677 2500 23130.068443 2500 25380.079865 2500 25220.059723 2500 24890.049873 2500 25640.087639 2500 25940.05175 2500 22740.043396 2500 24510.046532 2500 23210 2500 2356
Re: [R] Learning to speak R: simple data processing
?which.max On Fri, Jan 11, 2013 at 7:59 AM, ej ejohn...@earthlink.net wrote: apply(m, 1, max) __ 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] polr model, out-of-sample probabilities
On 11/01/2013 01:59, Alphan Kirayoglu wrote: Hi, Is there a function to calculate probabilities for new out-of-sample data once we fit a model using the in-sample data? predict(model, newdata=... ) seems to require the new data to be the same size as the original data used to fit the model. It does not, so why do you assume so? In short, I would like to fit a model and then pass out-of-sample data to calculate probabilities. That is what predict() is for. Your subject line says 'polr model'. Perhaps this was about my package MASS, but in any case please talk to the credit-where-credit-is-due department. And you could have seen this with a simple modification of the example on the help page for MASS::polr Regards, Alp [[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. -- Brian D. Ripley, rip...@stats.ox.ac.uk Professor of Applied Statistics, 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 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
