Re: [R] predict.drm not generating confidence intervals
I interpreted the to be that predict.drc was expecting a third
argument, curveid, which had no default, and that creating a
dataframe like this was going to solve the problem.
newdt - data.frame( conc= seq(0.2, 9, 0.01) , CURVE=1)
prd.p - predict(fit, newdata=newdt, curveid=CURVE,
interval=confidence)
head(prd.p)
PredictionLowerUpper
[1,] 7.790812 7.400131 8.181492
[2,] 7.790476 7.400161 8.180791
[3,] 7.790106 7.400185 8.180028
[4,] 7.789702 7.400203 8.179202
[5,] 7.789262 7.400213 8.178311
[6,] 7.788784 7.400215 8.177352
--
David
On Nov 28, 2010, at 5:19 PM, Peter Ehlers wrote:
Brant,
See below.
On 2010-11-28 12:25, David Winsemius wrote:
Puzzled. Why are the data you offer to predict() for the independent
variable, conc, all NA's? Is there something reversed or inverted
about how drc functions handle formulas.
-- David. On Nov 28, 2010, at 2:33 PM, Brant Inman wrote:
R-helpers,
I recently submitted a help request for the predict.drm function
found in the drc package. I am still having issues with the
function and I am submitting reproducible code hoping that
somebody
can help me figure out what is going on.
library(drc)
# Fit a 4 parameter logistic model to ryegrass dataset
fit- drm(rootl ~ conc, data = ryegrass, fct = LL.4())
summary(fit)
# Generate a fake dataset for prediction
newdt- data.frame( matrix(c(seq(0.2, 9, 0.01), rep(NA,881)),
ncol=2))
colnames(newdt)- c('rootl', 'conc')
# Generate prediction intervals and confidence intervals
prd.p- predict(fit, newdata=newdt, interval='prediction')
prd.c- predict(fit, newdata=newdt, interval='confidence')
# Check output
head(prd.p)
Prediction Lower Upper
[1,] 7.790812NANA
[2,] 7.790476NANA
[3,] 7.790106NANA
[4,] 7.789702NANA
[5,] 7.789262NANA
[6,] 7.788784NANA
head(prd.c)
Prediction Lower Upper
[1,] 7.790812NANA
[2,] 7.790476NANA
[3,] 7.790106NANA
[4,] 7.789702NANA
[5,] 7.789262NANA
[6,] 7.788784NANA
There appears to be a problem with the predict.drc function.
This
code previous generated confidence and prediction intervals in
columns 2 and 3 of the prediction matrices but now fails for
reasons
unknown to me. Anyone have an idea of what is going on and how I
can remedy the situation?
As David points out, you seem to have your dose and response
mixed up. But even assuming that you mean to predict rootl
from conc, the predict.drc function needs to be used in a
certain way (and that is not clear in the documentation).
Try this:
fit - drm(rootl ~ conc, data = ryegrass, fct = LL.4())
newdat - data.frame(conc = seq(0, 30, len=101), apples=1)
ypred - predict(fit, newdata=newdat, interval=confidence)
head(ypred, 3)
PredictionLowerUpper
[1,] 7.792958 7.399614 8.186303
[2,] 7.785771 7.400046 8.171495
[3,] 7.736545 7.380584 8.092505
Wondering about the 'apples'? That's just to show that
predict.drc doesn't care what you call your second column.
It seems that predict.drc requires either *no* newdata
or a *two*-variable data.frame in which the first column
(it, too, can have *any* name) contains the dose at which
to predict and the second column contains the 'curveid',
aka the grouping variable (if you have one). In your
case, the grouping reduces to a single group, i.e. a
vector of all '1's.
This could all be made considerably clearer in the documentation.
Peter Ehlers
Brant Inman
David Winsemius, MD
West Hartford, CT
__
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] Array help
On Nov 28, 2010, at 10:56 PM, bfhancock wrote: Josh, the data set is called StatTemps and is in the PASWR package. I want to make an array that involves only the 8 a.m. and a separate array that involves only the 9 a.m. so i can get info on the temperatures in those groups. So I still want it in the format of StatTemps but in two arrays that are based on 8 a.m. or 9 a.m.I have been messing with this for a while. And no it's not homework, I am just trying to learn R so I am more appealing out in the field eventually. They book I am using is confusing! Hopefully what I am trying to do isn't confusing. I want to do an array that holds the info from 1:6 12:22 for 8am and then 7:11 23:34 for 9am. I was able easily make two arrays based on sex by just putting in StatTemps[1:11,,] StatTemps[12:34,,] and assumed I could just go StatTemps[1:612:22,,] and StatTemps[7:1123:34,,] but that didn't work. Any ideas? Thanks so much! It may not be an array (since time values don't play very well with that data structure.) If there is time as an index, it may be a more complex object such as a time-series or zoo. Check with str(). -- David Winsemius, MD West Hartford, CT __ 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] Plot data inside matrix
On Nov 29, 2010, at 6:22 AM, alcesgabbo wrote: Hi, I have this problem: I have this matrix: Doubtful that is is a matrix. In R matrices are all of the same type of object. This looks more like a zoo object since it has a time index. How was it created and what does str() show? result property procProperty 2010-10-01 07:32:00 40 Asensor1 2010-10-01 17:32:00 15 Asensor3 2010-10-02 07:32:00 32 Asensor2 2010-10-03 04:33:21 20 Bsensor1 2010-10-03 04:33:21 33 Bsensor2 2010-10-03 14:33:21 12 Asensor3 2010-10-05 07:32:00 31 Bsensor1 2010-10-05 07:32:00 15 Bsensor2 2010-10-06 17:32:00 4Asensor3 I would like to plot this matrix in this way: create in this case 2 plots (one for each property: A and B ) for each plot there will be 3 lines (one for each procProperty: sensor1,sensor2,sensor3) composed by the result. So, what do you want: a dotplot, a barchart , a time-series or what? How can I do this with few commands?? Possibly depending on the correctness of my assumptions: xyplot.zoo -- David Winsemius, MD West Hartford, CT __ 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] selecting only corresponding categories from a confusion matrix
On Nov 29, 2010, at 8:32 AM, drflxms wrote: Dear R colleagues, as a result of my calculations regarding the inter-observer- variability in bronchoscopy, I get a confusion matrix like the following: 0 1 1001 1010 11 0609 11 54 36 6 1 1 260 2 1014 008 4 1004 000 0 1000 23 7 12 10 5 1001 0 040 0 1010 4 003 0 1011 1 010 2 11 0 033 1 1101 000 0 1100 2 000 0 1110 1 000 0 The first column represents the categories found among observers, the top row represents the categories found by the reference (goldstandard). I am looking for a way (general algorithm) to extract a data.frame with only the corresponding categories among observers and reference from the above confusion matrix. Corresponding means in this case, that a category has been chosen by both: observers and reference. In this example corresponding categories would be simply all categories that have been chosen by the reference (0,1,1001,1010,11), but generally there might also occur categories which are found by the reference only (and not among observers - in the first column). So the solution-dataframe for the above example would look like: 0 1 1001 1010 11 0609 11 54 36 6 1 1 260 2 1001 0 040 0 1010 4 003 0 11 0 033 1 I wasn't able to follow the confusing, er, confusion matrix explanation but it appears from a comparison of the input and output that you just want row indices that are the column names: mtx[colnames(mtx), ] 0 1 1001 1010 11 0609 11 54 36 6 1 1 260 2 1001 0 040 0 1010 4 003 0 11 0 033 1 # and the omitted mtx[!rownames(mtx) %in% colnames(mtx), ] 0 1 1001 1010 11 10 14 008 4 100 4 000 0 1000 23 7 12 10 5 1011 1 010 2 110 1 000 0 1100 2 000 0 1110 1 000 0 # and their number: NROW(mtx[!rownames(mtx) %in% colnames(mtx), ]) [1] 7 All the categories found among observers only, were omitted. If the solution algorithm would include a method to list the omitted categories and to count their number as well as the number of omitted cases, it would be just perfect for me. David Winsemius, MD West Hartford, CT __ 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] Troubles in plotting to a postscript file (not to png)
On Nov 29, 2010, at 9:00 AM, pilchat wrote:
Hi guys,
to make it easier, here is a simple case with the same issues. I use
the short function below to make the attached PS file.
Things to fix:
-) the greek letter lambda is not printed, while mu is printed
(see the plot command)
-) the annotation inside the plot area: the +- symbol and (K)
overlap with the substitute for tmed and tstd respectively (see the
text command). Also, how can I set the number of decimal digits for
tmed and tstd? (option(digits=4) does not work )
I would have thought one would do any formatting (of digits) outside
the text( ...substitute(...),...) setting.
Moreover, I'd like to make the characters thicker. Is there any way?
Which characters? There is a bold() option within plotmath.
Finally, once I close the R session without saving it (I answer n
when quitting), the content of the ps file is erased.
Now _that_ is weird. A file should have been created in your default
directory and closing R should not have made it go away.
Do I miss something in writing the function?
Perhaps. (But you certainly missed something in writing the question.)
When I remove the family=ComputerModern from the postscript call, I
start seeing lambda. And the other spacing weirness also
resolves. I am on a Mac and ComputerModern is not one of the
pdfFonts() on my machine. The list of available fonts varies widely
across various OSes and devices about which you have given us no clues.
--
David.
Thanks
Gaetano
plot_example=function()
{
setPS()
postscript (file='plot_example.ps',width=5,height=5,horizontal =
FALSE, paper = special,family =
ComputerModern,encoding=TeXtext.enc)
tmed-1.23456789
tstd-1.23456789
plot
(c
(0,1
),c
(0,1
),xlab=expression(paste(lambda,mu,T)),main=,sub=(a))#lambda
not printed
text(.0,.8,substitute( T[disk,med] == tmed %+-% tstd
(K),list(tmed=tmed,tstd=tstd)),pos=4,cex=1.5 )#overlapping symbols
and numbers
dev.off()
}
plot_example.ps__
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, MD
West Hartford, CT
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] data.frame and formula classes of aggregate
On Nov 29, 2010, at 9:35 AM, David Freedman wrote:
Hi - I apologize for the 2nd post, but I think my question from a
few weeks
ago may have been overlooked on a Friday afternoon.
I might be missing something very obvious, but is it widely known
that the
aggregate function handles missing values differently depending if a
data
frame or a formula is the first argument ?
I'm not sure if it is widely known, but it is certainly suggested by
the documentation for aggregate, since aggregate.data.frame has
different defaults than aggregate.formula. See the Usage section at
the very top of ?aggregate.
For example,
(d- data.frame(sex=rep(0:1,each=3),
wt=c(100,110,120,200,210,NA),ht=c(10,20,NA,30,40,50)))
x1- aggregate(d, by = list(d$sex), FUN = mean);
names(x1)[3:4]- c('mean.dfcl.wt','mean.dfcl.ht')
x2- aggregate(cbind(wt,ht)~sex,FUN=mean,data=d);
names(x2)[2:3]- c('mean.formcl.wt','mean.formcl.ht')
cbind(x1,x2)[,c(2,3,6,4,7)]
The output from the data.frame class has an NA if there are missing
values
in the group for the variable with missing values. But, the formula
class
output seems to delete the entire row (missing and non-missing
values) if
there are any NAs. Wouldn't one expect that the 2 forms (data frame
vs
formula) of aggregate would give the same result?
thanks very much
david freedman, atlanta
--
David Winsemius, MD
West Hartford, CT
__
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] FW: how to use by() ?
On Nov 29, 2010, at 1:36 PM, Jim Moon wrote:
Thank you for the suggestion, Bill. The result is not quite what I
would like. Here's sample code for you or anyone else who may be
interested:
Al1 = c('A','C','C','C')
Al2 = c('G','G','G','T')
Freq1 = c(0.0078,0.0567,0.9434,0.9908)
MAF = c(0.0078,0.0567,0.0566,0.0092)
m1 = data.frame(Al1=Al1, Al2=Al2,Freq1=Freq1,MAF=MAF,major_allele='')
m1
Al1 Al2 Freq1MAF major_allele
1 A G 0.0078 0.0078
2 C G 0.0567 0.0567
3 C G 0.9434 0.0566
4 C T 0.9908 0.0092
Using the suggestion involving with() (I swapped Al1 and Al2 from
before, but this does not affect the nature of the output):
m1$major_allele - with(m1, ifelse(Freq1==MAF, Al2, Al1));m1
Al1 Al2 Freq1MAF major_allele
I suspect that you have just been bitten by the data.frame-
stringsAsFactors=TRUE crocodile. Since you are comparing floating
point numbers, you are also wading in rivers where floating-point
crocodiles are also hungry and searching out their next victim.
--
David.
1 A G 0.0078 0.00781
2 C G 0.0567 0.05671
3 C G 0.9434 0.05662
4 C T 0.9908 0.00922
The output I desire is:
Al1 Al2 Freq1MAF major_allele
1 A G 0.0078 0.0078G
2 C G 0.0567 0.0567G
3 C G 0.9434 0.0566C
4 C T 0.9908 0.0092C
Jim
-Original Message-
From: William Dunlap [mailto:wdun...@tibco.com]
Sent: Monday, November 29, 2010 10:02 AM
To: Jim Moon
Subject: RE: [R] how to use by() ?
m1$major_allele - with(m1, ifelse(Freq1==MAF, Al1, Al2))
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 Jim Moon
Sent: Monday, November 29, 2010 9:44 AM
To: r-help@r-project.org
Subject: [R] how to use by() ?
Hello, All!
How might one accomplish this using the by() function?
m1 is a data frame.
# populate column m1$major_allele
for ( i in 1:length(m1$major_allele)) {
if ( m1$Freq1[i] == m1$MAF[i]){
m1$major_allele[i] = m1$Al1[i]
}
else{
m1$major_allele[i] = m1$Al2[i]
}
}
Jim
[[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.
David Winsemius, MD
West Hartford, CT
__
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] reference text variables as column name to plot
On Nov 30, 2010, at 11:40 AM, Graves, Gregory wrote:
Having given myself carpal tunnel looking for answer to this ...
I have a dataset
Named what?
each column of which has 12 rows in it. I created a
variable 'z' as follows:
z=1:24 s
Why did you put twice as many elements in z as there are in a column?
Since I have a large number of these plots to make, and they are a bit
complex, I want to want to reference the column I want to plot via a
variable containing the name of that column. As follows:
similar='1978'
s=paste('Y',similar,sep='')
variable s now contains 'Y1978' which is the name of one of the
columns.
However, when I try to plot
plot(z,s,type='l')
Try this ... assuming that there really are 12 item length columns in
a dataframe named, dfm.
plot(1:12, dfm[[s]], type=l)
dataframes are lists that can be accessed by the names of columns
which are interpreted. Don't assume that you can get such
interpretation with the $ operator.
I get a 'x and y lengths differ' error because variable s is being
recognized as 'Y1978' length=1, rather than the contents of the column
Y1978 length=12.
I tried all the usual tricks I know like s.
Huh? is a logical operator.
How do you get R to
reference a variable as a column name?
Thank you.
Gregory A. Graves, Lead Scientist
Everglades REstoration COoordination and VERification (RECOVER)
Restoration Sciences Department
South Florida Water Management District
Phones: DESK: 561 / 682 - 2429
CELL: 561 / 719 - 8157
[[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, MD
West Hartford, CT
__
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] And about sum() ...Re: Minor warning about seq
On Nov 30, 2010, at 1:49 PM, Prof. John C Nash wrote: I spent more time than I should have debugging a script because I wanted x-seq(0,100)*0.1 but typed x-seq(O:100)*0.1 seq(0:100) yields 1 to 101, Clearly my own brain to fingers fumble, but possibly one others may want to avoid it. I discovered some flakey code (that ran without syntactic error and sometimes even without semantic error) that was intended to add a risk estimate from a number of completed years of exposure to a fractional year: sum(q[1:4] + frac*q[5]) ... looked to my brain as sensible, but was actually adding the fractional risk to each of the q[1:4]'s before summation occurred. (Should have been ... sum(q[1:4] , frac*q[5]) -- David Winsemius, MD West Hartford, CT __ 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] warning creating an as.array method in a package
On Nov 30, 2010, at 5:43 PM, Michael Friendly wrote: [Env: R 2.11.1, Win Xp, using Eclipse/StatET] In a package I'm working on, I want to create as.matrix() and as.array() methods for a particular kind of object (log odds ratios). These are returned in a loddsratio object as the $coefficients component, a vector, but really reflect an underlying (R-1)x(C-1)xstrata array, whose attributes are contained in other components. I define coef, dim and dimnames methods, and then want an as.array method, ## dim methods dimnames.loddsratio - function(object, ...) object$dimnames dim.loddsratio - function(object, ...) object$dim coef.loddsratio - function(object, log = object$log, ...) if(log) object$coefficients else exp(object$coefficients) as.array.loddsratio - function (x, log=x$log, ...) drop(array(coef(x, log = log), dim = dim(x), dimnames=dimnames(x))) I believe everything is declared properly in the NAMESPACE file, e.g., ... S3method(dim, loddsratio) S3method(dimnames, loddsratio) S3method(print, loddsratio) S3method(vcov, loddsratio) S3method(as.matrix, loddsratio) S3method(as.array, loddsratio) All my tests work correctly when run in the R console, but R CMD check gives me a perplexing warning: * checking whether package 'vcdExtra' can be installed ... WARNING Found the following significant warnings: Warning: found an S4 version of 'as.array' so it has not been imported correctly See 'C:/eclipse/vcdExtra.Rcheck/00install.out' for details. It's just a warning, so maybe I can ignore it, but I can't figure out where this might have come from. Has anyone seen this before? Where might an S4 version of as.array be found? Have you looked in the .out file in which you were told there were details? I'm guessing from this behavior on my system that it may be Matrix. showMethods(as.array) Function: as.array (package base) x=ANY x=Matrix -- David Winsemius, MD West Hartford, CT __ 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] Can't Destroy Dim Names
On Nov 30, 2010, at 6:07 PM, clangkamp wrote: Dear Jim I think the target is to get from a Named chr to a just chr str(mat) Named chr [1:32268] yQAAA jQAAQ UQAAg FQAAw 1QABA ... - attr(*, names)= chr [1:32268] CA CC CG CT ... I have presumably the same problem str(DC1a) num [1:18, 1:48, 1:35] 3124.4 3049.2 227.8 41.4 76 ... - attr(*, dimnames)=List of 3 ..$ Figure : Named chr [1:18] CDS1 ... .. ..- attr(*, names)= chr [1:18] 1 ..$ Code: Named chr [1:48] AGR .. ..- attr(*, names)= chr [1:48] 1 36 71 106 ... ..$ variable: Named chr [1:35] X30.09.2009 .. ..- attr(*, names)= chr [1:35] 1 2 3 4 ... DC1_SM-abind(DC1a, DC1_PLCF_SM1, along=1, new.names=) str(DC1_SM) num [1:24, 1:48, 1:35] 3124.4 3049.2 227.8 41.4 76 ... - attr(*, dimnames)=List of 3 ..$ : chr [1:24] CDS1 ... ..$ : chr [1:48] AGR ..$ : chr [1:35] X30.09.2009 names(dimnames(DC1_PLCF_SM1))-names(dimnames(DC1a)) The point is to kill the lines with the bit .. ..- attr(*, names)= chr [1:35] 1 2 3 4 ... and change the Named chr into a plain chr. It is not at all clear to me that the problem posed a year and a half ago is the same as the one you perceive you are facing. In any event you are welcome to mangle your object (which you have not offered for testing) by turning a named dimension name vector into an unnamed one: ?unname ?Extract DCtest - array(1:27, c(3,3,3)) dimnames(DCtest) - list(dim1 =c(a=a,b=b,c=c), #named vector dim2=letters[4:6],#unnamed vectors dim3= letters[7:9]) str(DCtest) int [1:3, 1:3, 1:3] 1 2 3 4 5 6 7 8 9 10 ... - attr(*, dimnames)=List of 3 ..$ dim1: Named chr [1:3] a b c .. ..- attr(*, names)= chr [1:3] a b c ..$ dim2: chr [1:3] d e f ..$ dim3: chr [1:3] g h i dimnames(DCtest)[1] $dim1 a b c a b c dimnames(DCtest)[[1]] a b c a b c So use the [[- function to replace the named vector with an unnamed one: dimnames(DCtest)[[1]] - unname( dimnames(DCtest)[[1]] ) str(DCtest) int [1:3, 1:3, 1:3] 1 2 3 4 5 6 7 8 9 10 ... - attr(*, dimnames)=List of 3 ..$ dim1: chr [1:3] a b c ..$ dim2: chr [1:3] d e f ..$ dim3: chr [1:3] g h i - Christian Langkamp christian.langkamp-at-gmxpro.de -- View this message in context: http://r.789695.n4.nabble.com/Can-t-Destroy-Dim-Names-tp876633p3066413.html David Winsemius, MD West Hartford, CT __ 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] Pass an operator to function
On Nov 30, 2010, at 9:54 PM, randomcz wrote:
Hi guys,
How to pass an operator to a function. For example,
test - function(a, , b)
{
return(ab) #the operator is passed as an argument
}
I think you have just requested the definition of do.call() although
you infix positioning is a bit non-standard:
?do.call
do_this - function(a, fn=, b) {do.call(fn, list(a , b))}
do_this(a=1, b=4)
[1] FALSE
do_this(a=1, b=0)
[1] TRUE
--
David Winsemius, MD
West Hartford, CT
__
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 mean, min lists and numbers
On Jul 12, 2010, at 11:10 AM, g...@ucalgary.ca wrote:
I would like to sum/mean/min a list of lists and numbers to return the
related lists.
You will advance in your understanding faster if you adopt the correct
terminology:
-1+2*c(1,1,0)+2+c(-1,10,-1) returns c(2,13,0) but
... which is NOT a list, it is a vector.
sum(1,2*c(1,1,0),2,c(-1,10,-1)) returns 15 not a list.
Using the suggestions of Gabor Grothendieck,
Reduce('+',list(-1,2*c(1,1,0),2,c(-1,10,-1))) returns what we want,
c(2,13,0).
However, it seems that this way does not work to mean/min.
If you want a running cumulative mean of a vector, i.e,
c( mean(vec[1]), mean(vec[1:2]), ,,, mean(vec) ):
vec - sample(1:20)
sapply(1:length(vec), function(x) mean(vec[1:x])
So, how to mean/min a list of lists and numbers to return a list?
Not a list and not working on a list of lists. A vector.
--
David Winsemius, MD
West Hartford, CT
__
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 mean, min lists and numbers
On Jul 12, 2010, at 11:19 AM, Duncan Murdoch wrote:
On 12/07/2010 11:10 AM, g...@ucalgary.ca wrote:
I would like to sum/mean/min a list of lists and numbers to return
the
related lists.
-1+2*c(1,1,0)+2+c(-1,10,-1) returns c(2,13,0) but
sum(1,2*c(1,1,0),2,c(-1,10,-1)) returns 15 not a list.
Using the suggestions of Gabor Grothendieck,
Reduce('+',list(-1,2*c(1,1,0),2,c(-1,10,-1))) returns what we want,
c(2,13,0).
However, it seems that this way does not work to mean/min.
So, how to mean/min a list of lists and numbers to return a list?
Thanks,
You need to be careful of terminology: c(1,1,0) is not a list, it's
a vector. What you want is to apply functions componentwise to
lists of vectors.
One way to do that is to bind them into a matrix, and use apply.
For example:
M - cbind(-1, c(1,1,0), c(-1,10,-1))
apply(M, 1, mean)
As usual Duncan's understanding is better than mine. Just so you know,
there are also utility functions row-oriented functions which are
conisderably faster when they are the correct solution:
?rowSums
?rowMeans
rowMeans(cbind(-1, c(1,1,0), c(-1,10,-1)) )
[1] -0.333 3.333 -0.667
apply(cbind(-1, c(1,1,0), c(-1,10,-1)), 1, mean)
[1] -0.333 3.333 -0.667
... and there is a parallel version of a minimum functions, pmin
that would have given you results for arguments that are just the
vectors of varying length you were working with:
pmin(2, c(2,2,0) ,-1 , c(-1,10,-1))
#[1] -1 -1 -1 #(Done with argument recycling.)
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.
David Winsemius, MD
West Hartford, CT
__
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] in continuation with the earlier R puzzle
On Jul 12, 2010, at 10:09 AM, Raghu wrote:
When I just run a for loop it works. But if I am going to run a for
loop
every time for large vectors I might as well use C or any other
language.
The reason R is powerful is becasue it can handle large vectors
without each
element being manipulated? Please let me know where I am wrong.
for(i in 1:length(news1o)){
+ if(news1o[i]s2o[i])
+ s[i]-1
+ else
+ s[i]--1
+ }
Perhaps:
s - 2*( news1o s2o[1:length(news1o)] ) - 1
...which I think will throw errors under pretty much the same
conditions that would cause errors in that loop.
--
David Winsemius, MD
West Hartford, CT
__
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 select the column header with \Sexpr{}
On Jul 12, 2010, at 5:45 PM, Felipe Carrillo wrote:
Thanks for the quick reply Duncan.
I don't think I have explained myself well, I have a dataset named
report and
my column headers are run1,run2,run3,run4 and so on.
I know how to access the data below those columns with
\Sexpr{report[1,1]}
\Sexpr{report[1,2]} and so on, but I can't access my column headers
with \Sexpr{} because I can't find the way to reference
run1,run2,run3 and run4.
Sorry if I am not explain myself really well.
Wouldn't this just be:
\Sexpr{names(report)} # ? or perhaps you want specific items in that
vector?
Sexpr{names(report)[1]}, Sexpr{names(report)[2]}, etc
--
David.
- Original Message
From: Duncan Murdoch murdoch.dun...@gmail.com
To: Felipe Carrillo mazatlanmex...@yahoo.com
Cc: r-h...@stat.math.ethz.ch
Sent: Mon, July 12, 2010 2:18:15 PM
Subject: Re: [R] How to select the column header with \Sexpr{}
On 12/07/2010 5:10 PM, Felipe Carrillo wrote:
Hi:
Since I work with a few different fish runs my column headers change
everytime
I start a new Year. I have been using \Sexpr{} for my row and
columns and now
I am trying to use with my report column headers. \Sexpr{1,1} is
row 1 column 1,
what can I use for headers? I tried \Sexpr{0,1} but sweave didn't
like
it..Thanks in advance
for any hints
\Sexpr takes an R expression, and inserts the first element of the
result into
your text. Using just 0,1 (not including the quotes) is not a
valid R
expression.
You need to use paste() or some other function to construct the
label you want
to put in place, e.g. \Sexpr{paste(0,1,sep=,)} will give you 0,1.
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.
David Winsemius, MD
West Hartford, CT
__
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] Comparison of two very large strings
On Jul 12, 2010, at 6:03 PM, harsh yadav wrote:
Hi,
I have a function in R that compares two very large strings for
about 1
million records.
The strings are very large URLs like:-
http://query.nytimes.com/gst/sitesearch_selector.html?query=US+Visa+Lawstype=nytx=25y=8
.
..
or of larger lengths.
The data-frame looks like:-
id url
1
http://query.nytimes.com/gst/sitesearch_selector.html?query=US+Visa+Lawstype=nytx=25y=8
.
..
2 http://query.nytimes.com/search/sitesearch?query=US+Visa+Lawssrchst=cse
3
http://www.google.com/search?hl=enq=us+student+visa+changes+9/11+washington+poststart=10sa=N
.
..
4
http://www.google.com/search?hl=enq=us+student+visa+changes+9/11+washington+poststart=10sa=N
5
http://www.google.com/url?sa=Ustart=11q=http://app1.chinadaily.com.cn/star/2004/0610/fo4-1.htmlei=uUKwSe7XN9CCt
and so on for about 1 million records.
Here is the function that I am using to compare the two strings:-
stringCompare - function(currentURL, currentId){
j - currentId - 1
while(j=1)
previousURL -
previousURLLength - nchar(previousURL)
#Compare smaller with bigger
if(nchar(currentURL) = previousURLLength){
matchPhrase - substr(previousURL,1,nchar(currentURL))
if(matchPhrase == currentURL){
return(TRUE)
}
}else{
matchPhrase - substr(currentURL,1,previousURLLength)
if(matchPhrase == previousURL){
return(TRUE)
}
}
j - j -1
}
return(FALSE)
}
Couldn't you just store the url vector after running through nchar
and then do the comparison in a vectorized manner?
test - rd.txt('id url
1 http://query.nytimes.com/gst/sitesearch_selector.html?query=US+Visa+Lawstype=nytx=25y=8
2 http://query.nytimes.com/search/sitesearch?query=US+Visa+Lawssrchst=cse
3 http://www.google.com/search?hl=enq=us+student+visa+changes+9/11+washington+poststart=10sa=N
4 http://www.google.com/search?hl=enq=us+student+visa+changes+9/11+washington+poststart=10sa=N
5 http://www.google.com/url?sa=Ustart=11q=http://app1.chinadaily.com.cn/star/2004/0610/fo4-1.htmlei=uUKwSe7XN9CCt
', stringsAsFactors=FALSE)
copyUrls - test[,url]
sizeUrls - nchar(copyUrls)
lengU - length(sizeUrls)
sizidx - pmax(sizeUrls[1:(lengU-1)], sizeUrls[2:(lengU)])
substr(copyUrls[2:lengU], 1, sizidx) == substr(copyUrls[1:(lengU-1)],
1, sizidx)
#[1] FALSE FALSE TRUE FALSE
Here, I compare the URL at a given row with all the previous URLs in
the
data-frame. I compare the smaller of the two given URls with the
larger one
(upto the length of the smaller).
When I run the above function for about 1 million records, the
execution
becomes really slow, which otherwise is fast if I remove the
string comparison step.
Any ideas how it can be implemented in a fast and efficient way.
Thanks and Regards,
Harsh Yadav
[[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, MD
West Hartford, CT
__
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] Comparison of two very large strings
On Jul 12, 2010, at 6:46 PM, David Winsemius wrote:
On Jul 12, 2010, at 6:03 PM, harsh yadav wrote:
Hi,
I have a function in R that compares two very large strings for
about 1
million records.
The strings are very large URLs like:-
http://query.nytimes.com/gst/sitesearch_selector.html?query=US+Visa+Lawstype=nytx=25y=8
.
..
or of larger lengths.
The data-frame looks like:-
id url
1
http://query.nytimes.com/gst/sitesearch_selector.html?query=US+Visa+Lawstype=nytx=25y=8
.
..
2 http://query.nytimes.com/search/sitesearch?query=US+Visa+Lawssrchst=cse
3
http://www.google.com/search?hl=enq=us+student+visa+changes+9/11+washington+poststart=10sa=N
.
..
4
http://www.google.com/search?hl=enq=us+student+visa+changes+9/11+washington+poststart=10sa=N
5
http://www.google.com/url?sa=Ustart=11q=http://app1.chinadaily.com.cn/star/2004/0610/fo4-1.htmlei=uUKwSe7XN9CCt
and so on for about 1 million records.
Here is the function that I am using to compare the two strings:-
stringCompare - function(currentURL, currentId){
j - currentId - 1
while(j=1)
previousURL -
previousURLLength - nchar(previousURL)
#Compare smaller with bigger
if(nchar(currentURL) = previousURLLength){
matchPhrase - substr(previousURL,1,nchar(currentURL))
if(matchPhrase == currentURL){
return(TRUE)
}
}else{
matchPhrase - substr(currentURL,1,previousURLLength)
if(matchPhrase == previousURL){
return(TRUE)
}
}
j - j -1
}
return(FALSE)
}
Couldn't you just store the url vector after running through
nchar and then do the comparison in a vectorized manner?
test - rd.txt('id url
1 http://query.nytimes.com/gst/sitesearch_selector.html?query=US+Visa+Lawstype=nytx=25y=8
2 http://query.nytimes.com/search/sitesearch?query=US+Visa+Lawssrchst=cse
3 http://www.google.com/search?hl=enq=us+student+visa+changes+9/11+washington+poststart=10sa=N
4 http://www.google.com/search?hl=enq=us+student+visa+changes+9/11+washington+poststart=10sa=N
5 http://www.google.com/url?sa=Ustart=11q=http://app1.chinadaily.com.cn/star/2004/0610/fo4-1.htmlei=uUKwSe7XN9CCt
', stringsAsFactors=FALSE)
copyUrls - test[,url]
sizeUrls - nchar(copyUrls)
lengU - length(sizeUrls)
sizidx - pmax(sizeUrls[1:(lengU-1)], sizeUrls[2:(lengU)])
substr(copyUrls[2:lengU], 1, sizidx) == substr(copyUrls[1:
(lengU-1)], 1, sizidx)
#[1] FALSE FALSE TRUE FALSE
Let me hasten to admit that when I tried to fix what I thought was an
error in that program, I go the same result. It seemed as though I
should have been getting errors by choosing the maximum string
length. Changing the pmax to pmin did not alter the results ... to my
puzzlement ... until I further noticed that urls #3 and #4 were of the
same length. When I extend the lengths, then only the version using
pmin works properly.
--
David.
Here, I compare the URL at a given row with all the previous URLs
in the
data-frame. I compare the smaller of the two given URls with the
larger one
(upto the length of the smaller).
When I run the above function for about 1 million records, the
execution
becomes really slow, which otherwise is fast if I remove the
string comparison step.
Any ideas how it can be implemented in a fast and efficient way.
Thanks and Regards,
Harsh Yadav
[[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, MD
West Hartford, CT
__
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, MD
West Hartford, CT
__
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] Continuing on with a loop when there's a failure
On Jul 12, 2010, at 6:18 PM, Josh B wrote:
Hi R sages,
Here is my latest problem. Consider the following toy example:
x - read.table(textConnection(y1 y2 y3 x1 x2
indv.1 bagels donuts bagels 4 6
indv.2 donuts donuts donuts 5 1
indv.3 donuts donuts donuts 1 10
indv.4 donuts donuts donuts 10 9
indv.5 bagels donuts bagels 0 2
indv.6 bagels donuts bagels 2 9
indv.7 bagels donuts bagels 8 5
indv.8 bagels donuts bagels 4 1
indv.9 donuts donuts donuts 3 3
indv.10 bagels donuts bagels 5 9
indv.11 bagels donuts bagels 9 10
indv.12 bagels donuts bagels 3 1
indv.13 donuts donuts donuts 7 10
indv.14 bagels donuts bagels 2 10
indv.15 bagels donuts bagels 9 6), header = TRUE)
I want to fit a logistic regression of y1 on x1 and x2. Then I want
to run a
logistic regression of y2 on x1 and x2. Then I want to run a
logistic regression
of y3 on x1 and x2. In reality I have many more Y columns than
simply y1,
y2, and y3, so I must design a loop. Notice that y2 is invariant
and thus it
will fail. In reality, some y columns will fail for much more subtle
reasons.
Simply screening my data to eliminate invariant columns will not
eliminate the
problem.
What I want to do is output a piece of the results from each run of
the loop to
a matrix. I want the to try each of my y columns, and not give up
and stop
running simply because a particular y column is bad. I want it to
give me NA
or something similar in my results matrix for the bad y columns, but
I want it
to keep going give me good data for the good y columns.
For instance:
results - matrix(nrow = 1, ncol = 3)
colnames(results) - c(y1, y2, y3)
for (i in 1:2) {
mod.poly3 - lrm(x[,i] ~ pol(x1, 3) + pol(x2, 3), data=x)
results[1,i] - anova(mod.poly3)[1,3]
}
If I run this code, it gives up when fitting y2 because the y2 is
bad. It
doesn't even try to fit y3. Here's what my console shows:
results
y1 y2 y3
[1,] 0.6976063 NA NA
As you can see, it gave up before fitting y3, which would have worked.
How do I force my code to keep going through the loop, despite the
rotten apples
it encounters along the way?
?try
http://cran.r-project.org/doc/FAQ/R-FAQ.html#How-can-I-capture-or-ignore-errors-in-a-long-simulation_003f
(Doesn't only apply to simulations.)
Exact code that gets the job done is what I am
interested in. I am a post-doc -- I am not taking any classes. I
promise this is
not a homework assignment!
--
David Winsemius, MD
West Hartford, CT
__
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] Xyplot or Tin-R problem?
On Jul 12, 2010, at 6:26 PM, YANG, Richard ChunHSi wrote:
I ran the following script from xyplot Examples using Tin-R on
Windows and saw no plot produced.
EE - equal.count(ethanol$E, number=9, overlap=1/4)
xyplot(NOx ~ C | EE, data=ethanol,
prepanel = function(x,y) prepanel.loess(x, y, span=1),
xlab=Compression Ratio, ylab=NOx (micrograms/J),
panel = function(x,y) {
panel.grid()(h = -1, v=2)
panel.xyplot(x,y)
panel.loess(x,y, span=1)
},
aspect = xy)
The Rgui showed
source(.trPaths[5])
Without any error msg. Did I miss anything? Please enlighten me.
I got the example to work fine but had no plotting with your version
and cannot see the difference in the code. I assigned them to t1 nd t2
and ...
all.equal(t1, t2)
[1] Component 5: target, current do not match when deparsed
[2] Component 29: target, current do not match when deparsed
Looking at str applied to both does not illuminate me. I have seen
problems on my Mac with examples copied from the help page and I
suspect there is some invisible character sitting in a copy-pasted
version that out mail-clients are not displaying. What happens if you
try:
examples(xyplot) #???
--
David Winsemius, MD
West Hartford, CT
__
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] Continuing on with a loop when there's a failure
On Jul 13, 2010, at 8:47 AM, Josh B wrote:
Thanks again, David.
...but, alas, I still can't get it work! Here's what I'm trying now:
for (i in 1:2) {
mod.poly3 - try(lrm(x[,i] ~ pol(x1, 3) + pol(x2, 3), data=x))
results[1,i] - anova(mod.poly3)[1,3]
}
You need to do some programming. You did not get an error from the lrm
but rather from the anova call because you tried to give the results
of the try function to anova without first checking to see if an error
had occurred.
--
David.
Here's what happens (from the console):
Error in fitter(X, Y, penalty.matrix = penalty.matrix, tol = tol,
weights = weights, :
NA/NaN/Inf in foreign function call (arg 1)
Error in UseMethod(anova) :
no applicable method for 'anova' applied to an object of class
try-error
...so I still can't make my results matrix. Could I ask you for some
specific code to make this work? I'm not that familiar with the
syntax for try or tryCatch, and the help files for them are pretty
bad, in my humble opinion.
I should clarify that I actually don't care about the failed runs
per se. I just want R to keep going in spite of them and give me my
results matrix.
From: David Winsemius dwinsem...@comcast.net
To: Josh B josh...@yahoo.com
Cc: R Help r-help@r-project.org
Sent: Mon, July 12, 2010 8:09:03 PM
Subject: Re: [R] Continuing on with a loop when there's a failure
On Jul 12, 2010, at 6:18 PM, Josh B wrote:
Hi R sages,
Here is my latest problem. Consider the following toy example:
x - read.table(textConnection(y1 y2 y3 x1 x2
indv.1 bagels donuts bagels 4 6
indv.2 donuts donuts donuts 5 1
indv.3 donuts donuts donuts 1 10
indv.4 donuts donuts donuts 10 9
indv.5 bagels donuts bagels 0 2
indv.6 bagels donuts bagels 2 9
indv.7 bagels donuts bagels 8 5
indv.8 bagels donuts bagels 4 1
indv.9 donuts donuts donuts 3 3
indv.10 bagels donuts bagels 5 9
indv.11 bagels donuts bagels 9 10
indv.12 bagels donuts bagels 3 1
indv.13 donuts donuts donuts 7 10
indv.14 bagels donuts bagels 2 10
indv.15 bagels donuts bagels 9 6), header = TRUE)
I want to fit a logistic regression of y1 on x1 and x2. Then I
want to run a
logistic regression of y2 on x1 and x2. Then I want to run a
logistic regression
of y3 on x1 and x2. In reality I have many more Y columns than
simply y1,
y2, and y3, so I must design a loop. Notice that y2 is
invariant and thus it
will fail. In reality, some y columns will fail for much more
subtle reasons.
Simply screening my data to eliminate invariant columns will not
eliminate the
problem.
What I want to do is output a piece of the results from each run
of the loop to
a matrix. I want the to try each of my y columns, and not give up
and stop
running simply because a particular y column is bad. I want it to
give me NA
or something similar in my results matrix for the bad y columns,
but I want it
to keep going give me good data for the good y columns.
For instance:
results - matrix(nrow = 1, ncol = 3)
colnames(results) - c(y1, y2, y3)
for (i in 1:2) {
mod.poly3 - lrm(x[,i] ~ pol(x1, 3) + pol(x2, 3), data=x)
results[1,i] - anova(mod.poly3)[1,3]
}
If I run this code, it gives up when fitting y2 because the y2 is
bad. It
doesn't even try to fit y3. Here's what my console shows:
results
y1 y2 y3
[1,] 0.6976063 NA NA
As you can see, it gave up before fitting y3, which would have
worked.
How do I force my code to keep going through the loop, despite the
rotten apples
it encounters along the way?
?try
http://cran.r-project.org/doc/FAQ/R-FAQ.html#How-can-I-capture-or-ignore-errors-in-a-long-simulation_003f
(Doesn't only apply to simulations.)
Exact code that gets the job done is what I am
interested in. I am a post-doc -- I am not taking any classes. I
promise this is
not a homework assignment!
--
David Winsemius, MD
West Hartford, CT
David Winsemius, MD
West Hartford, CT
__
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] Continuing on with a loop when there's a failure
On Jul 13, 2010, at 9:04 AM, David Winsemius wrote:
On Jul 13, 2010, at 8:47 AM, Josh B wrote:
Thanks again, David.
...but, alas, I still can't get it work!
(BTW, it did work.)
Here's what I'm trying now:
for (i in 1:2) {
mod.poly3 - try(lrm(x[,i] ~ pol(x1, 3) + pol(x2, 3), data=x))
results[1,i] - anova(mod.poly3)[1,3]
}
You need to do some programming.
(Or I suppose you could wrap both the lrm and the anova calls in try.)
You did not get an error from the lrm but rather from the anova call
because you tried to give the results of the try function to anova
without first checking to see if an error had occurred.
--
David.
Here's what happens (from the console):
Error in fitter(X, Y, penalty.matrix = penalty.matrix, tol = tol,
weights = weights, :
NA/NaN/Inf in foreign function call (arg 1)
Error in UseMethod(anova) :
no applicable method for 'anova' applied to an object of class
try-error
...so I still can't make my results matrix. Could I ask you for
some specific code to make this work? I'm not that familiar with
the syntax for try or tryCatch, and the help files for them are
pretty bad, in my humble opinion.
I should clarify that I actually don't care about the failed runs
per se. I just want R to keep going in spite of them and give me my
results matrix.
From: David Winsemius dwinsem...@comcast.net
To: Josh B josh...@yahoo.com
Cc: R Help r-help@r-project.org
Sent: Mon, July 12, 2010 8:09:03 PM
Subject: Re: [R] Continuing on with a loop when there's a failure
On Jul 12, 2010, at 6:18 PM, Josh B wrote:
Hi R sages,
Here is my latest problem. Consider the following toy example:
x - read.table(textConnection(y1 y2 y3 x1 x2
indv.1 bagels donuts bagels 4 6
indv.2 donuts donuts donuts 5 1
indv.3 donuts donuts donuts 1 10
indv.4 donuts donuts donuts 10 9
indv.5 bagels donuts bagels 0 2
indv.6 bagels donuts bagels 2 9
indv.7 bagels donuts bagels 8 5
indv.8 bagels donuts bagels 4 1
indv.9 donuts donuts donuts 3 3
indv.10 bagels donuts bagels 5 9
indv.11 bagels donuts bagels 9 10
indv.12 bagels donuts bagels 3 1
indv.13 donuts donuts donuts 7 10
indv.14 bagels donuts bagels 2 10
indv.15 bagels donuts bagels 9 6), header = TRUE)
I want to fit a logistic regression of y1 on x1 and x2. Then I
want to run a
logistic regression of y2 on x1 and x2. Then I want to run a
logistic regression
of y3 on x1 and x2. In reality I have many more Y columns than
simply y1,
y2, and y3, so I must design a loop. Notice that y2 is
invariant and thus it
will fail. In reality, some y columns will fail for much more
subtle reasons.
Simply screening my data to eliminate invariant columns will not
eliminate the
problem.
What I want to do is output a piece of the results from each run
of the loop to
a matrix. I want the to try each of my y columns, and not give up
and stop
running simply because a particular y column is bad. I want it to
give me NA
or something similar in my results matrix for the bad y columns,
but I want it
to keep going give me good data for the good y columns.
For instance:
results - matrix(nrow = 1, ncol = 3)
colnames(results) - c(y1, y2, y3)
for (i in 1:2) {
mod.poly3 - lrm(x[,i] ~ pol(x1, 3) + pol(x2, 3), data=x)
results[1,i] - anova(mod.poly3)[1,3]
}
If I run this code, it gives up when fitting y2 because the y2 is
bad. It
doesn't even try to fit y3. Here's what my console shows:
results
y1 y2 y3
[1,] 0.6976063 NA NA
As you can see, it gave up before fitting y3, which would have
worked.
How do I force my code to keep going through the loop, despite
the rotten apples
it encounters along the way?
?try
http://cran.r-project.org/doc/FAQ/R-FAQ.html#How-can-I-capture-or-ignore-errors-in-a-long-simulation_003f
(Doesn't only apply to simulations.)
Exact code that gets the job done is what I am
interested in. I am a post-doc -- I am not taking any classes. I
promise this is
not a homework assignment!
--
David Winsemius, MD
West Hartford, CT
David Winsemius, MD
West Hartford, CT
__
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, MD
West Hartford, CT
__
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] Continuing on with a loop when there's a failure
On Jul 13, 2010, at 9:24 AM, Josh B wrote:
In my opinion the try and tryCatch commands are written and
documented rather poorly. Thus I am not sure what to program exactly.
Didn't you see the silent parameter? Its seems to be documented
fairly clearly to me.
The testing of try at the console is not going to be very
illuminating, since it really only has value within a function that
you want to continue despite an error. try() WILL provide that
facility but _you_ need to decide what you do with the information it
returns, which in the case of its use with the default silent=FALSE is
just the error message itself.
For instance, I could query mod.poly3 and use an if/then statement
to proceed,
So why didn't you? A good result would be signaled by: lrm %in
class(mod.poly3)
--
David.
but querying mod.poly3 is weird. For instance, here's the output
when it fails:
mod.poly3 - try(lrm(x[,2] ~ pol(x1, 3) + pol(x2, 3), data=x))
Error in fitter(X, Y, penalty.matrix = penalty.matrix, tol = tol,
weights = weights, :
NA/NaN/Inf in foreign function call (arg 1)
mod.poly3
[1] Error in fitter(X, Y, penalty.matrix = penalty.matrix, tol =
tol, weights = weights, : \n NA/NaN/Inf in foreign function call
(arg 1)\n
attr(,class)
[1] try-error
...and here's the output when it succeeds:
mod.poly3 - try(lrm(x[,1] ~ pol(x1, 3) + pol(x2, 3), data=x))
mod.poly3
Logistic Regression Model
lrm(formula = x[, 1] ~ pol(x1, 3) + pol(x2, 3), data = x)
Frequencies of Responses
bagels donuts
10 5
Obs Max Deriv Model L.R. d.f. P C
15 4e-04 3.37 6 0.7616 0.76
Dxy Gamma Tau-a R2 Brier g
0.52 0.52 0.248 0.279 0.183 1.411
gr gp
4.1 0.261
Coef S.E.Wald Z P
Intercept -5.68583 5.23295 -1.09 0.2772
x1 1.87020 2.14635 0.87 0.3836
x1^2 -0.42494 0.48286 -0.88 0.3788
x1^3 0.02845 0.03120 0.91 0.3618
x2 3.49560 3.54796 0.99 0.3245
x2^2 -0.94888 0.82067 -1.16 0.2476
x2^3 0.06362 0.05098 1.25 0.2121
...so what exactly would I query to design my if/then statement?
From: David Winsemius dwinsem...@comcast.net
To: David Winsemius dwinsem...@comcast.net
Cc: Josh B josh...@yahoo.com; R Help r-help@r-project.org
Sent: Tue, July 13, 2010 9:09:04 AM
Subject: Re: [R] Continuing on with a loop when there's a failure
On Jul 13, 2010, at 9:04 AM, David Winsemius wrote:
On Jul 13, 2010, at 8:47 AM, Josh B wrote:
Thanks again, David.
...but, alas, I still can't get it work!
(BTW, it did work.)
Here's what I'm trying now:
for (i in 1:2) {
mod.poly3 - try(lrm(x[,i] ~ pol(x1, 3) + pol(x2, 3), data=x))
results[1,i] - anova(mod.poly3)[1,3]
}
You need to do some programming.
(Or I suppose you could wrap both the lrm and the anova calls in try.)
You did not get an error from the lrm but rather from the anova
call because you tried to give the results of the try function to
anova without first checking to see if an error had occurred.
--David.
Here's what happens (from the console):
Error in fitter(X, Y, penalty.matrix = penalty.matrix, tol = tol,
weights = weights, :
NA/NaN/Inf in foreign function call (arg 1)
Error in UseMethod(anova) :
no applicable method for 'anova' applied to an object of class
try-error
...so I still can't make my results matrix. Could I ask you for
some specific code to make this work? I'm not that familiar with the
syntax for try or tryCatch, and the help files for them are pretty
bad, in my humble opinion.
I should clarify that I actually don't care about the failed runs
per se. I just want R to keep going in spite of them and give me my
results matrix.
From: David Winsemius dwinsem...@comcast.net
To: Josh B josh...@yahoo.com
Cc: R Help r-help@r-project.org
Sent: Mon, July 12, 2010 8:09:03 PM
Subject: Re: [R] Continuing on with a loop when there's a failure
On Jul 12, 2010, at 6:18 PM, Josh B wrote:
Hi R sages,
Here is my latest problem. Consider the following toy example:
x - read.table(textConnection(y1 y2 y3 x1 x2
indv.1 bagels donuts bagels 4 6
indv.2 donuts donuts donuts 5 1
indv.3 donuts donuts donuts 1 10
indv.4 donuts donuts donuts 10 9
indv.5 bagels donuts bagels 0 2
indv.6 bagels donuts bagels 2 9
indv.7 bagels donuts bagels 8 5
indv.8 bagels donuts bagels 4 1
indv.9 donuts donuts donuts 3 3
indv.10 bagels donuts bagels 5 9
indv.11 bagels donuts bagels 9 10
indv.12 bagels donuts bagels 3 1
indv.13 donuts donuts donuts 7 10
indv.14 bagels donuts bagels 2 10
indv.15 bagels donuts bagels 9 6), header = TRUE)
I want to fit a logistic regression of y1 on x1 and x2. Then I
want to run a
logistic regression of y2 on x1 and x2. Then I want to run a
logistic regression
of y3 on x1 and x2. In reality I have many more Y
Re: [R] Continuing on with a loop when there's a failure
On Jul 13, 2010, at 10:26 AM, David Winsemius wrote:
On Jul 13, 2010, at 9:24 AM, Josh B wrote:
In my opinion the try and tryCatch commands are written and
documented rather poorly. Thus I am not sure what to program exactly.
Didn't you see the silent parameter? Its seems to be documented
fairly clearly to me.
The testing of try at the console is not going to be very
illuminating, since it really only has value within a function that
you want to continue despite an error. try() WILL provide that
facility but _you_ need to decide what you do with the information
it returns, which in the case of its use with the default
silent=FALSE is just the error message itself.
For instance, I could query mod.poly3 and use an if/then statement
to proceed,
So why didn't you? A good result would be signaled by:
rather: lrm %in% class(mod.poly3)
--
David.
but querying mod.poly3 is weird. For instance, here's the output
when it fails:
mod.poly3 - try(lrm(x[,2] ~ pol(x1, 3) + pol(x2, 3), data=x))
Error in fitter(X, Y, penalty.matrix = penalty.matrix, tol = tol,
weights = weights, :
NA/NaN/Inf in foreign function call (arg 1)
mod.poly3
[1] Error in fitter(X, Y, penalty.matrix = penalty.matrix, tol =
tol, weights = weights, : \n NA/NaN/Inf in foreign function call
(arg 1)\n
attr(,class)
[1] try-error
...and here's the output when it succeeds:
mod.poly3 - try(lrm(x[,1] ~ pol(x1, 3) + pol(x2, 3), data=x))
mod.poly3
Logistic Regression Model
lrm(formula = x[, 1] ~ pol(x1, 3) + pol(x2, 3), data = x)
Frequencies of Responses
bagels donuts
10 5
Obs Max Deriv Model L.R. d.f. P C
15 4e-04 3.37 6 0.7616 0.76
Dxy Gamma Tau-a R2 Brier g
0.52 0.52 0.248 0.279 0.183 1.411
gr gp
4.1 0.261
Coef S.E.Wald Z P
Intercept -5.68583 5.23295 -1.09 0.2772
x1 1.87020 2.14635 0.87 0.3836
x1^2 -0.42494 0.48286 -0.88 0.3788
x1^3 0.02845 0.03120 0.91 0.3618
x2 3.49560 3.54796 0.99 0.3245
x2^2 -0.94888 0.82067 -1.16 0.2476
x2^3 0.06362 0.05098 1.25 0.2121
...so what exactly would I query to design my if/then statement?
From: David Winsemius dwinsem...@comcast.net
To: David Winsemius dwinsem...@comcast.net
Cc: Josh B josh...@yahoo.com; R Help r-help@r-project.org
Sent: Tue, July 13, 2010 9:09:04 AM
Subject: Re: [R] Continuing on with a loop when there's a failure
On Jul 13, 2010, at 9:04 AM, David Winsemius wrote:
On Jul 13, 2010, at 8:47 AM, Josh B wrote:
Thanks again, David.
...but, alas, I still can't get it work!
(BTW, it did work.)
Here's what I'm trying now:
for (i in 1:2) {
mod.poly3 - try(lrm(x[,i] ~ pol(x1, 3) + pol(x2, 3), data=x))
results[1,i] - anova(mod.poly3)[1,3]
}
You need to do some programming.
(Or I suppose you could wrap both the lrm and the anova calls in
try.)
You did not get an error from the lrm but rather from the anova
call because you tried to give the results of the try function to
anova without first checking to see if an error had occurred.
--David.
Here's what happens (from the console):
Error in fitter(X, Y, penalty.matrix = penalty.matrix, tol =
tol, weights = weights, :
NA/NaN/Inf in foreign function call (arg 1)
Error in UseMethod(anova) :
no applicable method for 'anova' applied to an object of class
try-error
...so I still can't make my results matrix. Could I ask you for
some specific code to make this work? I'm not that familiar with
the syntax for try or tryCatch, and the help files for them are
pretty bad, in my humble opinion.
I should clarify that I actually don't care about the failed
runs per se. I just want R to keep going in spite of them and give
me my results matrix.
From: David Winsemius dwinsem...@comcast.net
To: Josh B josh...@yahoo.com
Cc: R Help r-help@r-project.org
Sent: Mon, July 12, 2010 8:09:03 PM
Subject: Re: [R] Continuing on with a loop when there's a failure
On Jul 12, 2010, at 6:18 PM, Josh B wrote:
Hi R sages,
Here is my latest problem. Consider the following toy example:
x - read.table(textConnection(y1 y2 y3 x1 x2
indv.1 bagels donuts bagels 4 6
indv.2 donuts donuts donuts 5 1
indv.3 donuts donuts donuts 1 10
indv.4 donuts donuts donuts 10 9
indv.5 bagels donuts bagels 0 2
indv.6 bagels donuts bagels 2 9
indv.7 bagels donuts bagels 8 5
indv.8 bagels donuts bagels 4 1
indv.9 donuts donuts donuts 3 3
indv.10 bagels donuts bagels 5 9
indv.11 bagels donuts bagels 9 10
indv.12 bagels donuts bagels 3 1
indv.13 donuts donuts donuts 7 10
indv.14 bagels donuts bagels 2 10
indv.15 bagels donuts bagels 9 6), header = TRUE)
I want to fit a logistic regression of y1 on x1 and x2. Then I
want to run a
logistic regression of y2 on x1 and x2. Then I want to run a
logistic
Re: [R] how to extract information from anova results
On Jul 13, 2010, at 10:35 AM, Luis Borda de Agua wrote: Hi, I have used the instruction aov in the following manner: res - aov(qwe ~ asd) when I typed res I get: _ Call: aov(formula = qwe ~ asd) Terms: asd Residuals Sum of Squares 0.0708704 0.5255957 Deg. of Freedom 1 8 Residual standard error: 0.2563191 Estimated effects may be unbalanced _ I need to access the value of the Sum of Squares (i.e. I want another variable to be equal to it, e.g myvar - Sum.of.Squares) . I tried names(res) to see which values are accessible, but I couldn't find the Sum of Squares. I had a similar problem when I tried to access the p.value which can be readily SEEN using summary(res). In general, is there an easy way to access the values generated by an R function? When you typed res, the interpreter determined that it was of type aov and dispatched it to the print method for objects of that class. The list of print methods is accessed with: methods(print) and it's a long list. print.aov is asterisked so you either need to look at the function with: getAnywhere(print.aov) or perhaps more directly assign summary(res) to an object and access its SS values. -- David. Thank you, LBA __ 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, MD West Hartford, CT __ 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] Continuing on with a loop when there's a failure
On Jul 18, 2010, at 4:25 AM, Josh B wrote:
Hello Peter,
I tried your suggestion, but I was still not able to get it to work.
Would you
mind looking at my code again? Here's what I'm trying:
x - read.table(textConnection(y1 y2 y3 x1 x2
indv.1 bagels donuts bagels 4 6
indv.2 donuts donuts donuts 5 1
indv.3 donuts donuts donuts 1 10
indv.4 donuts donuts donuts 10 9
indv.5 bagels donuts bagels 0 2
indv.6 bagels donuts bagels 2 9
indv.7 bagels donuts bagels 8 5
indv.8 bagels donuts bagels 4 1
indv.9 donuts donuts donuts 3 3
indv.10 bagels donuts bagels 5 9
indv.11 bagels donuts bagels 9 10
indv.12 bagels donuts bagels 3 1
indv.13 donuts donuts donuts 7 10
indv.14 bagels donuts bagels 2 10
indv.15 bagels donuts bagels 9 6), header = TRUE)
closeAllConnections()
results - matrix(nrow = 1, ncol = 3)
colnames(results) - c(y1, y2, y3)
require(rms) # or Design
for (i in 1:2) {
mod.poly3 - try(lrm(x[,i] ~ pol(x1, 3) + pol(x2, 3), data=x),
silent=TRUE)
if(class(mod.poly3)[1] != 'try-error')
{results[1,i] - anova(mod.poly3)[1,3]}
}
results
y1 y2 y3
[1,] 0.6976063 NA NA
...and here's the output:
results
y1 y2 y3
[1,] NA NA NA
The results matrix is empty!
From: Peter Konings peter.l.e.koni...@gmail.com
Sent: Tue, July 13, 2010 5:45:17 PM
Subject: Re: [R] Continuing on with a loop when there's a failure
Hi Josh,
Test the class of the resulting object. If it is 'try-error' fill
your result
with NA or do some other error handling.
result - try(somemodel)
if(class(result) == 'try-error')
{
# some error handling
} else {
# whatever happens if the result is ok
}
HTH
Peter.
In my opinion the try and tryCatch commands are written and
documented rather
poorly. Thus I am not sure what to program exactly.
For instance, I could query mod.poly3 and use an if/then statement
to proceed,
but querying mod.poly3 is weird. For instance, here's the output
when it fails:
mod.poly3 - try(lrm(x[,2] ~ pol(x1, 3) + pol(x2, 3), data=x))
Error in fitter(X, Y, penalty.matrix = penalty.matrix, tol = tol,
weights =
weights, :
NA/NaN/Inf in foreign function call (arg 1)
mod.poly3
[1] Error in fitter(X, Y, penalty.matrix = penalty.matrix, tol =
tol, weights
=
weights, : \n NA/NaN/Inf in foreign function call (arg 1)\n
attr(,class)
[1] try-error
...and here's the output when it succeeds:
mod.poly3 - try(lrm(x[,1] ~ pol(x1, 3) + pol(x2, 3), data=x))
mod.poly3
Logistic Regression Model
lrm(formula = x[, 1] ~ pol(x1, 3) + pol(x2, 3), data = x)
Frequencies of Responses
bagels donuts
10 5
Obs Max Deriv Model L.R. d.f. P C
15 4e-04 3.37 6 0.7616 0.76
Dxy Gamma Tau-a R2 Brier g
0.52 0.52 0.248 0.279 0.183 1.411
gr gp
4.1 0.261
Coef S.E.Wald Z P
Intercept -5.68583 5.23295 -1.09 0.2772
x1 1.87020 2.14635 0.87 0.3836
x1^2 -0.42494 0.48286 -0.88 0.3788
x1^3 0.02845 0.03120 0.91 0.3618
x2 3.49560 3.54796 0.99 0.3245
x2^2 -0.94888 0.82067 -1.16 0.2476
x2^3 0.06362 0.05098 1.25 0.2121
...so what exactly would I query to design my if/then statement?
From: David Winsemius dwinsem...@comcast.net
To: David Winsemius dwinsem...@comcast.net
Sent: Tue, July 13, 2010 9:09:04 AM
Subject: Re: [R] Continuing on with a loop when there's a failure
On Jul 13, 2010, at 9:04 AM, David Winsemius wrote:
On Jul 13, 2010, at 8:47 AM, Josh B wrote:
Thanks again, David.
[[elided Yahoo spam]]
(BTW, it did work.)
Here's what I'm trying now:
for (i in 1:2) {
mod.poly3 - try(lrm(x[,i] ~ pol(x1, 3) + pol(x2, 3), data=x))
results[1,i] - anova(mod.poly3)[1,3]
}
You need to do some programming.
(Or I suppose you could wrap both the lrm and the anova calls in
try.)
You did not get an error from the lrm but rather from the anova
call because
you tried to give the results of the try function to anova without
first
checking to see if an error had occurred.
--David.
Here's what happens (from the console):
Error in fitter(X, Y, penalty.matrix = penalty.matrix, tol = tol,
weights =
weights, :
NA/NaN/Inf in foreign function call (arg 1)
Error in UseMethod(anova) :
no applicable method for 'anova' applied to an object of class
try-error
...so I still can't make my results matrix. Could I ask you for
some
specific
code to make this work? I'm not that familiar with the syntax for
try or
tryCatch, and the help files for them are pretty bad, in my
humble opinion.
I should clarify that I actually don't care about the failed runs
per se. I
just want R to keep going in spite of them and give me my results
matrix.
From: David Winsemius dwinsem...@comcast.net
Cc: R Help r-help@r-project.org
Sent: Mon, July 12, 2010 8:09:03 PM
Subject: Re: [R] Continuing
Re: [R] sort file names in numerical order
Another option: require(gtools) ?mixedsort mixedsort(fileNames) [1] A1 A2 A10 B1 B2 B10 -- David On Jul 18, 2010, at 5:16 AM, Duncan Mackay wrote: Hi Yes it is possible- one way is: fileNames[order(sprintf(%02s, sub([[:upper:]],, fileNames)))] [1] A1 B1 A2 B2 A10 B10 Regards Duncan Duncan Mackay Department of Agronomy and Soil Science University of New England ARMIDALE NSW 2351 Email home: mac...@northnet.com.au At 06:47 18/07/2010, you wrote: I get some file names by list.files(). These names are in alphabetical order. I want to change it to logical numeric order. Example: fileNames - c(A10, A1, A2, B1, B2, B10) sort(fileNames) [1] A1 A10 A2 B1 B10 B2 I want to have: A1 A2 A10 B1 B2 B10 Is this possible? Greetings. I see that you've gotten an elegant solution to your problem. I've appended a poor-man's solution, which I generated more for my own edification than for yours. -- Mike ## modified file names, changed to exhibit sorting on letters fileNames - c(A10, B10, A1, A2, B1, B2); fileNames ## use regular expressions to pick off letters, numbers fnLet - gsub((^[^0-9]+).*$, \\1, fileNames); fnLet fnNum - gsub(^[^0-9]+(.*)$, \\1, fileNames); fnNum ## need to force numeric sorting of the numbers fnNum - as.numeric(fnNum) ## find the order of the numbers, for later use in subsetting fnNumOrd - order(fnNum); fnNumOrd ## pack all the relevant information into a data frame, then select from that fnPieces - data.frame(fnLet, fnNum, fileNames, stringsAsFactors=FALSE); fnPieces ## partial sort by numbers only (gives new df) dfPartialSort - fnPieces[fnNumOrd, ]; dfPartialSort ## find the order of the names, then select from new df with that fnLetOrd - order(dfPartialSort[, 1]); fnLetOrd dfFullSort - dfPartialSort[fnLetOrd, ]; dfFullSort ## the file names have gone along for the ride, so pick them off now sortedFileNames - dfFullSort$fileNames; sortedFileNames __ 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. David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] loop troubles
On Jul 17, 2010, at 9:09 PM, Joe P King wrote:
Hi all, I appreciate the help this list has given me before. I have a
question which has been perplexing me. I have been working on doing a
Bayesian calculating inserting studies sequentially after using a
non-informative prior to get a meta-analysis type result. I created a
function using three iterations of this, my code is below. I insert
prior
mean and precision (I add precision manually because I want it to be
zero
but if I include zero as prior variance I get an error cant divide
by zero.
Now my question is instead of having the three iterations below, is
there
anyway to create a loop, the only problem is I want to have elements
from
the previous posterior to be the new prior and now I cant figure out
how to
do the code below in a loop.
Perhaps if you set the problem up with vectors, it would become more
obvious how to do it with indexing or loops:
n - c(5, 10, 15)
ybar - c(12,12,14) # and so on...
--
David
The data below is dummy data, I used a starting
mu of 1, and starting precision of 0.
bayes.analysis.treat-function(mu0,p0){
n1 = 5
n2 = 10
n3 = 15
ybar1 = 12
ybar2 = 13
ybar3 = 14
sd1 = 2
sd2 = 3
sd3 = 4
#posterior 1
var1 = sd1^2 #sample variance
p1 = n1/var1 #sample precision
p1n = p0+p1
mu1 = ((p0)/(p1n)*mu0)+((p1)/(p1n))*ybar1
sigma1 = 1/sqrt(p1n)
#posterior 2
var2 = sd2^2 #sample variance
p2 = n2/var2 #sample precision
p2n = p1n+p2
mu2 = ((p1n)/(p2n)*mu1)+((p2)/(p2n))*ybar2
sigma2 = 1/sqrt(p2n)
#posterior 3
var3 = sd3^2 #sample variance
p3 = n3/var3 #sample precision
p3n = p2n+p3
mu3 = ((p2n)/(p3n)*mu2)+((p3)/(p3n))*ybar3
sigma3 = 1/sqrt(p3n)
posterior-cbind(rbind(mu1,mu2,mu3),rbind(sigma1,sigma2,sigma3))
rownames(posterior)-c(Posterior 1, Posterior 2, Posterior 3)
colnames(posterior)-c(Mu, SD)
return(posterior)}
---
Joe King, M.A.
Ph.D. Student
University of Washington - Seattle
206-913-2912
mailto:j...@joepking.com j...@joepking.com
---
Never throughout history has a man who lived a life of ease left a
name
worth remembering. --Theodore Roosevelt
[[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, MD
West Hartford, CT
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] loop troubles
On Jul 18, 2010, at 9:12 AM, Joe P King wrote:
I tried that, this is what I tried, but I only get it to do one
iteration
and then it wont cycle back. I am not sure how to tell it how to
look at the
previous precision to add to the current new sample precision to get
the new
precision.
And I'm not exactly sure what you are trying to do since I found your
problem statement insufficiently precise, but if you have a
vector, p and an index in a loop is i, then the previous p is
just p[i-1]. Sometimes you can get constructs like the following to
accomplish an indexed assignment:
p[2:length(n)] -a*p[1:(length(n)-1)]
This is my first try at a loop so I am not sure if I am doing
anything right, sorry about the elementary question.
n=c(5,10,15)
ybar=c(12,13,14)
sd=c(2,3,4)
mutotals-matrix(nrow=length(n), ncol=1)
mu=1;p0=0
for (i in 1:length(n)){
var = sd[i]^2 #sample variance
p = n[i]/var #sample precision
At this point you seem to be using p as a constant
p0[i]= i
Comments might help to figure out why you are doing each of htese steps.
pn = p0[i]+p[i]
Wouldn't pn need to be indexed as well? And now you seem to be using
p as a vector but haven't extended it with c() or some other
function, so on the second iteration it should fail.
mu.out = ((p0[i])/(pn)*mu[i])+((p)/(pn))*ybar[i]
sigma[i] = 1/sqrt(pn[i])
mutotals[i,]-mu.out[i]
}
I get:
Error in sigma[i] = 1/sqrt(pn[i]) : object 'sigma' not found
You are not reporting the error messages this throws (which should
have made it apparent that you needed to do some further setup of the
data structures you are using, not just with p, pn, but also
sigma. )
Joe King
206-913-2912
j...@joepking.com
Never throughout history has a man who lived a life of ease left a
name
worth remembering. --Theodore Roosevelt
-Original Message-
From: David Winsemius [mailto:dwinsem...@comcast.net]
Sent: Sunday, July 18, 2010 5:36 AM
To: Joe P King
Cc: r-help@r-project.org
Subject: Re: [R] loop troubles
On Jul 17, 2010, at 9:09 PM, Joe P King wrote:
Hi all, I appreciate the help this list has given me before. I have a
question which has been perplexing me. I have been working on doing a
Bayesian calculating inserting studies sequentially after using a
non-informative prior to get a meta-analysis type result. I created a
function using three iterations of this, my code is below. I insert
prior
mean and precision (I add precision manually because I want it to be
zero
but if I include zero as prior variance I get an error cant divide
by zero.
Now my question is instead of having the three iterations below, is
there
anyway to create a loop, the only problem is I want to have elements
from
the previous posterior to be the new prior and now I cant figure out
how to
do the code below in a loop.
Perhaps if you set the problem up with vectors, it would become more
obvious how to do it with indexing or loops:
n - c(5, 10, 15)
ybar - c(12,12,14) # and so on...
--
David
The data below is dummy data, I used a starting
mu of 1, and starting precision of 0.
bayes.analysis.treat-function(mu0,p0){
n1 = 5
n2 = 10
n3 = 15
ybar1 = 12
ybar2 = 13
ybar3 = 14
sd1 = 2
sd2 = 3
sd3 = 4
#posterior 1
var1 = sd1^2 #sample variance
p1 = n1/var1 #sample precision
p1n = p0+p1
mu1 = ((p0)/(p1n)*mu0)+((p1)/(p1n))*ybar1
sigma1 = 1/sqrt(p1n)
#posterior 2
var2 = sd2^2 #sample variance
p2 = n2/var2 #sample precision
p2n = p1n+p2
mu2 = ((p1n)/(p2n)*mu1)+((p2)/(p2n))*ybar2
sigma2 = 1/sqrt(p2n)
#posterior 3
var3 = sd3^2 #sample variance
p3 = n3/var3 #sample precision
p3n = p2n+p3
mu3 = ((p2n)/(p3n)*mu2)+((p3)/(p3n))*ybar3
sigma3 = 1/sqrt(p3n)
posterior-cbind(rbind(mu1,mu2,mu3),rbind(sigma1,sigma2,sigma3))
rownames(posterior)-c(Posterior 1, Posterior 2, Posterior 3)
colnames(posterior)-c(Mu, SD)
return(posterior)}
---
Joe King, M.A.
Ph.D. Student
University of Washington - Seattle
206-913-2912
mailto:j...@joepking.com j...@joepking.com
---
Never throughout history has a man who lived a life of ease left a
name
worth remembering. --Theodore Roosevelt
[[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, MD
West Hartford, CT
__
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, MD
West Hartford, CT
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman
Re: [R] Identify points (was Plot error)
On Jul 18, 2010, at 9:43 AM, James Platt wrote: This is exactly what I want as I will have several thousand data points on the final graph i make, so the scroll over option is ideal. I've read the TeachingDemos pdf, I'm working on a Mac so Cannot use HWidentify. Not true. I have installed and loaded the TeachingDemos package but im having trouble loading the tkrplot package after installing it won't load. library(tcltk) In the Mac GUI Package Installer, I cannot even find a package by that name. However, then running the HWidentify example after require(Teaching Demos), tcltk does get loaded as well as X11, making me think you need to provide further information such as that requested in the Posting Guide. You may also want to review the material in the R for Mac OS X FAQ. -- David. sessionInfo() R version 2.11.1 Patched (2010-06-14 r52281) x86_64-apple-darwin9.8.0 locale: [1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] tcltk splines stats graphics grDevices utils [7] datasets methods base other attached packages: [1] tkrplot_0.0-18TeachingDemos_2.6 gtools_2.6.2 [4] Design_2.3-0 Hmisc_3.8-1 survival_2.35-8 [7] lattice_0.18-8 loaded via a namespace (and not attached): [1] cluster_1.12.3 grid_2.11.1tools_2.11.1 Loading Tcl/Tk interface ... Error : .onLoad failed in loadNamespace() for 'tcltk', details: call: dyn.load(file, DLLpath = DLLpath, ...) error: unable to load shared library '/Library/Frameworks/ R.framework/Resources/library/tcltk/libs/i386/tcltk.so': dlopen(/Library/Frameworks/R.framework/Resources/library/tcltk/libs/ i386/tcltk.so, 10): Library not loaded: /usr/local/lib/libtcl8.5.dylib Referenced from: /Library/Frameworks/R.framework/Resources/library/ tcltk/libs/i386/tcltk.so Reason: image not found Error: package/namespace load failed for 'tcltk' On 17 Jul 2010, at 18:12, Peter Ehlers wrote: On 2010-07-17 9:50, James Platt wrote: The other question I have: Is there any way to link the data point on the graph to the name of a row i.e in my table: name value_1 value_2 bill 14 ben 2 2 jane 3 1 I click on the data point at 2,2 and it would read out ben Check out HWidentify or HTKidentify in pkg:TeachingDemos. -Peter Ehlers David Winsemius, MD West Hartford, CT __ 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] NA preserved in logical call - I don't understand this behavior because NA is not equal to 0
On Jul 18, 2010, at 12:02 PM, stephen sefick wrote:
I am confused by the behavior of the below piece of code. The NAs are
making it past the logical call ==0. I am sure that I am missing
something. I just don't understand this behavior. Thanks for your
help in advance.
code###
left - (structure(list(measurment_num = c(2, 2.2, 2.4, 2.6, 2.8,
2.82,
3, NA, NA, NA), bankfull_depths_m = c(1.29, 1.28, 1.23, 0.18,
-0.05, 0, -0.09, NA, NA, NA)), .Names = c(measurment_num,
bankfull_depths_m
), row.names = c(10L, 11L, 12L, 13L, 14L, 20L, 15L, 16L, 17L,
18L), class = data.frame))
if(sum(left[,bankfull_depths_m]==0, na.rm=TRUE) == 1){
left_min - left[left[,bankfull_depths_m]==0, measurment_num]
}
left
Why aren't you looking at left_min? That is what the code is supposed
to be changing. I do not get an error and left_min[1] (now) =='s 2.82.
Regarding why there are the 3 NA's, ... you need to look at the
documentation for Extract in the section very appropriately entitled:
NAs in indexing
I have been bitten by that behavior more times than I care to admit
and I am not the only person that thinks that aspect of the R language
is badly designed. See Aniko Szabo's comments:
http://radfordneal.wordpress.com/2008/09/21/design-flaws-in-r-3-%E2%80%94-zero-subscripts/
The problem is in data.frame[ and any NA in a logical vector will
return a row of NA's. This can be avoid by wrapping which() around the
logical vector which seems entirely wasteful or using subset().
David Winsemius, MD
West Hartford, CT
__
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] Continuing on with a loop when there's a failure
On Jul 18, 2010, at 12:28 PM, Josh B wrote:
Thanks very much again David for your helpful answers. However, the
code STILL does not appear to be working properly!
Even though the third time through the loop *should* work, it
appears that R has given up after the second time through the loop.
What I mean is: although y2 causes the lrm function to fail, y3 is a
kosher variable. If the loop continues on, it should give data for
the model with y3.
Right, but your loop index only went to 2!!! Your design, not mine.
--
David.
But if you look at the matrix called results, it returns NA for
the third spot corresponding to the model of y3:
results
y1 y2 y3
[1,] 0.6976063 NA NA
If you run y3 in isolation, rather than through the loop, you can
see that it should work and contribute data to the matrix called
results:
mod.poly3 - lrm(x[,3] ~ pol(x1, 3) + pol(x2, 3), data=x)
anova(mod.poly3)[1,3]
[1] 0.6976063
Any ideas?
From: David Winsemius dwinsem...@comcast.net
To: Josh B josh...@yahoo.com
Cc: Peter Konings peter.l.e.koni...@gmail.com; R Help r-help@r-project.org
Sent: Sun, July 18, 2010 3:33:07 PM
Subject: Re: [R] Continuing on with a loop when there's a failure
On Jul 18, 2010, at 4:25 AM, Josh B wrote:
Hello Peter,
I tried your suggestion, but I was still not able to get it to
work. Would you
mind looking at my code again? Here's what I'm trying:
x - read.table(textConnection(y1 y2 y3 x1 x2
indv.1 bagels donuts bagels 4 6
indv.2 donuts donuts donuts 5 1
indv.3 donuts donuts donuts 1 10
indv.4 donuts donuts donuts 10 9
indv.5 bagels donuts bagels 0 2
indv.6 bagels donuts bagels 2 9
indv.7 bagels donuts bagels 8 5
indv.8 bagels donuts bagels 4 1
indv.9 donuts donuts donuts 3 3
indv.10 bagels donuts bagels 5 9
indv.11 bagels donuts bagels 9 10
indv.12 bagels donuts bagels 3 1
indv.13 donuts donuts donuts 7 10
indv.14 bagels donuts bagels 2 10
indv.15 bagels donuts bagels 9 6), header = TRUE)
closeAllConnections()
results - matrix(nrow = 1, ncol = 3)
colnames(results) - c(y1, y2, y3)
require(rms) # or Design
for (i in 1:2) {
mod.poly3 - try(lrm(x[,i] ~ pol(x1, 3) + pol(x2, 3), data=x),
silent=TRUE)
if(class(mod.poly3)[1] != 'try-error')
{results[1,i] - anova(mod.poly3)[1,3]}
}
results
y1 y2 y3
[1,] 0.6976063 NA NA
...and here's the output:
results
y1 y2 y3
[1,] NA NA NA
The results matrix is empty!
From: Peter Konings peter.l.e.koni...@gmail.com
Sent: Tue, July 13, 2010 5:45:17 PM
Subject: Re: [R] Continuing on with a loop when there's a failure
Hi Josh,
Test the class of the resulting object. If it is 'try-error' fill
your result
with NA or do some other error handling.
result - try(somemodel)
if(class(result) == 'try-error')
{
# some error handling
} else {
# whatever happens if the result is ok
}
HTH
Peter.
In my opinion the try and tryCatch commands are written and
documented rather
poorly. Thus I am not sure what to program exactly.
For instance, I could query mod.poly3 and use an if/then
statement to proceed,
but querying mod.poly3 is weird. For instance, here's the output
when it fails:
mod.poly3 - try(lrm(x[,2] ~ pol(x1, 3) + pol(x2, 3), data=x))
Error in fitter(X, Y, penalty.matrix = penalty.matrix, tol = tol,
weights =
weights, :
NA/NaN/Inf in foreign function call (arg 1)
mod.poly3
[1] Error in fitter(X, Y, penalty.matrix = penalty.matrix, tol =
tol, weights
=
weights, : \n NA/NaN/Inf in foreign function call (arg 1)\n
attr(,class)
[1] try-error
...and here's the output when it succeeds:
mod.poly3 - try(lrm(x[,1] ~ pol(x1, 3) + pol(x2, 3), data=x))
mod.poly3
Logistic Regression Model
lrm(formula = x[, 1] ~ pol(x1, 3) + pol(x2, 3), data = x)
Frequencies of Responses
bagels donuts
10 5
Obs Max Deriv Model L.R. d.f. P C
15 4e-04 3.37 60.7616 0.76
Dxy Gamma Tau-aR2 Brier g
0.52 0.52 0.248 0.279 0.183 1.411
grgp
4.1 0.261
CoefS.E.Wald Z P
Intercept -5.68583 5.23295 -1.09 0.2772
x11.87020 2.14635 0.87 0.3836
x1^2 -0.42494 0.48286 -0.88 0.3788
x1^3 0.02845 0.03120 0.91 0.3618
x23.49560 3.54796 0.99 0.3245
x2^2 -0.94888 0.82067 -1.16 0.2476
x2^3 0.06362 0.05098 1.25 0.2121
...so what exactly would I query to design my if/then statement?
From: David Winsemius dwinsem...@comcast.net
To: David Winsemius dwinsem...@comcast.net
Sent: Tue, July 13, 2010 9:09:04 AM
Subject: Re: [R] Continuing on with a loop when there's a failure
On Jul 13, 2010, at 9:04 AM, David Winsemius wrote:
On Jul 13, 2010, at 8:47 AM, Josh B wrote:
Thanks again, David.
[[elided Yahoo spam]]
(BTW, it did work.)
Here's what I'm
Re: [R] NA preserved in logical call - I don't understand this behavior because NA is not equal to 0
On Jul 18, 2010, at 2:52 PM, Hadley Wickham wrote: The problem is in data.frame[ and any NA in a logical vector will return a row of NA's. This can be avoid by wrapping which() around the logical vector which seems entirely wasteful or using subset(). The basic philosophy that causes this behaviour is sensible in my opinion: missing values must always be handled explicitly. Sure it's a bit of a hassle, but it being explicit prevents you from making major errors in an analysis. How does this apply to the treatment of logical indexing in the function data.frame[ ? The returned value is not expected to have the same number of rows as the source object and the profusion of duplicate NA rows serves what purpose? Wouldn't it be more consistent to have the constructs: dfrm[which(logical_vector), ] dfrm[logical_vector , ] return the same object? In my opinion, the flaw is that R doesn't apply this philosophy entirely consistently: I think factor and table should default to exclude = NULL. Hadley -- Assistant Professor / Dobelman Family Junior Chair Department of Statistics / Rice University http://had.co.nz/ David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] loop troubles
On Jul 18, 2010, at 4:21 PM, Joe P King wrote:
This is the latest code I have been trying, but it when I try to
turn the
vectors of SD into a vector of variances, it turns into a scalar.
n=c(NULL,13,89,50)
ybar=c(NULL,1.58,1.26,0.80)
sd=c(NULL,2.19,4.18,5.47)
mu=1;sigma=1;p0=0;var1=NULL;p=NULL;pn=NULL
If you are going to passing values to the above using indices then
they should be defined as vectors of the appropriate length. It would
not be absolutely necessary to to so , but if you don't, then you need
to extend them using the c() function, which you do not seem to have
caught onto yet.
for (i in 2:length(n)){
var1[i] = sd[i]^2 #sample variance
p[i] = n[i]/var1[i] #sample precision
p0=p[i-1]
pn[i] = p0[i]+p[i]
mu[i] = ((p0[i])/(pn[i])*mu[i])+((p[i])/(pn[i]))*ybar[i]
sigma[i] = 1/sqrt(pn[i])
mutotals[i,]-mu[i]
}
Joe King
206-913-2912
j...@joepking.com
Never throughout history has a man who lived a life of ease left a
name
worth remembering. --Theodore Roosevelt
-Original Message-
From: Joe P King [mailto:j...@joepking.com]
Sent: Sunday, July 18, 2010 6:13 AM
To: 'r-help@r-project.org'
Subject: RE: [R] loop troubles
I tried that, this is what I tried, but I only get it to do one
iteration
and then it wont cycle back. I am not sure how to tell it how to
look at the
previous precision to add to the current new sample precision to get
the new
precision. This is my first try at a loop so I am not sure if I am
doing
anything right, sorry about the elementary question.
n=c(5,10,15)
ybar=c(12,13,14)
sd=c(2,3,4)
mutotals-matrix(nrow=length(n), ncol=1)
mu=1;p0=0
for (i in 1:length(n)){
var = sd[i]^2 #sample variance
p = n[i]/var #sample precision
p0[i]= i
pn = p0[i]+p[i]
mu.out = ((p0[i])/(pn)*mu[i])+((p)/(pn))*ybar[i]
sigma[i] = 1/sqrt(pn[i])
mutotals[i,]-mu.out[i]
}
Joe King
206-913-2912
j...@joepking.com
Never throughout history has a man who lived a life of ease left a
name
worth remembering. --Theodore Roosevelt
-Original Message-
From: David Winsemius [mailto:dwinsem...@comcast.net]
Sent: Sunday, July 18, 2010 5:36 AM
To: Joe P King
Cc: r-help@r-project.org
Subject: Re: [R] loop troubles
On Jul 17, 2010, at 9:09 PM, Joe P King wrote:
Hi all, I appreciate the help this list has given me before. I have a
question which has been perplexing me. I have been working on doing a
Bayesian calculating inserting studies sequentially after using a
non-informative prior to get a meta-analysis type result. I created a
function using three iterations of this, my code is below. I insert
prior
mean and precision (I add precision manually because I want it to be
zero
but if I include zero as prior variance I get an error cant divide
by zero.
Now my question is instead of having the three iterations below, is
there
anyway to create a loop, the only problem is I want to have elements
from
the previous posterior to be the new prior and now I cant figure out
how to
do the code below in a loop.
Perhaps if you set the problem up with vectors, it would become more
obvious how to do it with indexing or loops:
n - c(5, 10, 15)
ybar - c(12,12,14) # and so on...
--
David
The data below is dummy data, I used a starting
mu of 1, and starting precision of 0.
bayes.analysis.treat-function(mu0,p0){
n1 = 5
n2 = 10
n3 = 15
ybar1 = 12
ybar2 = 13
ybar3 = 14
sd1 = 2
sd2 = 3
sd3 = 4
#posterior 1
var1 = sd1^2 #sample variance
p1 = n1/var1 #sample precision
p1n = p0+p1
mu1 = ((p0)/(p1n)*mu0)+((p1)/(p1n))*ybar1
sigma1 = 1/sqrt(p1n)
#posterior 2
var2 = sd2^2 #sample variance
p2 = n2/var2 #sample precision
p2n = p1n+p2
mu2 = ((p1n)/(p2n)*mu1)+((p2)/(p2n))*ybar2
sigma2 = 1/sqrt(p2n)
#posterior 3
var3 = sd3^2 #sample variance
p3 = n3/var3 #sample precision
p3n = p2n+p3
mu3 = ((p2n)/(p3n)*mu2)+((p3)/(p3n))*ybar3
sigma3 = 1/sqrt(p3n)
posterior-cbind(rbind(mu1,mu2,mu3),rbind(sigma1,sigma2,sigma3))
rownames(posterior)-c(Posterior 1, Posterior 2, Posterior 3)
colnames(posterior)-c(Mu, SD)
return(posterior)}
David Winsemius, MD
West Hartford, CT
__
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: package plotrix
On Jul 18, 2010, at 1:32 PM, mau...@alice.it wrote: -Messaggio originale- Da: Uwe Ligges [mailto:lig...@statistik.tu-dortmund.de] Inviato: dom 18/07/2010 18.58 A: mau...@alice.it Cc: j...@bitwrit.com.au; r-h...@stat.math.ethz.ch Oggetto: Re: [R] package plotrix On 18.07.2010 17:59, mau...@alice.it wrote: I installed package plotrix because reading its vignette it looks like it can help me solve a legend problem. The package instaleed correctly on my Mac OS/X 10.5.8 But I cannot reproduce the examples centered on function lgendg. You mean legendg. Yes. Sorry for my typo. library(plotrix) plot(0.5,0.5,xlim=c(0,1),ylim=c(0,1),type=n, + main=Test of grouped legend function) legendg(0.5,0.8,c(one,two,three),pch=list(1,2:3,4:6), + col=list(2,3:4,5:7)) Error: could not find function legendg That works for me. Try to upgrade bopth R and plotrix if you have not already. I am running R 2.9.2 on my Mac OS/X 10.5.8 I always feel uncomfortable with upgrades. Shall I uninstall R 2.9.2 in advance of instaling the ew bversion for Mac ? No need to do so. The appropriate 2.11.1 folders will be created in the R.framework folder and the links in the Current folder will be updated. I ave no reason for keeping two versions. My problem with the regular R legend function is that I cannot indicate in the legend the line plotted by the command abline Here is the code. Attached is the plot. Any suggestion is welcome. Thank you. Maura x11(width=10,heigh=8) plot(NN,LB,xlim=c(1,300),ylim=c(0,3), main=Intrinsic Dimensionality of 1D Helix,font.main=2,cex.main=2,col.main=red, xlab=Number of Nearest-Neighbors,ylab=Estimated Dimension,col.lab=red,cex.lab=1,font.lab=2, xaxt=n,yaxt=n,pch=19,col=red4) axis (1 ,at = NN [1 : 40 ],labels = NN [1 :40],cex.lab=1,cex.axis=0.8,col.axis=blue,lwd.ticks=0.5,col.ticks =gray3) axis (2 ,at = c (0,0.5,1,1.5,2,2.5,3,3.5,4 ),labels=c(0,0.5,1,1.5,2,2.5,3,3.5,4),cex.lab=1,cex.axis=0.8, col.axis=blue,lwd.ticks=0.5,col.ticks =gray3) lines(NN,McG,pch=18,col=green) lines(NN,PBJD,pch=20,col=magenta1) points(NN,GRPR,pch=18,col=cyan) abline(h=1,pch= 1,lty=1,col=black) legend(topright,legend = c(Levina-Bickel,MacKay- Ghahramani,Pettis-Bailey-Jain-Dubes,Grassberger- Procaccia,Helix Dimension), text.width = strwidth(Pettis-Bailey-Jain-Dubes),pch = c(19,18,20,18,1), col = c(red4,green,magenta1,cyan,black),text.col = black,lty=c(-1,-1,-1,1), bg =snow) We do not have the data nor do we see the figure, hence this is not reproducible for us. Actually I attached the TIFF plot that I am attaching again. Re-read the Posting Guide. It specifies what can be attached and TIFF files are NOT in the list of acceptable types. Uwe might have gotten a copy but the rest of us did not. -- David Winsemius, MD West Hartford, CT __ 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] specifying column names in a vector of characters and the use?
On Jul 18, 2010, at 10:09 PM, Seth wrote:
Hi,
What I would like to do is have a data.frame with column names and
have
these column names stored as strings in another vector. Then I
would like
to be able to access the data.fram columns via referencing the
vector of
names. The code below shows the last few executions that failed to
retrieve
the values for column named X1. Seth
table.1-cbind(c(1,2,3,2,2),c(0,9,0,7,9),c(7,5,9,8,8))
table.1
[,1] [,2] [,3]
[1,]107
[2,]295
[3,]309
[4,]278
[5,]298
table.1-data.frame(table.1)
table.1
X1 X2 X3
1 1 0 7
2 2 9 5
3 3 0 9
4 2 7 8
5 2 9 8
hold-c(X1,X2,X3)
hold
[1] X1 X2 X3
table.1$X1
[1] 1 2 3 2 2
hold[1]
[1] X1
table.1$hold[1] # FROM HERE DOWN ARE MY ATTEMPTS TO ACCESS X1
NULL
Try instead:
table.1[ , hold[1] ]
The $ formalism does not evaluate its argument, but the [ function
does.
--
David.
table.1$(hold[1])
Error: unexpected '(' in table.1$(
table.1$get(hold[1])
Error: attempt to apply non-function
table.1$(get(hold[1]))
Error: unexpected '(' in table.1$(
--
View this message in context:
http://r.789695.n4.nabble.com/specifying-column-names-in-a-vector-of-characters-and-the-use-tp2293494p2293494.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] Grouping and stacking bar plot for categorical variables
On Jul 19, 2010, at 9:36 AM, Simon Kiss wrote: Hi all, I have a series of cateogiral variables that look just like this: welfare=sample(c(less, same, more), 1000, replace=TRUE) education=sample(c(less, same, more), 1000, replace=TRUE) defence=sample(c(less, same, more), 1000, replace=TRUE) egp=sample(c(salariat, routine non-manual, self-employed, farmers, skilled labour, foremen, unskilled labour, social and cultural specialists), 1000, replace=TRUE) welfare, education and defence are responses to a series of questions about whether or not the respondent supports, less, the same or more spending on an issue. egp is a class category. What I would like is a barplot that is both stacked and grouped. The x-axis categories should be the egp class category. Within each class category I would like a cluster of stacked bars that show the distribution of spending support for each issue. Can anyone suggest something? Learn to search: RSiteSearch(stacked barchart) Once you are there you can also ask for prior years' r-help searching which will provide a large number of worked examples since this has been a frequently asked (and answered) question. -- David. Yours, Simon Kiss * Simon J. Kiss, PhD SSHRC and DAAD Post-Doctoral Fellow John F. Kennedy Institute of North America Studies Free University of Berlin Lansstraße 7-9 14195 Berlin, Germany Cell: +49 (0)1525-300-2812, Web: http://www.jfki.fu-berlin.de/index.html __ 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, MD West Hartford, CT __ 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] invalid type error
On Jul 19, 2010, at 3:43 PM, Erik Iverson wrote: jd6688 wrote: myDF = data.frame(id=c(A10,A20),d1=c(.3,.3),d2=c(.4,.4),d3=c(-. 2,.5),d4=c(-.3,.6),d5=c(.5,-.2),d6=c(.6,-.4),d7=c(-.9,-.5),d8=c(-. 8,-.6)) doit=function(x)c(x[1],sum_LK_positive=sum(x[-1] [x[-1]0]),sum_LK_negative=sum(x[-1][x[-1]0])) myDF id d1 d2 d3 d4 d5 d6 d7 d8 1 A10 0.3 0.4 -0.2 -0.3 0.5 0.6 -0.9 -0.8 2 A20 0.3 0.4 0.5 0.6 -0.2 -0.4 -0.5 -0.6 t(apply(myDF,1,doit)) Error in sum(x[-1][x[-1] 0]) : invalid 'type' (character) of argument: I changed the id=c(100,101) in myDF, it worked. are there any way to have this working if the id=c(a10,a20)? doit - function(x)c(sum_LK_positive=sum(x[x0]),sum_LK_negative=sum(x[x0])) cbind(myDF, t(apply(myDF[-1],1,doit))) And the reason that works (which I'm sure Erik knows) is that it does not pass the entire row which would result in coercion of the row vector to the lowest common type ... which in this case would be character. And I think the OP might have wanted this: cbind(myDF[,1,drop=FALSE], t(apply(myDF[-1],1,doit))) id sum_LK_positive sum_LK_negative 1 A10 1.8-2.2 2 A20 1.8-1.7 And the reason the ,drop=FALSE is needed is that otherwise the vector and the matrix would have been non-conformable after the column vector became a row vector. -- David Winsemius, MD West Hartford, CT __ 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] Indexing by logical vectors
On Jul 19, 2010, at 7:16 PM, Christian Raschke wrote: Dear R-Listers, My question concerns indexing vectors by logical vectors that are based on the original vector. Consider the following simple example to hopefully make clear what I mean: a - rnorm(10) a[a0] - NA However, I am now working with multiple data frames that I received, where each of them has nicely descriptive, yet long names(). In my scripts there are many instances where operations similar to the one above are required. Again a simple example: some.data.frame - data.frame(some.long.variable.name=rnorm(10), some.other.long.variable.name=rnorm(10)) some.data.frame$some.other.long.variable.name[some.data.frame $some.other.long.variable.name 0] - NA The fact that the names are so long makes things not very readable in the script and hard to debug. Is there a way in R to refer to the self of whatever is being indexed? I am looking for something like some.data.frame$some.other.long.variable.name[.self 0] - NA There is an alternative, is.na()- which I think is a bit more readable: is.na($some.other.long.variable.name) - some.data.frame $some.other.long.variable.name 0 But do _not_ do: with(some.data.frame, is.na(some.other.long.variable.name) - some.other.long.variable.name 0 ) -- David. that would accomplish the same result as above. Or is there another concise, but less messy way to do this? I prefer not attaching the data.frames and partial matching makes things even more messy since many names() are very similar. I know I could just rename everything, but I'd like to learn if there is and easy or obvious way to do this in R that I have missed so far. I would appreciate any advice, and I apologize if this topic has been discussed before. sessionInfo() R version 2.11.0 (2010-04-22) x86_64-redhat-linux-gnu locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=C LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base -- Christian Raschke Department of Economics and ISDS Research Lab (HSRG) Louisiana State University cras...@lsu.edu __ 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, MD West Hartford, CT __ 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] ACCTGMX to 1223400 in R?
On Jul 19, 2010, at 5:31 PM, John1983 wrote: Hi, I am a newbie in R and was working on some DNA data represented as strings of A,C,T and G (also wild-character like M and X). I use the Bioconductor package in R. Well, I guess it's sort of a meta package, but it is really more of a subculture. It also has its own mailing list. Currently I need to convert a string of the form ACCTGMX to 1223400 i.e. A is replaced by 1, C with 2, T with 3, G with 4 and any other character with a 0. I checked with 'replace' and also with a function called 'copySubstitute' found in the Biobase package but this is only for files. The data here is a string (ACCTGMX ) and we need to convert it to yet another string (1223400). Now I use the strsplit function to split ACCTGM into A C C T G M and then use 'which' to assign the corresponding numbers. Is there a faster way to do this or some function I can make use of? tst - rep( ACCTGMX, 5) newtst - gsub(A, 1, tst) newtst - gsub(C, 2, newtst) newtst - gsub(T, 3, newtst) newtst - gsub(G, 4, newtst) newtst - gsub([[:alpha:]], 0, newtst) newtst [1] 1223400 1223400 1223400 1223400 1223400 There is also a rollaply function in teh zoo and an strapply function in the gsubfn package that might be even more powerful, but I am insufficiently talented to give you a one-liner using them. Please advise. Thank you. -- -- David Winsemius, MD West Hartford, CT __ 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] Constrain density to 0 at 0?
On Jul 19, 2010, at 9:57 PM, Farley, Robert wrote: I'm plotting some trip length frequencies using the following code: plot( density(zTestData$Distance, weights=zTestData$Actual), xlim=c(0,10), main=Test TLFD, xlab=Distance, col=6 ) lines(density(zTestData$Distance, weights=zTestData$FlatWeight), col=2) lines(density(zTestData$Distance, weights=zTestData$BrdWeight ), col=3) which works fine except the distances are all positive, but the densities don't drop to 0 until around -2 or -3. Is there a way for me to force the density plot to 0 at 0? Yes. (Assuming it can be zero, given the data.) Read the help page for density more carefully. Especially the bw and from arguments. -- 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.
Re: [R] trouble getting table of coeffs with quantreg with fixed effects
On Jul 19, 2010, at 9:37 PM, Steve McDonald wrote: I'm a new user, so my apologies for what is likely a dumb question... I am having a hard time getting a table of regression results when using Koenker's code for quantile regression with fixed effects (http://www.econ.uiuc.edu/~roger/research/panel/rq.fit.panel.R ). I use the example data parameters that Koenker provides (see below). m - 3 n - 10 s - rep(1:n,rep(m,n)) x - exp(rnorm(n*m)) X - cbind(1,x) u - x*rnorm(m*n) + (1-x)*rf(m*n,3,3) a - rep(rnorm(n),rep(m,n)) y - a + u fit - rq.fit.panel(X,y,s) When I use summary(fit), I get the following output... Length Class Mode coef 16 -none- numeric ierr 1 -none- numeric it1 -none- numeric time 1 -none- numeric When I use list(fit), I can't tell which coefficients correspond to which parameters... [[1]]$coef [1] 7.701321 -2.608782 7.664906 -1.866059 7.894542 [6] -1.878973 -6.117778 -1.570840 -5.380230 -5.815231 [11] -7.191127 -5.671467 -1.712052 -5.055224 -5.623140 [16] -5.617846 Is there a way to get a table of results or a way of discerning what these coeffs refer to? Thanks in advance. Did you read the documentation in the rq.fit.panel function? 2. On return the coefficient vector has m*p + n elements where m is the number + # quantiles being estimated, p is the number of colums of X, and n is the + # number of distinct values of s. The first m*p coefficients are the + # slope estimates, and the last n are the fixed effects -- David Winsemius, MD West Hartford, CT __ 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] Indexing by logical vectors
On Jul 20, 2010, at 1:12 AM, Christian Raschke wrote:
On Tue, 2010-07-20 at 10:12 +1000, bill.venab...@csiro.au wrote:
As far as I know the answer to your question is No, but there are
things you can do to improve the readability of your code. One
thing I find useful is to avoid using $ as much as possible and
to favour things like with() and within().
Thank you for your answer. I had not looked at within() for this until
now.
The first thing you might do is think about choosing shorter names,
of course. If that's not possible, you could try something like
this.
ensureNN - function(x) { # ensure non-negative
is.na(x[x 0]) - TRUE
x
}
This approach would essentially require a different function for the
different operations to be performed on the vector. I suppose that
assigning NA based on a certain condition is probably the most common
use, but in the end I have other cases, where the logical vector is
obtained from other operations or where the value that is assigned is
different case by case; for example,
levels(something.long)[levels(something.long) %in% LETTERS[1:3]] -
Z
So given that your general answer above to my inquiry was No, I will
keep experimenting and I'll also have another look at with() and
within().
You might want to look at the sqldf package. I have noted over the
year or two since it was released that it is sometimes possible to do
rather amazing operations with minimal code. The sort of operations
you anticipate (transformations dependent on logical criteria) seem
to be a good candidate for a database oriented syntax.
--
David.
Thanks again!
some.data.frame - within(some.data.frame, {
some.long.variable.name - ensureNN(some.long.variable.name)
some.other.long.variable.name -
ensureNN(some.other.long.variable.name)
})
Of course if you wanted to do this to all variables in a data frame
you could do
some.data.frame - data.frame(lapply(some.data.frame, ensureNN))
and it all happens, no questions asled. (I can see a generic
function emerging here, perhaps...)
W.
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org
] On Behalf Of Christian Raschke
Sent: Tuesday, 20 July 2010 9:16 AM
To: r-help@r-project.org
Subject: [R] Indexing by logical vectors
Dear R-Listers,
My question concerns indexing vectors by logical vectors that are
based
on the original vector. Consider the following simple example to
hopefully make clear what I mean:
a - rnorm(10)
a[a0] - NA
However, I am now working with multiple data frames that I received,
where each of them has nicely descriptive, yet long names(). In my
scripts there are many instances where operations similar to the one
above are required. Again a simple example:
some.data.frame - data.frame(some.long.variable.name=rnorm(10),
some.other.long.variable.name=rnorm(10))
some.data.frame$some.other.long.variable.name[some.data.frame
$some.other.long.variable.name
0] - NA
The fact that the names are so long makes things not very readable in
the script and hard to debug. Is there a way in R to refer to the
self
of whatever is being indexed? I am looking for something like
some.data.frame$some.other.long.variable.name[.self 0] - NA
that would accomplish the same result as above. Or is there another
concise, but less messy way to do this? I prefer not attaching the
data.frames and partial matching makes things even more messy since
many
names() are very similar. I know I could just rename everything,
but I'd
like to learn if there is and easy or obvious way to do this in R
that I
have missed so far.
I would appreciate any advice, and I apologize if this topic has been
discussed before.
sessionInfo()
R version 2.11.0 (2010-04-22)
x86_64-redhat-linux-gnu
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=C LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
__
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, MD
West Hartford, CT
__
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 time in R
On Jul 20, 2010, at 7:33 AM, Sarah Chisholm wrote: Hi, I have a problem with the time formatting in R. I have entered time in the format MM:SS.xyz and R has automatically classified this as a factor, but I need it numerically. However when I use as.numeric() it gives me totally different numbers. Is there any way I can tell R to read thes input as a number? There are ways using As methods of doing so, but I think those is a bit more complex than a learner at you level should be asked to implement. Guessing that you used one of the read functions to bring in the data (or possibly used hte data.frame function) , you are probably experiencing the effects of the default stringAsFactors=TRUE behavior. Try setting stringsAsFactors to FALSE. Then you can use the standard Date functions: ?Dates -- David Winsemius, MD West Hartford, CT __ 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 time in R
On Jul 20, 2010, at 8:41 AM, David Winsemius wrote: On Jul 20, 2010, at 7:33 AM, Sarah Chisholm wrote: Hi, I have a problem with the time formatting in R. I have entered time in the format MM:SS.xyz and R has automatically classified this as a factor, but I need it numerically. However when I use as.numeric() it gives me totally different numbers. Is there any way I can tell R to read thes input as a number? There are ways using As methods of doing so, but I think those is a bit more complex than a learner at you level should be asked to implement. Guessing that you used one of the read functions to bring in the data (or possibly used hte data.frame function) , you are probably experiencing the effects of the default stringAsFactors=TRUE behavior. Try setting stringsAsFactors to FALSE. Then you can use the standard Date functions: ?Dates I should have said: ? DateTimeClasses -- David Winsemius, MD West Hartford, CT __ 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] Means from selected columns in a data frame
On Jul 20, 2010, at 9:13 AM, Marcus Drescher wrote: Hi all, I have a dataframe with survey data. Now I want to calculate means from several but not all columns (e.g. a1, a2, a3) and save them in a new separate column (e.g. a). Well that would be possible unless you are in Minitab or Excel where you can put stuff of variable length where ever. Like: a = mean(a1, a2, a3) Can someone help me or give me the right key word to look for? ?colMeans Perhaps: meansvec - colMeans(dfrm[ , c(a1, a2, a3)] ) If on the other hand, you wanted rowMeans (since your problems specification was rather ambiguous) that function is also available. -- David Winsemius, MD West Hartford, CT __ 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] Means from selected columns in a data frame
On Jul 20, 2010, at 9:25 AM, David Winsemius wrote: On Jul 20, 2010, at 9:13 AM, Marcus Drescher wrote: Hi all, I have a dataframe with survey data. Now I want to calculate means from several but not all columns (e.g. a1, a2, a3) and save them in a new separate column (e.g. a). Well that would be possible unless you are in Minitab or Excel where you ^not^ can put stuff of variable length where ever. But I probably misunderstood your intent, anyway Like: a = mean(a1, a2, a3) Can someone help me or give me the right key word to look for? ?colMeans Perhaps: meansvec - colMeans() If on the other hand, you wanted rowMeans (since your problems specification was rather ambiguous) that function is also available. I suppsoe your request that these results be put in another column may have meant you wanted row means and therefore you could try: dfrm$a_mns - rowMeans( dfrm[ , c(a1, a2, a3)] ) David Winsemius, MD West Hartford, CT __ 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] Convert Unix (Epoch) timestamp to DD/MM/YY and HH:MM:SS
On Jul 20, 2010, at 9:26 AM, Jim Hargreaves wrote: Dear List, After much searching with no success, I would like to ask how I can convert a unix/POSIX time (seconds since Jan 01, 1970) into a string like 01/01/1970 00:00 This is probably easily done with a system(date...) but it would be great if I could do it in R. By default DateTime objects are printed similarly to your specification: Sys.time() [1] 2010-07-20 09:34:36 EDT as.numeric( Sys.time() ) # the internal representation [1] 1279632906 The default format for the default print function is -MM-DD HH:MM TZ. You can find the format codes with: ?strptime print(Sys.time(), format=%m/%d/%Y %H:%M) # your request is ambiguous w.r.t month and day order [1] 07/20/2010 09:39 EDT -- David Winsemius, MD West Hartford, CT __ 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] Nesting functions in loops that result in error messages breaking the loop
On Jul 20, 2010, at 9:44 AM, Patrick McKann wrote:
Hello all,
I am trying to write a program in R in which I call a function
multiple
times within a loop. The problem is that sometimes the function
breaks down
while calling another function, and produces an error message that
breaks my
loop and the program stops. I would like to keep the loop running
when this
function breaks down, and just move on to the next iteration in the
loop.
Is there any way to buffer the output of a function within the loop,
so that
I can note that the function produced an error message, without the
error
message breaking the loop and stopping my program? Let me know if
this
question does not make sense.
?try
This is a worked example that Josh Barnett posted a couple of days ago
after some (lengthy) coaching:
x - read.table(textConnection(y1 y2 y3 x1 x2
indv.1 bagels donuts bagels 4 6
indv.2 donuts donuts donuts 5 1
indv.3 donuts donuts donuts 1 10
indv.4 donuts donuts donuts 10 9
indv.5 bagels donuts bagels 0 2
indv.6 bagels donuts bagels 2 9
indv.7 bagels donuts bagels 8 5
indv.8 bagels donuts bagels 4 1
indv.9 donuts donuts donuts 3 3
indv.10 bagels donuts bagels 5 9
indv.11 bagels donuts bagels 9 10
indv.12 bagels donuts bagels 3 1
indv.13 donuts donuts donuts 7 10
indv.14 bagels donuts bagels 2 10
indv.15 bagels donuts bagels 9 6), header = TRUE)
closeAllConnections()
results - matrix(nrow = 1, ncol = 3)
colnames(results) - c(y1, y2, y3)
require(rms) # or Design
for (i in 1:3) {
mod.poly3 - try(lrm(x[,i] ~ pol(x1, 3) + pol(x2, 3), data=x),
silent=TRUE)
if(class(mod.poly3)[1] != 'try-error')
{results[1,i] - anova(mod.poly3)[1,3]}
David Winsemius, MD
West Hartford, CT
__
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] Servreg $loglik
On Jul 20, 2010, at 11:20 AM, Charles Annis, P.E. wrote: Dear R-experts: I am using survreg() to estimate the parameters of a Weibull density having right-censored observations. Some observations are weighted. To do that I regress the weighed observations against a column of ones. When I enter the data as 37 weighted observations, the parameter estimates are exactly the same as when I enter the data as the corresponding 70 unweighted observations. This is to be expected, of course. I don't understand, however, why the reported loglikelihood is parameter.estimates$loglik [1] -120.4699 -120.4699 for the 37 weighted observations, but parameter.estimates$loglik [1] -135.1527 -135.1527 for the 70 unweighted observations. This has come up on r-help many times before (and probably on other lists as well), despite not being an R question at all. It is commonplace in modeling grouped data to see likelihoods reported differently from the result obtained when modeling ungrouped data representations with the same frequencies. The only valid statistical process is to compare differences in the likelihoods (or log(L) ), since the likelihood (or log(L) ) is only defined up to an arbitrary constant. You need to be comparing the result to some sort of null model for it to have any meaning. (... or perhaps that is your null model and you need to be looking at the impact of adding a covariate or two.) -- David. (For the record, my computations of the loglikelihood, using the dweibull() function for the observations and the pweibull() function for the censored observations, is -135.1527 for both 37 weighted and 70 unweighted.) I am using the data from Meeker and Escobar, _Statistical Methods for Reliability Data_, Wiley (1998), Table C.1, shown below: Hours Status Num.Parts 450Failure 1 460R-Censored 1 1150Failure 2 1560R-Censored 1 1600Failure 1 1660R-Censored 1 1850R-Censored 5 2030R-Censored 3 2070Failure 2 2080Failure 1 2200R-Censored 1 3000R-Censored 4 3100Failure 1 3200R-Censored 1 3450Failure 1 3750R-Censored 2 4150R-Censored 4 4300R-Censored 4 4600Failure 1 4850R-Censored 4 5000R-Censored 3 6100R-Censored 3 6100Failure 1 6300R-Censored 1 6450R-Censored 2 6700R-Censored 1 7450R-Censored 1 7800R-Censored 2 8100R-Censored 2 8200R-Censored 1 8500R-Censored 3 8750R-Censored 2 8750Failure 1 9400R-Censored 1 9900R-Censored 1 10100 R-Censored 3 11500 R-Censored 1 I am running R version 2.11.1 (2010-05-31) on a HP Windows 7 box with 8 gig RAM. Thank you for your help. Charles Annis, P.E. charles.an...@statisticalengineering.com 561-352-9699 http://www.StatisticalEngineering.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, MD West Hartford, CT __ 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] If help
On Jul 20, 2010, at 1:14 PM, Heiman, Thomas J. wrote:
Hi Y'all,
I have some data in a table with 2 columns. There are two values:
Reduction and No Reduction. I am trying to make a new variable
change which recode the combinations from column 1 and 2 into a
single number. Here is a snippet from the table:
[1,] NoReduction NoReduction
[2,] ReductionReduction
[3,] NoReduction NoReduction
[4,] NoReduction NoReduction
[5,] ReductionReduction
[6,] ReductionReduction
[7,] ReductionReduction
[8,] NoReduction NoReduction
[9,] NoReduction NoReduction
[10,] NoReduction NoReduction
This is the code that I have written so far..
for (i in 1:nrow(change20082009))
if(change20082009[i,1]=='No Reduction' change20082009[i,2]=='No
Reduction') ){change20082009[i,3] - 1} else
if(change20082009[i,1]=='No Reduction' change20082009[i,
2]=='Reduction'){change20082009[i,3] - -1} else
if(change20082009[i,1]=='Reduction' change20082009[i,2]=='No
Reduction') {change20082009[i,3] - 2} else
if(change20082009[i,1]=='Reduction' change20082009[i,
2]=='Reduction') {change20082009[i,3] - 0}
)
I can't seem to get the code above to work..Any suggestions (I am
sure it is really basic)? Is there a better way to do this?
You are going to have a problem if you think that NoReduction == 'No
Reduction'
--
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.
Re: [R] Registered / trademark signs
Works on OSX 10.5.8 + R2.11.1
--
David.
On Jul 20, 2010, at 1:36 PM, Henrique Dallazuanna wrote:
I don't know if this works in OSX, but in XP:
plot(0)
title('\u00ae - \u2122')
On Tue, Jul 20, 2010 at 2:22 PM, Dennis Fisher
fis...@plessthan.com wrote:
Colleagues,
What is the easiest means to embed a:
® (registered)
or
™ (trademark)
sign in text in a graphic. I would like to use mtext and avoid
plotmath,
if possible. Ideally, the sign should be superscripted but I can
easily
sacrifice that.
Optimally, I need a solution that works in both OS X and Windows
(≥ XP) and
with R versions ≥ 2.11
Thanks in advance.
__
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] square brackets in expression in plots
On Jul 20, 2010, at 12:42 PM, Jannis wrote:
Dears,
do you know whether it is possible to include any square parantheses
(brackets) in an expression to use it as an axis label?
e.g. I would like to have a label like (with the sub and
superscripts correct):
speed [m s^-1]
Not exactly clear what you mean by correct.
I know how to combine an expression with text via paste, but as I
run soimething like:
a='m*s^{-1}'
plot(1:10,main=parse(text=a))
?plotmath
This seemed to work fine with main, sub and xlab:
a='speed [m*s^-1]'
plot(1:10,main=parse(text=a))
But maybe you wanted the square-brackets?
a='speed~bgroup([, m*s^-1, ])'
plot(1:10,main=parse(text=a))
I found now way of doing this. I use the parse thing as I have all
these units stored as strings that represent expressions.
Cheers for any help!
Jannis
__
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, MD
West Hartford, CT
__
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] Limited output
On Jul 20, 2010, at 1:05 PM, confusedcius wrote: Hi there, I entered over 2000 lines of data into R from SQLite How? and saved it in R as a data.frame, How? but the data.frame only gives the first 500 lines. Is there any way to either increase the default limit or to determine which of the original lines should be in the output (e.g. maybe the last 500 instead of the first 500)? There is no limit on the size of data.frames, Well, at least not anything near 500, anyway. I suspect that the limit is set by the maximum size of the (positive) integer data-type = 2*10^9 What does this return? options()$max.print # [1] 9 # which is the default And what does this return? str(your.data.frame) -- David Winsemius, MD West Hartford, CT __ 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] Constrain density to 0 at 0?
On Jul 20, 2010, at 2:45 PM, Farley, Robert wrote: Read the help page for density more carefully. Especially the bw and from arguments. Yes, it was my inability to make sense of the help page that motivated my email. My distances range from 0.4 to 7.6 but the density plot ranges from -2 to 10. from=0 seems to hide the negative portion of the density without setting it to 0. I've tried adjust, but it requires a value of 0.2, which results in a very, very spiky plot. I was hoping for a process that could forbid impossible (negative) values, yet allow the density plot to blend the individual measurements. Is there a variable bw that could be set small at the extrema, and larger in the range of the data? The best answer may depend on what your (unstated) goals are. Have you considered fitting a log-Normal or Gamma distribution to the data? You could plot the histogram and then overlay a smooth fitted distribution curve. ?hist require(MASS) ?fitdistr -- David. -Original Message- From: David Winsemius [mailto:dwinsem...@comcast.net] Sent: Monday, July 19, 2010 19:31 To: Farley, Robert Cc: r-help@r-project.org Subject: Re: [R] Constrain density to 0 at 0? On Jul 19, 2010, at 9:57 PM, Farley, Robert wrote: I'm plotting some trip length frequencies using the following code: plot( density(zTestData$Distance, weights=zTestData$Actual), xlim=c(0,10), main=Test TLFD, xlab=Distance, col=6 ) lines(density(zTestData$Distance, weights=zTestData$FlatWeight), col=2) lines(density(zTestData$Distance, weights=zTestData$BrdWeight ), col=3) which works fine except the distances are all positive, but the densities don't drop to 0 until around -2 or -3. Is there a way for me to force the density plot to 0 at 0? Yes. (Assuming it can be zero, given the data.) Read the help page for density more carefully. Especially the bw and from arguments. -- 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. David Winsemius, MD West Hartford, CT __ 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] blank pdf output when called in a loop
On Jul 20, 2010, at 7:21 PM, Joshua Wiley wrote:
Hi Nicolas,
You nee to explicitly wrap it in print() when it is inside a loop (if
I'm not mistaken also when inside a function). With lattice loaded,
you can find the specific print methods by methods(print) .
The interpreter handles finding the proper print methods. You just
need to have a clear list in your mind regarding what is base graphics
(e.g.plot) and what is grid graphics (e.g. xyplot and qplot). This is
part of the FAQ which points out that this also applies to any
source()-ed code:
http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-do-lattice_002ftrellis-graphics-not-work_003f
(which points out that this applies to ggplot2 functions as well.)
To further confuse the issue, not all functions named plot are
dispatched to base graphics. Some are dispatched to lattice graphics.
--
David.
From your
example:
pdf(temp1.pdf, width=8, height=8)
for(k in 1) {
print(wireframe(volcano))
}
dev.off()
HTH,
Josh
On Tue, Jul 20, 2010 at 4:02 PM, Nicolas STRANSKY
stran...@broadinstitute.org wrote:
Hi,
I'm hitting a strange problem where pdf plots that I'm trying to make
are blank, only when produced from within a loop. The pdf contains
0 page.
I've narrowed the problem to this minimal script that invariably
produces an empty pdf with my setup:
pdf(/local/scratch/1.pdf, width=8, height=8)
for (k in 1) {
wireframe(volcano)
}
dev.off()
The odd thing is that
pdf(/local/scratch/1.pdf, width=8, height=8)
wireframe(volcano)
dev.off()
works fine!
Am I doing something wrong here? I've tried on two different systems,
Linux or Mac.
Thanks,
Nico
R version 2.11.0 (2010-04-22)
x86_64-unknown-linux-gnu
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
LC_TIME=en_US.UTF-8LC_COLLATE=en_US.UTF-8 LC_MONETARY=C
[6] LC_MESSAGES=en_US.UTF-8LC_PAPER=en_US.UTF-8 LC_NAME=C
LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lattice_0.18-5 custom_1.1
loaded via a namespace (and not attached):
[1] grid_2.11.0
--
Nicolas STRANSKY, Ph.D.
Computational Biologist, Cancer Program
Broad Institute of MIT and Harvard
5CC-1339 | 7 Cambridge Center | Cambridge, MA 02142, USA
Phone: +1 617 714 7564 |n...@broadinstitute.org
__
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.
--
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.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, MD
West Hartford, CT
__
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] best way to apply a list of functions to a dataset ?
as far as style or simple stability-enhancing improvements would be handy. regards, Glen __ 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. David Winsemius, MD West Hartford, CT __ 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] date from weeknumber
On Jul 20, 2010, at 6:37 PM, H Rao wrote: Hi, Is there a function to get the last(or first) day of the week, given the week number of the year? For eg, week number for 7/20 is 29 as obtained by format(Sys.Date(),%U), is there a function which returns 7/25 - the last day of week # 29 require(tis) nthMon - function(x) as.Date(currentMonday(xTi=as.Date(2010-01-01)) +7*(x-1)) nthMon(2) [1] 2010-01-11 Or: nthMonYr - function(n, Yr) as.Date(currentMonday(xTi=as.Date(paste(Yr,-01-01,sep=)))+7*(n-1)) nthMonYr(2,2010) [1] 2010-01-11 TIA, Rao. [[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] Obtaining the unmerged cases from one of the two data set
On Jul 21, 2010, at 5:33 AM, Mangalani Peter Makananisa wrote: Dear R Gurus, I saw no reason to copy Rob Hyndman. I did not see that this involves any of the packages he maintains. I am having two dummy csv data sets A and B containing 19 and 15 cases/observations respectively. From the two data set 13 cases are intersection. From one of the two (any) data set, How do I then retrieve the unmerged data ? let's take A for example, six cases must appear in our results. See the R codes below. ?setdiff Perhaps: setdiff( (NAME(A), NAME(B) ) You can also do a merge that is an outer join that includes all the NAME information and then extract the records with SALARY and .B.SALARY data. Untested in absence of working example: ?merge mer - merge(A,B, all=TRUE) mer[ mer$NAME %in% setdiff(NAME(A), NAME(B) ), ] -- David. A = read.csv(C:/Documents and Settings/S1033067/Desktop/A.csv, header = TRUE, dec =,, sep = ,) names(A) [1] NAME SALARY dim(A)[1] [1] 19 B = read.csv(C:/Documents and Settings/S1033067/Desktop/B.csv, header = TRUE, dec =,, sep = ,) names(B) [1] NAME B.SALARY dim(B)[1] [1] 15 common = merge(A,B) names(common) [1] NAME SALARY B.SALARY dim(common)[1] [1] 13 David Winsemius, MD West Hartford, CT __ 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] Chi-square distribution probability density function:
On Jul 21, 2010, at 6:52 AM, Knut Krueger wrote: Hi to all I found an formular of an ** ***p-Value Calculator for the Chi-Square test* *http://www.danielsoper.com/statcalc/calc11.aspx* *with the formula* *http://www.danielsoper.com/statkb/topic11.aspx* *what's the gamma function of this formula in r?* *df=5* *ch2=25.50878* *the following code does not give the result 0.001 for the values above * *p= ((0.5^(df/2))/gamma(df))*(ch2^((df/2)-1))* (2.718281828459^(- ch2/2)) or is there any other error? Check your implementation of that formula. You made an error in the gamma argument. See also: ?Chisquare ?gamma -- David Winsemius, MD West Hartford, CT __ 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] Contingency Table Analysis - specific cell to specific cell comparisons
On Jul 21, 2010, at 8:07 AM, Tsunhin John Wong wrote: Dear R users, I have a question of how to do some specific cell to cell comparisons on a R x C contingency table. The table is a 3 x 5 table with frequency / count data. langcons.table - table(lang, cons) langcons.table[cbind(lang,cons)] - freq langcons.table Adj Int Oth Pas Tra C 69 221 17 3 198 E 56 214 33 31 174 J 36 291 8 9 164 I know how to do an independent model test using Poisson in glm glm.out1 - glm(freq~lang+cons, family=poisson, data=langcons.data) summary(glm.out1) And then fit the saturated model glm.out2 - glm(freq~lang*cons, family=poisson, data=langcons.data) summary(glm.out2) However, the results are difficult to interpret: C and Adj are used to as a baseline. And I can only see main effects and interactions and *always according to the baseline*. Coefficients: Estimate Std. Error z value Pr(|z|) (Intercept) lang1 lang2 cons1 cons2 cons3 cons4 lang1:cons1 lang2:cons1 lang1:cons2 lang2:cons2 lang1:cons3 lang2:cons3 lang1:cons4 lang2:cons4 If anyone know, please suggest me some way to do specific cell to cell comparison on such a contingency table. Even if you are daunted by the task of plugging the covariates into the formula, exp(intercept+sum(beta_N*var_n)), you can always use the predict function to create an estimate for all (or a specific set) of the covariates. They come out on the log(rate) scale so would need to be exponentiated. Consult your instructor for further help. Say, to compare pairs of cells: along a column: 3 vs 31, 9 vs 31, 3 vs 9 along a row: 36 vs 9 or even across column and row: 36 vs 31, and 36 vs 3 David Winsemius, MD West Hartford, CT __ 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] Error message: attempt to set rownames on object with no dimensions
On Jul 21, 2010, at 9:39 AM, Igor Blanco wrote: Hi there, I am trying to analyze some data using the lme package. In my case there is a variable called Newmarker that can only take 2 values (number 1 or number 2). I have used the as.factor to remark this fact but an error message appears stating: attempt to set rownames on object with no dimensions. Here is the code I have used: *Dataset - read.table(Data2.txt, header=TRUE, sep=, na.strings=NA, dec=., strip.white=TRUE) Dataset$Newmarker - as.factor(Dataset$Newmarker) Why do you use as.factor? It's probably already a factor. What does str(Dataset) show? attach(Dataset) You do not need to do this, and it creates potential for confusion. The data= argument will allow the column names to be accessed by the lme() function. lme.1 - lme(Newmarker ~ Age, data=Dataset, random= ~1|ID) summary(lme.1)* If I don't use the as.factor I have no problems but I don't think this is the correct way to proceed... The other variables I use can have any value. -- David Winsemius, MD West Hartford, CT __ 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] Interactions in GAMMs
On Jul 21, 2010, at 11:17 AM, Karen Moore wrote: Hi, I've an issue adding an interaction to a GAMM: My model was of form: # Package? Probably: require(mgcv) gamm1 - gamm(TOTSR ~ fROT + s(PH) + s(LOI) + s(ASP) + s(SQRT_ELEV) + CANCOV + s(SQRT_TOTCWD) + s(WELLF) + s(WELLN) + s(OLDWDLD) + s(DISTWOOD) + s(Annprec) + s(OLDWDLD:DISTWOOD) + (1|fSITE), family = poisson, data = BIOFOR2) with interaction of s(OLDWDLD:DISTWOOD). However I got this error message below but don't know what it means? The documentation for s() says the unnamed arguments must be a list of covariates and I do not think OLDWDLD:DISTWOOD constitutes a list. (On the other hand it appears that the term list is being used a bit loosely here, since the examples do not explicitly enclose the covariates in the list() function and when I did so with the example on the te help page I get an error.) (my data is composed of info for 150 plots) #I Warning messages: #2: In OLDWDLD:DISTWOOD : # numerical expression has 150 elements: only the first used #3: In OLDWDLD:DISTWOOD : # numerical expression has 150 elements: only the first used Can anyone offer advice on how to include this interaction in GAMM model? Do you have Wood's book? Ch/sect 6.7 has worked examples using the te() function which I believe is the GAM(M) counterpart to linear interactions. He also makes the point that random arguments to gamm() need to be in the list form. (The documentation suggests this has not changes since the book was published.) -- David Winsemius, MD West Hartford, CT __ 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] [SPAM](Aktien) Re: Chi-square distribution probability density function:
On Jul 21, 2010, at 12:36 PM, Knut Krueger wrote: David Winsemius schrieb: On Jul 21, 2010, at 6:52 AM, Knut Krueger wrote: Hi to all I found an formular of an ** ***p-Value Calculator for the Chi-Square test* *http://www.danielsoper.com/statcalc/calc11.aspx* *with the formula* http://www.danielsoper.com/statkb/topic11.aspx *what's the gamma function of this formula in r?* *df=5* *ch2=25.50878* *the following code does not give the result 0.001 for the values above * *p= ((0.5^(df/2))/gamma(df))*(ch2^((df/2)-1))* (2.718281828459^(- ch2/2)) or is there any other error? Check your implementation of that formula. You made an error in the gamma argument. Sorry a copy and paste error but gamma(df/2) does not solve the problem I am not able to get the implementation of the function working :-( See also: ?Chisquare ?gamma there is also a different result using this function from the helpf files and the function /f_n(x) = 1 / (2^(n/2) Gamma(n/2)) x^(n/2-1) e^(-x/2) That is (perhaps) because you are confusing a density function with a cumulative probability function. The expression above should give the same result as a (proper) R transcription of the formula on the page you offered above and should be what dchisq referenced on the ? Chisquare page returns as well. Sure enough... dchisq(25.50878, 5) [1] 4.951e-05 p= ((0.5^(df/2))/gamma(df/2))*(ch2^((df/2)-1))* (2.718281828459^(- ch2/2)) p [1] 4.951e-05 Those formulae only differ in how they use grouping of the constant terms. (In R the function is gamma, not Gamma.) You should be able to get the P(X^2|x25.508, df=5) from an integration of that function. Knut / David Winsemius, MD West Hartford, CT __ 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] Chi-square distribution probability density function:
On Jul 21, 2010, at 1:18 PM, Knut Krueger wrote: David Winsemius schrieb: And exactly why did you think I offered ?Chisquare I was completely on the wrong way, and tried to find a solution with the formula instead to substitute the formula. So I tried to implement pchisq into the formula - and of course I got wrong values ... p - function(x) ((0.5^(5/2))/gamma(5/2))*(x^((5/2)-1))* (2.718281828459^(-x/2)) integrate(p, 25.508, Inf) 0.000 with absolute error 2.7e-05 Check result pchisq(25.508,5) [1] 0. 1-pchisq(25.508,5) [1] 0.000 Thank's Knut David Winsemius, MD West Hartford, CT __ 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] Command that is conditional upon file retrieval: is it possible?
On Jul 21, 2010, at 5:33 PM, AndrewPage wrote:
Hi all,
I'm currently working on an R program where I have to access an FTP
server
to download some of the data I need. However, the people who post
up the
files I access are at times inconsistent with regards to time
posted, if
they post at all, etc Here's some of the code I use:
library(RCurl)
url1 = paste(ftp://user:passw...@a.great.website.com/;, file, num1,
.csv, sep = )
data1 = getURL(url1)
write(data1, file = paste(inMyFolder, num1, .csv, sep = ))
Sometimes this process works perfectly, and sometimes I get an error
message
like this attached to data1 = getURL(url1):
Error in curlPerform(curl = curl, .opts = opts, .encoding
= .encoding) :
RETR response: 550
That's because that particular file isn't on the FTP server yet.
Now...
let's just imagine that there's another way for me to access the file
elsewhere, and I can drag it into the working directory with the
same name
as the file I'm telling R to write immediately after finding it on
the FTP
server.
So here's my question: is it possible to write a command that will
write the
file if there isn't an error message going along with data1 =
getURL(url1), but won't write the file if it can't find it
?try
Untested in absence of workable example...
Return - try(getURL(url1), silent=TRUE) # silent keeps error from
showing up on console
if(grep(550, Return[1]) { # not sure about exact structure of error
message
# but could do nothing
} else { # assign data
data1 - Return
if((550, Return[1]) { # go to the other site
Return2 - try( getURL(url2), silent =TRUE ) # and so on
As of right
now, if I got the error message, dragged the file into the working
directory
and ran the program again, R will overwrite my good file with an
empty one
because in all cases, I'm telling it to write a file with that name
that
includes the information in data1.
Thanks in advance,
Search help archives for try examples or perhaps try-error or tryCatch
which might be more specific.
Andrew
--
--
David Winsemius, MD
West Hartford, CT
__
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] post hoc test for lme using glht ?
On Jul 21, 2010, at 6:36 PM, SNAFU1 wrote: Hi, I have a fairly simple repeated measures-type data set I've been attempting to analyze using the lme function in the nlme package. Repeated searches here and other places lead me to believe I have specified my model correctly. However, I am having trouble with post-hoc tests. From what I gather, other people are successfully using the glht function from the multcomp package to perform post-hoc tests. I've tried multiple iterations, but can't seem to get it to work. Here is (a subset of) what I have been trying. Any help would be greatly appreciated: model.3-lme(fixed=Totnum~Wk*Pop, random=~1|nUID, data=nona) #Wk has 6 levels, 0-5; Pop has 2 levels; nUID is a unique ID given to each individual (multiple individuals from each population were re- measured each week) anova(model.3) numDF denDF F-value p-value (Intercept) 1 592 649.7753 .0001 Week5 592 302.9706 .0001 Pop 1 222 70.7268 .0001 Week:Pop5 592 36.8576 .0001 Come on... are you cutting and pasting fro disjoint sections of your console? The formula says Wk and the parameter name in the anova result is Week. What gives? summary(glht(model.3, linfct=mcp(Wk = Tukey))) Error in mcp2matrix(model, linfct = linfct) : Variable(s) ‘Wk’ of class ‘integer’ is/are not contained as a factor in ‘model’. is.factor(Wk) Well, this object should not even exist at the top level environment unless you attached it and didn't tell us about it. Try: Wk - NULL nona$Wk - factor(nona$Wk) Then re-run your code. [1] FALSE Week-factor(Wk) is.factor(Week) [1] TRUE model.3-lme(fixed=Totnum~Week*Pop, random=~1|nUID) summary(glht(model.3, linfct=mcp(Week = Tukey))) Error in contrMat(table(mf[[nm]]), type = types[pm]) : less than two groups Thanks, Rob -- David Winsemius, MD West Hartford, CT __ 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] xtable
Or class(influencia) and then see it is in this vector:
methods(xtable)
[1] xtable.anova* xtable.aov* xtable.aovlist*
[4] xtable.coxph* xtable.data.frame* xtable.glm*
[7] xtable.lm* xtable.matrix* xtable.prcomp*
[10] xtable.summary.aov* xtable.summary.aovlist* xtable.summary.glm*
[13] xtable.summary.lm* xtable.summary.prcomp* xtable.table*
[16] xtable.ts*
On Jul 21, 2010, at 8:20 PM, John Kane wrote:
Try
str(influencia)
I don't think xtable is intended to print lists.
--- On Wed, 7/21/10, Silvano silv...@uel.br wrote:
From: Silvano silv...@uel.br
Subject: [R] xtable
To: r-help@r-project.org
Received: Wednesday, July 21, 2010, 4:15 PM
Hi,
How do I build a table from a regression model adjusted
using xtable?
Commands are:
modelo1 = lm(Y~X1 + X2)
influencia = influence.measures(modelo1)
require(xtable)
xtable(influencia)
but it isn't work.
Less informative words were never typed. Why did you not include the
error message? Was it because the answer would have been clear?
Error in UseMethod(xtable) :
no applicable method for 'xtable' applied to an object of class
infl
You could have extracted the first element of influencia and used
xtable on the unlisted values.
ctl - c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
trt - c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
group - gl(2,10,20, labels=c(Ctl,Trt))
weight - c(ctl, trt)
anova(lm.D9 - lm(weight ~ group))
summary(lm.D90 - lm(weight ~ group - 1))# omitting intercept
influencia = influence.measures(lm.D9)
str(influencia)
xtable(matrix(unlist(influencia[1]), ncol=6))
% latex table generated in R 2.11.1 by xtable 1.5-6 package
% Wed Jul 21 21:22:30 2010
\begin{table}[ht]
\begin{center}
\begin{tabular}{rrr}
\hline
1 2 3 4 5 6 \\
\hline
1 -0.44 0.31 -0.44 1.02 0.09 0.10 \\
2 0.27 -0.19 0.27 1.15 0.04 0.10 \\
3 0.07 -0.05 0.07 1.24 0.00 0.10 \\
4 0.57 -0.40 0.57 0.90 0.15 0.10 \\
5 -0.27 0.19 -0.27 1.16 0.04 0.10 \\
6 -0.21 0.15 -0.21 1.19 0.02 0.10 \\
7 0.07 -0.05 0.07 1.24 0.00 0.10 \\
8 -0.25 0.18 -0.25 1.17 0.03 0.10 \\
9 0.15 -0.10 0.15 1.22 0.01 0.10 \\
10 0.05 -0.04 0.05 1.24 0.00 0.10 \\
11 0.00 0.05 0.07 1.24 0.00 0.10 \\
12 0.00 -0.17 -0.24 1.17 0.03 0.10 \\
13 0.00 -0.09 -0.12 1.23 0.01 0.10 \\
14 -0.00 -0.40 -0.57 0.91 0.15 0.10 \\
15 -0.00 0.46 0.66 0.83 0.19 0.10 \\
16 -0.00 -0.30 -0.43 1.04 0.09 0.10 \\
17 -0.00 0.54 0.77 0.72 0.24 0.10 \\
18 -0.00 0.08 0.11 1.23 0.01 0.10 \\
19 0.00 -0.12 -0.17 1.21 0.01 0.10 \\
20 -0.00 0.01 0.01 1.25 0.00 0.10 \\
\hline
\end{tabular}
\end{center}
\end{table}
Might have been better to wait on the xtable processing until you had
added back the column names though.
*
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
***
David Winsemius, MD
West Hartford, CT
__
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 special characters problem
On Jul 22, 2010, at 9:39 AM, Bunny, lautloscrew.com wrote:
Sorry all,
for not posting a minimal example. I am running R on Mac OS X snow
leopard with Komodo edit / Sciviews-R. The problem is that the code
does not work any more if there is an umlaut in my R code. The error
message is some strange mixture of german and english, so this what
the error message is supposed to look like:
Invalid multibyte character in Parser line 195 . Allan´s example
works as a standalone, with his code I get the following error
message. Maybe this is a mac problem ...
Error in parse(text = chunk) :
Unexptected entry in x - data.frame(GeschÂ
I ran Allan Englehart's example and got the expected result on R
2.11.1 with Mac OS 10.5.8.
Seems possible that this is specific to Komodo Edit/Sciviews-R. I
tried getting that editor/system to work as an editing environment a
year ago and formed the opinion that there were too many got-cha's in
the specification of options and user environments. I gave up and went
back to a less complex (and in my hands, less fragile) strategy. I
would take your question to whatever support site is available for
Sciviews-R. (I suppose it could be a Leopard vs Snow Leopard issue,
but do not have the resources to investigate.)
--
David all thumbs with Unix Winsemius
sessionInfo()
R version 2.11.1 Patched (2010-06-14 r52281)
x86_64-apple-darwin9.8.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] grid splines stats graphics grDevices utils
datasets methods
[9] base
other attached packages:
[1] xtable_1.5-6mgcv_1.6-2 prettyR_1.8-1 sos_1.2-9
brew_0.1-1
[6] doBy_4.0.6 MASS_7.3-6 rms_3.0-0 Hmisc_3.8-1
survival_2.35-8
[11] panel_1.0.7 quantreg_4.50 SparseM_0.85lattice_0.18-8
loaded via a namespace (and not attached):
[1] cluster_1.12.3 Matrix_0.999375-40 nlme_3.1-96
tools_2.11.1
Thx for any help in advance
best
matt
On 22.07.2010, at 14:47, Allan Engelhardt wrote:
A standalone example is always helpful. The following works for
me, so I am probably not understanding your problem:
---[BEGIN: umlaut.Rnw]---
\documentclass[a4paper]{article}
\usepackage{ucs}
\usepackage[utf8x]{inputenc}
\begin{document}
test,echo=TRUE=
x - data.frame(Geschäftslage=1:10)
summary(x)
@
\end{document}
---[END: umlaut.Rnw]---
$ R CMD Sweave umlaut.Rnw
$ R CMD pdflatex umlaut.tex
$ gnome-open umlaut.pdf
Allan
On 22/07/10 13:19, Bunny, lautloscrew.com wrote:
Dear all,
I use Sweave to create my reports. Unfortunately my script crashes
whenever I my R code contains special characters like umlauts.
Is there a way to to escape special characters in Sweave... This
is the line that crashes Sweave:
gl_bybranch = ddply(new_wans,.(period,Branchen), function(X)
data.frame(Geschäftslage=mean(X$sentiment)))
Unfortunately I can't just rename it, because I it´s displayed in
the legend of graphics later on.
Thx for any suggestions!
best
matt
__
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.
David Winsemius, MD
West Hartford, CT
__
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] lme4 on Mac OS X
On Jul 22, 2010, at 11:34 AM, Nicholas Griffin wrote: I have been trying to get the lme4 package installed on Mac OS X... with no success. The Mac OS binary is not available on any CRAN, and I also can’t install the package from old source. Has anyone found a solution to this problem? I am happy to use nlme for now, but I tend to prefer to do my mixed model analyses with lmer(). I believe this has been addressed several times recently. The new settings for RSiteSearch no longer search the r-help listings but you can use it to get to Baron's repository and then change the parameters: http://search.r-project.org/cgi-bin/namazu.cgi?query=mac+lme4max=20result=normalsort=scoreidxname=Rhelp10 I have fixed my perception of this design deficiency by creating a new search function in my .Rprofile: rhelpSearch - function(string, restrict = c(Rhelp10, Rhelp08, Rhelp02, functions ), matchesPerPage = 100, ...) RSiteSearch(string=string, restrict = restrict, matchesPerPage = matchesPerPage, ...) David Winsemius, MD West Hartford, CT __ 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] t.test in for loop
On Jul 22, 2010, at 1:28 PM, Arne Schulz wrote:
Dear list,
I'd like to do several t-test within a for loop. Example follows:
data1 - rnorm(1:25)
data2 - rnorm(1:25)
vars - c(data1, data2)
for (i in vars) {
t.test(i)
}
Unfortunately, it does not work because of the quotes in the vars-
vector (running t.test(data1) by hand works).
How can I remove the quotes before performing the test?
a) don't quote them
b) put them in a list
c) print() the results.
--
David Winsemius, MD
West Hartford, CT
__
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] Bar Plot Bars Bleed off Plotting Area
On Jul 22, 2010, at 2:13 PM, Gregory Gilbert wrote: R Community, I have a stupid little barplot I am trying to construct (Windows XP, R.11.1, 32-bit). Whenever I run it my bars run below the horizontal axis. Can anyone (1) reproduce this and (2) offer a solution? I'm not sure I should be offering this, as I consider this practice to be graphical prevarication and statistically unethical, but ... ?barplot read help page so, add , xpd=FALSE, -- David. Rather simplistic code follows (I am new to the community) Thanks for your help! Greg Gilbert === CUT HERE === options(repos=http://lib.stat.cmu.edu/R/CRAN/;) install.packages(plotrix, dependencies=TRUE) update.packages() require(plotrix) ##DATA ##April September ## 1 46.6 50.0 ## 2 43.3 53.3 w - read.csv(C:/Documents and Settings/vhachagilbeg/Desktop/pope/women.csv, header=TRUE) attach(w) names(w) w barplot(as.matrix(w), ylab=Percentage, ylim=c(40, 55), yaxt=n, beside=TRUE, col=(rainbow(2))) text(x=1.5,y=48, 46.6%) text(x=2.5,y=44.3, 43.3%) text(x=4.5,y=51, 50.0%) text(x=5.5,y=54.3, 53.3%) box() axis(2, at=c(40, 45, 50, 55), labels=c(0, 45, 50, 55)) axis.break(2, 41) legend(1.5, 52.5, c(Mammograms, Pap Smears), bty=n, fill=(rainbow(2))) === END === [[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, MD West Hartford, CT __ 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 express bar(zeta) in main title in boxplot
On Jul 22, 2010, at 4:24 PM, Peter Ehlers wrote: On 2010-07-22 11:44, Marcus Liu wrote: Hi everyone, I am plotting a boxplot with main title as main = bquote(paste(.(ts.ind[s]), : , bar(zeta), Boxplot from 2001 to 2009, sep = )) but it doesn't work. The program said they cannot find the function bar. Does anyone know how to do it correctly? Thanks. A reproducible example with the exact error message would be good. Anyway, it seems pretty clear what you want and one solution is to _not_ use 'main='. For base graphics, I usually prefer to add titles with the title() function which will work here. a - pi boxplot(rnorm(200)) title(bquote(paste(.(a), : , bar(zeta), Boxplot from 2001 to 2009, sep = ))) It seems that setting main=... where ... contains bquote() works with plot(), but not with boxplot(). The help page for boxplot does not document a main argument, nor is it in the argument list for boxplot.default or its bxp function. The documentation for the ... argument does not suggest, to me anyway, that main would passed on to other graphical functions. -- David Winsemius, MD West Hartford, CT __ 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 express bar(zeta) in main title in boxplot
On Jul 22, 2010, at 5:01 PM, Peter Ehlers wrote: On 2010-07-22 14:40, David Winsemius wrote: On Jul 22, 2010, at 4:24 PM, Peter Ehlers wrote: On 2010-07-22 11:44, Marcus Liu wrote: Hi everyone, I am plotting a boxplot with main title as main = bquote(paste(.(ts.ind[s]), : , bar(zeta), Boxplot from 2001 to 2009, sep = )) but it doesn't work. The program said they cannot find the function bar. Does anyone know how to do it correctly? Thanks. A reproducible example with the exact error message would be good. Anyway, it seems pretty clear what you want and one solution is to _not_ use 'main='. For base graphics, I usually prefer to add titles with the title() function which will work here. a - pi boxplot(rnorm(200)) title(bquote(paste(.(a), : , bar(zeta), Boxplot from 2001 to 2009, sep = ))) It seems that setting main=... where ... contains bquote() works with plot(), but not with boxplot(). The help page for boxplot does not document a main argument, nor is it in the argument list for boxplot.default or its bxp function. The documentation for the ... argument does not suggest, to me anyway, that main would passed on to other graphical functions. Well, that's not quite so. From ?bxp: and main, cex.main, col.main, sub, cex.sub, col.sub, xlab, ylab, cex.lab, and col.lab are passed to title I stand (or sit) corrected. # To wit: y - rnorm(200) g - gl(2,100) boxplot(y ~ g, main=My title) # But I wasn't right, either: a - pi boxplot(y ~ g, main=paste(hello, bquote(.(a)), goodbye)) boxplot(y ~ g, main=expression(hello~~bar(zeta)~~goodbye)) I do usually prefer to use main= and then supply the title separately. But I'll wager that either Gabor or Uwe or ... will make main= work for the OP. Both are working on my machine. -- David Winsemius, MD West Hartford, CT __ 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] function return
On Jul 22, 2010, at 5:27 PM, Daniel Hocking wrote:
I am sorry if this question is vague or uninformed. I am just
learning R and struggling. I am using the book Hierarchical Modeling
and Inference in Ecology and they provide examples of R code. I have
the following code from the book but when I run it I don't get any
output. I cannot get the values of 'out' to show up. Basically, I
just want to see my estimates for b0, b1, b2, and b3. Any time that I
have used or see function( ) there have been arguments and I just call
for the value of the function. In this case there doesn't seem to be
a value for the function. Any help would be appreciated.
You have defined a function that has no arguments. You should be
typing this at the console:
panel3pt1.fn()
I get an error because there is no such file on my machine and you
have not provided a link to an accessible copy.
--
David.
Thanks,
Dan
`panel3pt1.fn` -
function(){
source(utilfns.Rd)
data- read.table(wtmatrix.txt,header=T,na.strings=T)
elev-data[,elev]
forest-data[,forest]
elev-scale(elev,center=TRUE)
forest-scale(forest,center=TRUE)
pamat-data[,c(y.1,y.2,y.3)]
z-pamat[,1]
M-length(z)
lik-function(parms){
b0-parms[1]
b1-parms[2]
b2-parms[3]
b3-parms[4]
ones-rep(1,M)
### Compute binomial success probabilities
probs-expit(b0*ones+b1*elev+b2*(elev^2)+b3*forest)
lik-rep(0,length(z))
### evaluate log of binomial pmf
tmp-log(dbinom(z,1,probs))
### substitute 0 for missing values
lik[!is.na(z)] - tmp[!is.na(z)]
lik- -1*sum(lik)
return(lik)
}
out - nlm(lik,c(0,0,0,0),hessian=TRUE)
return(out)
}
Daniel J. Hocking
122 James Hall
Department of Natural Resources the Environment
University of New Hampshire
Durham, NH 03824
dhock...@unh.edu
http://sites.google.com/site/danieljhocking/
www.hockingphotography.smugmug.com
Without data, all you are is just another person with an opinion.
-
[[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, MD
West Hartford, CT
__
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] function return
On Jul 22, 2010, at 5:34 PM, David Winsemius wrote:
On Jul 22, 2010, at 5:27 PM, Daniel Hocking wrote:
I am sorry if this question is vague or uninformed. I am just
learning R and struggling. I am using the book Hierarchical Modeling
and Inference in Ecology and they provide examples of R code. I have
the following code from the book but when I run it I don't get any
output. I cannot get the values of 'out' to show up. Basically, I
just want to see my estimates for b0, b1, b2, and b3. Any time
that I
have used or see function( ) there have been arguments and I just
call
for the value of the function. In this case there doesn't seem to be
a value for the function. Any help would be appreciated.
You have defined a function that has no arguments. You should be
typing this at the console:
panel3pt1.fn()
I get an error because there is no such file on my machine and you
have not provided a link to an accessible copy.
If you go to their website
http://www.mbr-pwrc.usgs.gov/pubanalysis/roylebook/chapters.htm
... and:
a) download the .Rd file and the similarly named .csv file referenced
by that function to your working directory
b) change the extension of the .csv file to match the function reference
c) add , sep=,to the read.table call
d) then do what I said above
panel3pt1.fn() # Voilà !!
$minimum
[1] 100.887223663
$estimate
[1] -0.542213281782 1.847343506809 -1.060864239423 0.646910762031
$gradient
[1] 2.34905428442e-05 3.14549972521e-05 6.5196111e-05
4.88853402203e-05
$hessian
[,1] [,2] [,3][,4]
[1,] 33.38181500113 15.26441629493 27.8670270859 6.61193837459
[2,] 15.26441629493 27.86699589010 30.5057286674 -9.66977848175
[3,] 27.86702708591 30.50572866744 52.5892374981 -12.50448633528
[4,] 6.61193837459 -9.66977848175 -12.5044863353 34.09027300450
$code
[1] 1
$iterations
[1] 16
--
David.
Thanks,
Dan
`panel3pt1.fn` -
function(){
source(utilfns.Rd)
data- read.table(wtmatrix.txt,header=T,na.strings=T)
elev-data[,elev]
forest-data[,forest]
elev-scale(elev,center=TRUE)
forest-scale(forest,center=TRUE)
pamat-data[,c(y.1,y.2,y.3)]
z-pamat[,1]
M-length(z)
lik-function(parms){
b0-parms[1]
b1-parms[2]
b2-parms[3]
b3-parms[4]
ones-rep(1,M)
### Compute binomial success probabilities
probs-expit(b0*ones+b1*elev+b2*(elev^2)+b3*forest)
lik-rep(0,length(z))
### evaluate log of binomial pmf
tmp-log(dbinom(z,1,probs))
### substitute 0 for missing values
lik[!is.na(z)] - tmp[!is.na(z)]
lik- -1*sum(lik)
return(lik)
}
out - nlm(lik,c(0,0,0,0),hessian=TRUE)
return(out)
}
Daniel J. Hocking
122 James Hall
Department of Natural Resources the Environment
University of New Hampshire
Durham, NH 03824
dhock...@unh.edu
http://sites.google.com/site/danieljhocking/
www.hockingphotography.smugmug.com
Without data, all you are is just another person with an opinion.
-
[[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, MD
West Hartford, CT
__
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, MD
West Hartford, CT
__
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] Filtering in R
On Jul 22, 2010, at 9:06 PM, stephen sefick wrote: Agian, Plead read the bottom of this email. Also, it looks like you should read An Introduction to R. It certainly does appear this poster needs to do both. On Thu, Jul 22, 2010 at 7:52 PM, jd6688 jdsignat...@gmail.com wrote: The dataframe is id salary 100500 101600 102700 103800 how can i generate a subsets if salary600? ?subset subset(DF, salary 600) -- David Winsemius, MD West Hartford, CT __ 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] Midpoint between coordinates
On Jul 23, 2010, at 7:13 AM, Mafalda Viana wrote: Dear R users, I need to find the coordinates for the point (midpoint) located half way between two pairs of coordinates (lon1,lat1 and lon2,lat2) assuming a straight line between them. What would be the best way? I tried to find an answer in the help archives but without success. I would greatly appreciate any help. df- data.frame(lon1=c(-4.568,-4.3980), lat1=c(59.235,56.369), lon2=c(-5.123,-4.698), lat2=c(60.258,59.197) ) Wouldn't that just be the arithmetic average of the values? Or to you have some need for a more accurate calculation based on spherical geometry? Or some thing that will handle some other coordinate weirdness? df$midlong - apply(df[,c(1,3)], 1, mean) df$midlat - apply(df[,c(2,4)], 1, mean) df lon1 lat1 lon2 lat2 midlong midlat 1 -4.568 59.235 -5.123 60.258 -4.8455 59.7465 2 -4.398 56.369 -4.698 59.197 -4.5480 57.7830 -- David Winsemius, MD West Hartford, CT __ 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] Midpoint between coordinates
On Jul 23, 2010, at 8:58 AM, Mafalda Viana wrote:
The arithmetic mean was my first approach and to nearby points it
doesn't make much difference. However, when the distance between the 2
points gets bigger this is no longer accurate enough. So yes, I was
thinking on spherical geometry, midpoint considering the great circle
distance or similar.
Are you up for some searching?
require(sos)
???distance spherical
found 65 matches; retrieving 4 pages
2 3 4
On first page of hits the geosphere package says it handles spherical
geometry. Perhaps the gcIntermediate function:
gcIntermediate {geosphere}
Get intermediate points on a Great Circle inbetween the two points
used to define the Great Circle
Usage
gcIntermediate(p1, p2, n=50)
Arguments
p1 Longitude/latitude of a single point, in degrees; can be a vector
of two numbers, a matrix of 2 columns (first one is longitude, second
is latitude) or a SpatialPoints* object
p2 As above
n The requested number of points on the Great Circle
---
There is also a geospatial Task View and SIG mailing list that is very
active.
--
David
Thank you
Mafalda
On 23 July 2010 13:30, David Winsemius dwinsem...@comcast.net wrote:
On Jul 23, 2010, at 7:13 AM, Mafalda Viana wrote:
Dear R users,
I need to find the coordinates for the point (midpoint) located half
way between two pairs of coordinates (lon1,lat1 and lon2,lat2)
assuming a straight line between them. What would be the best way? I
tried to find an answer in the help archives but without success. I
would greatly appreciate any help.
df- data.frame(lon1=c(-4.568,-4.3980), lat1=c(59.235,56.369),
lon2=c(-5.123,-4.698), lat2=c(60.258,59.197) )
Wouldn't that just be the arithmetic average of the values? Or to
you have
some need for a more accurate calculation based on spherical
geometry? Or
some thing that will handle some other coordinate weirdness?
df$midlong - apply(df[,c(1,3)], 1, mean)
df$midlat - apply(df[,c(2,4)], 1, mean)
df
lon1 lat1 lon2 lat2 midlong midlat
1 -4.568 59.235 -5.123 60.258 -4.8455 59.7465
2 -4.398 56.369 -4.698 59.197 -4.5480 57.7830
--
David Winsemius, MD
West Hartford, CT
--
Mafalda Viana
Department of Zoology
School of Natural Sciences
Trinity College Dublin
Dublin 2
Ireland
(+353) (0) 872829850
via...@tcd.ie
http://www.tcd.ie/Zoology/research/research/theoretical
David Winsemius, MD
West Hartford, CT
__
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] Figures in Latex
On Jul 23, 2010, at 9:43 AM, li li wrote: Hi all, I want to add 6 plots in the format of 2 columns and 3 rows as one figure in latex. The plots are in .eps file. I know how to add 2 plots side by side, but could not figure out how to do multiple rows. I know this may not be the right place to ask such a question. But I do not know who to ask, http://lmgtfy.com/?q=latex+users+group so just try my luck here. Thank you in advance. Hannah [[alternative HTML version deleted]] -- David Winsemius, MD West Hartford, CT __ 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] randomness using runif
On Jul 23, 2010, at 11:16 AM, Brigid Mooney wrote: I'm working on a problem where I'm introducing random error and have been using the built in function runif to provide that random error. However, I realized that I seem to be getting some unexpected behavior out of the function and was hoping someone could share some insight. I don't know the runif algorithm at all, but from the behavior I'm seeing, it seems that whenever I open a new R console, the function runif gets reset to some initial value. ?set.seed From the help page for set.seed: Note Initially, there is no seed; a new one is created from the current time when one is required. Hence, different sessions started at (sufficiently) different times will give different simulation results, by default. However, the seed might be restored from a previous session if a previously saved workspace is restored. So you probably have saved a value in your workspace. -- David. For example... In a NEW R console, enter the following: x1 - runif(1000, -1, 1) x2 - runif(1000, -1, 1) x1[1:5] x2[1:5] objectsToSave - c(x1, x2) filename - C:\\Documents\\x1x2file.Rdata save(list=objectsToSave, file=filename, envir = parent.frame()) Then in a different NEW R console, enter this: x3 - runif(1000, -1, 1) x4 - runif(1000, -1, 1) x3[1:5] x4[1:5] # For me, the values look identical to x1 and x2, but let's check by loading the x1x2 file and comparing them directly... filename - C:\\Documents\\x1x2file.Rdata load(filename) sum(x1==x3) sum(x2==x4) For my results, I get that x1=x3 for all 1000 elements in the vector, and x2=x4 for all 1000 elements in that vector. Does anyone have insight into what's going on here? Am I doing something wrong here, or is this a quirk of the runif algorithm? Is there a better function out there for seeding truly random error? For what it's worth, here's my R version info: platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major 2 minor 8.1 year 2008 month 12 day22 svn rev47281 language R version.string R version 2.8.1 (2008-12-22) Thanks for the help, Brigid __ 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, MD West Hartford, CT __ 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] converting a time to nearest half-hour
On Jul 23, 2010, at 11:20 AM, murali.me...@avivainvestors.com murali.me...@avivainvestors.com wrote: Hi folks, I've got a POSIXct datum as follows: Sys.time() [1] 2010-07-23 11:29:59 BST I want to convert this to the nearest half-hour, i.e., to 2010-07-23 11:30:00 BST (If the time were 11:59:ss, I want to convert to 12:00:00). How to achieve this? Couldn't you just coerce to numeric, divide by 60(sec)*30(half-hour minutes), round to integer, multiply by 60*30, coerce to POSIXct? _ David Winsemius, MD West Hartford, CT __ 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] decimal seperator
On Jul 23, 2010, at 11:32 AM, Jennifer Sabatier wrote: Hi R-List, I have a question regarding R-language formats, I think. I am producing a series of graphs (using plot, barplot, barchart, and bwplot, using either text or mtext to place values on the graphs) and tables for a Francophone country. In fact, I have already done so. However, while they are pleased with the results they've requested I convert all of my decimal points into the French format which uses commas rather than points. ?options OutDec: character string containing a single-byte character. The character to be used as the decimal point in output conversions, that is in printing, plotting and as.character but not deparsing. -- DAvid. I have no idea how to do this for the graphs, despite searching the help, other than to convert all of my statistics into strings and manually reset the decimal separators, which would take a long time. Is there some quick and easy way? It's only the graphs I need assistance with. For tables I simply output to excel and it's easy to change them there. Here's an example (I'm sure it's crude but I'm still new at R. To make the graph look right you have to expand the java window...which I'm sure you don't need to do if you know how to do this in a more elegant manner): sites - c(Kayes, Kita, Koulikoro, Fana, Sikasso, Koutiala, SgFam, SgHop, Bla, Mopti, Douentz, Tombc, Dire ,Gao ,Ansongo,Kidal,Tessalit,BkoCommI,BkoCommIII,BkoCommV) size - list (2.91,2.36,5.09,3.21,2.27,4.09,2.31,2.76,1.2,2.03,3.06,0.53,1.43,1.83,1,0.93,0,4.01,4.13,3.47 ) site_size - data.frame(cbind(sites, size)) newdata - (sapply(subset(site_size, select=c(size)), as.numeric)) rownames(newdata) - site_size$sites library(grid) plot(newdata, ylab = , xlab= , axes = FALSE)#, type=h, lwd=16) points(newdata, cex = 10, col = topo.colors(20), bg=topo.colors(20), pch=22) lines(newdata, type=h, lwd=40, col=topo.colors(20)) axis(1, at=seq(1, 20, by=1), labels = FALSE) text(seq(1, 20, by=1), par(usr)[3] - 0.2, labels = site2_labels, srt = 45, pos = 1, xpd = TRUE) reg.txt - as.character(c(Kayes Koulikoro Sikasso Segou Mopti Tombouctou Gao Kidal Bamako)) mtext(paste(reg.txt), side=3, font=4, cex=1, adj=0)#, outer=T) text(0, 5.35, Region:, cex = 1, font=4, xpd=T) text(0, -.5, Sites:, cex=1.2, font=1, xpd=T) abline(v=c(2.5, 4.5, 6.5, 9.5, 11.5, 13.5, 15.5, 17.5)) axis(2, at=3, labels = FALSE) mtext(paste(VIH Prevalence (%)), side=2, font=2, cex=1.2) text(1,3, labels=newdata[1], col=white, cex=1.5);text(2,2.45, labels=newdata[2], col=white, cex=1.5) text(3,5.2, labels=newdata[3], col=white, cex=1.5);text(4,3.27, labels=newdata[4], col=white, cex=1.5) text(5,2.35, labels=newdata[5], col=white, cex=1.5);text(6,4.2, labels=newdata[6], col=white, cex=1.5) text(7,2.36, labels=newdata[7], col=black, cex=1.5);text(8,2.85, labels=newdata[8], col=black, cex=1.5) text(9,1.28, labels=newdata[9], col=black, cex=1.5);text(10,2.1, labels=newdata[10], col=black, cex=1.5) text(11,3.13, labels=newdata[11], col=black, cex=1.5);text(12,.6, labels=newdata[12], col=black, cex=1.5) text(13,1.48, labels=newdata[13], col=black, cex=1.5);text(14,1.9, labels=newdata[14], col=black, cex=1.5) text(15,1.05, labels=newdata[15], col=black, cex=1.5);text(16,1, labels=newdata[16], col=black, cex=1.5) text(17,.1, labels=newdata[17], col=black, cex=1.5);text(18,4.1, labels=newdata[18], col=black, cex=1.5) text(19,4.2, labels=newdata[19], col=black, cex=1.5);text(20,3.55, labels=newdata[20], col=black, cex=1.5) David Winsemius, MD West Hartford, CT __ 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] converting a time to nearest half-hour
On Jul 23, 2010, at 11:35 AM, David Winsemius wrote: On Jul 23, 2010, at 11:20 AM, murali.me...@avivainvestors.com murali.me...@avivainvestors.com wrote: Hi folks, I've got a POSIXct datum as follows: Sys.time() [1] 2010-07-23 11:29:59 BST I want to convert this to the nearest half-hour, i.e., to 2010-07-23 11:30:00 BST (If the time were 11:59:ss, I want to convert to 12:00:00). How to achieve this? Couldn't you just coerce to numeric, divide by 60(sec)*30(half-hour minutes), round to integer, multiply by 60*30, coerce to POSIXct? When I tried my method I see that one also needs to add or subtract the proper number of seconds from Universal Time to get the output formatting correct. (Probably demonstrates that I do not have the proper understanding of the right place to employ a TZ specification.). David Winsemius, MD West Hartford, CT __ 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] converting a time to nearest half-hour
On Jul 23, 2010, at 12:04 PM, stephen sefick wrote:
If you have a zoo series this should work. If it doesn't then please
tell me because I think it works.
snap2min - function(zoo, min=00:15:00){
min15 - times(min)
a - aggregate(zoo, trunc(time(zoo), min15), function(x) mean(x,
na.rm=TRUE))
}
This works for producing 10 half-hour intervals of EDT times:
as.POSIXct(60*30*( round( as.numeric( Sys.time()+
60*30*(1:10))/ # the sequence creation
(60*30))) - # divide prior to rounding
5*60*60,# the TZ offset
origin=1970-01-01 )
[1] 2010-07-23 12:30:00 EDT 2010-07-23 13:00:00 EDT
[3] 2010-07-23 13:30:00 EDT 2010-07-23 14:00:00 EDT
[5] 2010-07-23 14:30:00 EDT 2010-07-23 15:00:00 EDT
[7] 2010-07-23 15:30:00 EDT 2010-07-23 16:00:00 EDT
[9] 2010-07-23 16:30:00 EDT 2010-07-23 17:00:00 EDT
hth
Stephen Sefick
On Fri, Jul 23, 2010 at 11:00 AM, David Winsemius
dwinsem...@comcast.net wrote:
On Jul 23, 2010, at 11:35 AM, David Winsemius wrote:
On Jul 23, 2010, at 11:20 AM, murali.me...@avivainvestors.com
murali.me...@avivainvestors.com wrote:
Hi folks,
I've got a POSIXct datum as follows:
Sys.time()
[1] 2010-07-23 11:29:59 BST
I want to convert this to the nearest half-hour, i.e., to
2010-07-23
11:30:00 BST
(If the time were 11:59:ss, I want to convert to 12:00:00).
How to achieve this?
Couldn't you just coerce to numeric, divide by 60(sec)*30(half-hour
minutes), round to integer, multiply by 60*30, coerce to POSIXct?
When I tried my method I see that one also needs to add or subtract
the
proper number of seconds from Universal Time to get the output
formatting
correct. (Probably demonstrates that I do not have the proper
understanding
of the right place to employ a TZ specification.).
David Winsemius, MD
West Hartford, CT
__
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.
--
Stephen Sefick
| Auburn University |
| Department of Biological Sciences |
| 331 Funchess Hall |
| Auburn, Alabama |
| 36849|
|___|
| sas0...@auburn.edu |
| http://www.auburn.edu/~sas0025 |
|___|
Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods. We are mammals, and have not exhausted the
annoying little problems of being mammals.
-K. Mullis
David Winsemius, MD
West Hartford, CT
__
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] Error produced by read.zoo: bad entries
?read.zoo You didn't specify the index column correctly. On Jul 23, 2010, at 12:36 PM, Dimitri Liakhovitski wrote: Hello! I have a data set similar to the data set monthly in the example below: monthly- data .frame (month = c (20090301,20090401,20090501,20100301,20100401,20090301,20090401,20090501,20100301,20100401 ),monthly.value=c(100,200,300,101,201,10,20,30,11,21),market=c(Market A,Market A, Market A,Market A, Market A,Market B, Market B,Market B,Market B, Market B)) monthly$month-as.character(monthly$month) monthly$month-as.Date(monthly$month,%Y%m%d) (monthly) str(monthly) I am trying to use read.zoo - like in 3 lines below: library(zoo) z - read.zoo(monthly, split = market) (z) With the artificially produced data set above, it works just fine. However, with my data it gives me an error: OrigData-read.csv(OrigData.csv) OrigData$Month-as.character(OrigData$Month) OrigData$Month-as.Date(OrigData$Month,%m/%d/%y) str(OrigData) ### The result of str(OrigData) is: 'data.frame': 440 obs. of 3 variables: $ Brand : Factor w/ 11 levels aBrand,bBrand,..: Month :Class 'Date' num [1:440] 13514 13545 13573 13604,... Value: int NA NA NA 100 100 100 100 100 100 99 ?read.zoo You didn't specify the index column correctly. In this case it needs to be = 2. Then I try: z - read.zoo(OrigData, split = Brand) And get the error: Error in read.zoo(OrigData, split = Brand) : index has 440 bad entries at data rows: 1 2 3 4 5 6 7 8 9 10 11 12 13 But the structure of my OrigData is exactly the same as of monthly. OK - OrigData always has a few NAs in Value coming first - but that's consistent for all brands. Any idea what might be wrong? Thanks a lot! Just in case -attaching the actual file. No. Not attached. -- David Winsemius, MD West Hartford, CT __ 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] Error produced by read.zoo: bad entries
But, but, but Did you read my message about the need to correctly specify index columns? The problem is that read.zoo is reading your first column as an index and it's actually the second column that should be used for that purpose. -- David. On Jul 23, 2010, at 1:01 PM, Dimitri Liakhovitski wrote: Strange, I did attach. Attaching again. Maybe the file just doesn't go through? I have: names(OrigData): [1] Brand Month Value I read ?read.zoo According to that index should be the column number. I thought it should be split = 1 in my case - because I am splitting by Brand. But neither split = 1 nor split =2 work. And split =Brand does not work either. Why? D. On Fri, Jul 23, 2010 at 12:52 PM, David Winsemius dwinsem...@comcast.net wrote: ?read.zoo You didn't specify the index column correctly. On Jul 23, 2010, at 12:36 PM, Dimitri Liakhovitski wrote: Hello! I have a data set similar to the data set monthly in the example below: monthly- data .frame (month = c (20090301,20090401,20090501,20100301,20100401,20090301,20090401,20090501,20100301,20100401 ),monthly .value=c(100,200,300,101,201,10,20,30,11,21),market=c(Market A,Market A, Market A,Market A, Market A,Market B, Market B,Market B,Market B, Market B)) monthly$month-as.character(monthly$month) monthly$month-as.Date(monthly$month,%Y%m%d) (monthly) str(monthly) I am trying to use read.zoo - like in 3 lines below: library(zoo) z - read.zoo(monthly, split = market) (z) With the artificially produced data set above, it works just fine. However, with my data it gives me an error: OrigData-read.csv(OrigData.csv) OrigData$Month-as.character(OrigData$Month) OrigData$Month-as.Date(OrigData$Month,%m/%d/%y) str(OrigData) ### The result of str(OrigData) is: 'data.frame': 440 obs. of 3 variables: $ Brand : Factor w/ 11 levels aBrand,bBrand,..: Month :Class 'Date' num [1:440] 13514 13545 13573 13604,... Value: int NA NA NA 100 100 100 100 100 100 99 ?read.zoo You didn't specify the index column correctly. In this case it needs to be = 2. Then I try: z - read.zoo(OrigData, split = Brand) And get the error: Error in read.zoo(OrigData, split = Brand) : index has 440 bad entries at data rows: 1 2 3 4 5 6 7 8 9 10 11 12 13 But the structure of my OrigData is exactly the same as of monthly. OK - OrigData always has a few NAs in Value coming first - but that's consistent for all brands. Any idea what might be wrong? Thanks a lot! Just in case -attaching the actual file. No. Not attached. -- David Winsemius, MD West Hartford, CT -- Dimitri Liakhovitski Ninah Consulting www.ninah.com OrigData.csv David Winsemius, MD West Hartford, CT __ 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] Error produced by read.zoo: bad entries
On Jul 23, 2010, at 1:39 PM, Dimitri Liakhovitski wrote: Very sorry - I mistunderstood and confused split with index.column - totally my fault. Ok, now I've run this line: z - read.zoo(OrigData, index.column = 2, split = Brand) And I am getting: Error in merge.zoo(` Plus` = c(NA, 98L, 95L, 97L, NA, 98L, 97L, 98L, NA, : series cannot be merged with non-unique index entries in a series In addition: There were 11 warnings (use warnings() to see them) I got the warnings but no error: z - read.zoo(OrigData, split = Brand, index.column=2) There were 11 warnings (use warnings() to see them) z Plus agrow chool gress Grib inKid kid omis plet pro romil [1,]NANANANA NANA NA NA NA NANA [2,]98999897 9696 100 97 97 9996 [3,]95 1009799 9297 100 97 99 10099 [4,]97999498 9195 99 98 98 9995 [5,]NANANANA NANA NA NA NA NANA [6,]98999897 9396 99 97 98 9996 [7,]97 1009898 9596 99 98 98 10097 [8,]98999499 9696 99 98 98 9997 [9,]NANANANA NANA NA NA NA NANA [10,]98999898 9596 99 98 98 9997 [11,]98999899 9796 99 98 97 9999 [12,]97 1009699 9595 99 99 97 10096 [13,]96 1009696 93 0 100 96 97 10096 [14,]989998 100 9496 100 98 97 9999 [15,]95 1009899 9395 99 99 99 9999 [16,]97999699 9495 98 98 90 9995 [17,]97 1009796 92 0 100 96 98 10095 [18,]96999898 9697 100 98 99 9898 [19,]98 1009898 9697 99 98 99 9998 [20,]98 1009796 95 0 100 96 98 9996 [21,]94 1009899 9297 99 98 98 9898 [22,]98999897 9696 99 97 98 9997 [23,]97 1009696 93 0 100 95 97 10095 [24,]97 1009897 9396 99 97 98 9795 [25,]98 1009697 9694 100 97 99 9996 [26,]98 1009896 95 0 100 96 98 9995 [27,]98 1009897 9396 96 97 98 9999 [28,]99 1009898 9296 100 98 99 9997 [29,]98 1009795 95 0 100 95 98 9995 [30,]99 10098 100 9898 99 100 99 10099 [31,]97999497 9595 99 97 98 9894 [32,]98999896 95 3 100 96 97 9996 [33,]97999899 9797 99 99 99 9999 [34,]96999596 9494 98 96 96 9893 [35,]98999897 9454 100 97 97 9996 [36,]95 1009799 9595 99 99 98 10099 [37,]98999898 9596 99 98 99 9997 [38,]98999897 9694 100 97 97 9896 [39,]95 10098 100 9597 100 99 99 10099 [40,]97 1009598 9396 99 98 98 9996 Since you didn't say what was expected, I am not in a position to know if this is success. And under warnings() it says: 1: In zoo(rval4[[i]], ix[[i]]) : some methods for “zoo” objects do not work if the index entries in ‘order.by’ are not unique On Fri, Jul 23, 2010 at 1:13 PM, David Winsemius dwinsem...@comcast.net wrote: But, but, but Did you read my message about the need to correctly specify index columns? The problem is that read.zoo is reading your first column as an index and it's actually the second column that should be used for that purpose. -- David. On Jul 23, 2010, at 1:01 PM, Dimitri Liakhovitski wrote: Strange, I did attach. Attaching again. Maybe the file just doesn't go through? I have: names(OrigData): [1] Brand Month Value I read ?read.zoo According to that index should be the column number. I thought it should be split = 1 in my case - because I am splitting by Brand. But neither split = 1 nor split =2 work. And split =Brand does not work either. Why? D. On Fri, Jul 23, 2010 at 12:52 PM, David Winsemius dwinsem...@comcast.net wrote: ?read.zoo You didn't specify the index column correctly. On Jul 23, 2010, at 12:36 PM, Dimitri Liakhovitski wrote: Hello! I have a data set similar to the data set monthly in the example below: monthly- data .frame (month = c (20090301,20090401,20090501,20100301,20100401,20090301,20090401,20090501,20100301,20100401 ),monthly .value=c(100,200,300,101,201,10,20,30,11,21),market=c(Market A,Market A, Market A,Market A, Market A,Market B, Market B,Market B,Market B, Market B)) monthly$month-as.character(monthly$month) monthly$month-as.Date(monthly$month,%Y%m%d) (monthly) str(monthly) I am trying
Re: [R] Error produced by read.zoo: bad entries
On Jul 23, 2010, at 1:50 PM, Dimitri Liakhovitski wrote: I am expecting to see the week names as row labels of z and the corresponding values (like in the monthly example). I am pretty sure - in order to get it one needs to install the latest version of zoo. I've done it just a couple of days ago. I am getting the error - and nothing is produced. Can it have to do with the fact that I am using the newer version of zoo? Again, my full code for that OrigData.csv file I sent is: Yep, updating to the current version of zoo on CRAN, zoo_1.6-4, now produces an error where before with the penultimate version, zoo_1.6-3, it did not. -- David. OrigData-read.csv(OrigData.csv) OrigData$Month-as.character(OrigData$Month) OrigData$Month-as.Date(OrigData$Month,%m/%d/%y) str(OrigData) 'data.frame': 440 obs. of 3 variables: $ Brand: Factor w/ 11 levels Plus,agrow,..: 2 2 2 2 2 2 2 2 2 2 ... $ Month:Class 'Date' num [1:440] 18262 18293 18322 18353 18383 ... $ Value: int NA NA NA 100 100 100 100 100 100 99 ... library(zoo) z - read.zoo(OrigData, index.column = 2, split = Brand) Error in merge.zoo(` Plus` = c(NA, 98L, 95L, 97L, NA, 98L, 97L, 98L, NA, : series cannot be merged with non-unique index entries in a series In addition: There were 11 warnings (use warnings() to see them) warnings() Warning messages: 1: In zoo(rval4[[i]], ix[[i]]) : some methods for “zoo” objects do not work if the index entries in ‘order.by’ are not unique 2: In zoo(rval4[[i]], ix[[i]]) : some methods for “zoo” objects do not work if the index entries in ‘order.by’ are not unique 3: In zoo(rval4[[i]], ix[[i]]) : etc. But it does not give me this error for my Monthly example - even when I introduce a few NAs there. And I get this message: Error in merge.zoo(` Plus` = c(NA, 98L, 95L, 97L, NA, 98L, 97L, 98L, NA, : series cannot be merged with non-unique index entries in a series In addition: There were 11 warnings (use warnings() to see them) On Fri, Jul 23, 2010 at 1:41 PM, David Winsemius dwinsem...@comcast.net wrote: On Jul 23, 2010, at 1:39 PM, Dimitri Liakhovitski wrote: Very sorry - I mistunderstood and confused split with index.column - totally my fault. Ok, now I've run this line: z - read.zoo(OrigData, index.column = 2, split = Brand) And I am getting: Error in merge.zoo(` Plus` = c(NA, 98L, 95L, 97L, NA, 98L, 97L, 98L, NA, : series cannot be merged with non-unique index entries in a series In addition: There were 11 warnings (use warnings() to see them) I got the warnings but no error: z - read.zoo(OrigData, split = Brand, index.column=2) There were 11 warnings (use warnings() to see them) z Plus agrow chool gress Grib inKid kid omis plet pro romil [1,]NANANANA NANA NA NA NA NANA [2,]98999897 9696 100 97 97 9996 [3,]95 1009799 9297 100 97 99 10099 [4,]97999498 9195 99 98 98 9995 [5,]NANANANA NANA NA NA NA NANA [6,]98999897 9396 99 97 98 9996 [7,]97 1009898 9596 99 98 98 10097 [8,]98999499 9696 99 98 98 9997 [9,]NANANANA NANA NA NA NA NANA [10,]98999898 9596 99 98 98 9997 [11,]98999899 9796 99 98 97 9999 [12,]97 1009699 9595 99 99 97 10096 [13,]96 1009696 93 0 100 96 97 10096 [14,]989998 100 9496 100 98 97 9999 [15,]95 1009899 9395 99 99 99 9999 [16,]97999699 9495 98 98 90 9995 [17,]97 1009796 92 0 100 96 98 10095 [18,]96999898 9697 100 98 99 9898 [19,]98 1009898 9697 99 98 99 9998 [20,]98 1009796 95 0 100 96 98 9996 [21,]94 1009899 9297 99 98 98 9898 [22,]98999897 9696 99 97 98 9997 [23,]97 1009696 93 0 100 95 97 10095 [24,]97 1009897 9396 99 97 98 9795 [25,]98 1009697 9694 100 97 99 9996 [26,]98 1009896 95 0 100 96 98 9995 [27,]98 1009897 9396 96 97 98 9999 [28,]99 1009898 9296 100 98 99 9997 [29,]98 1009795 95 0 100 95 98 9995 [30,]99 10098 100 9898 99 100 99 10099 [31,]97999497 9595 99 97 98 9894 [32,]98999896 95 3 100 96 97 9996 [33,]97999899 9797 99 99 99 9999 [34,]96999596 9494 98 96 96 9893 [35,]98999897 9454 100 97 97
Re: [R] Convert Row Names to data.frame column
On Jul 23, 2010, at 6:18 PM, chipmaney wrote: Here is an example dataset: ZoneCover.df- data.frame(Value=c(1,2)) row.names(ZoneCover.df) - c(Floodplain1.Tree, Floodplain1.Shrub) I want to Export the Row.Names to a column in the dataframe: ZoneCover.df$ID - names(ZoneCover.df) which yields this: ZoneCover.df ValueID Floodplain1.Tree 1 Floodplain1.Tree Floodplain1.Shrub 2 Floodplain1.Shrub QUESTION: How do I remove the .Tree and .Shrub extensions from the ZoneCover$ID values? ?gsub ?regex (The . character needs to be escaped with doubled \ but the second . is a regex for any character.) gsub(\\..*$,, c(Floodplain1.Tree, Floodplain1.Shrub)) [1] Floodplain1 Floodplain1 Use the same method with assignment to the column: ZoneCover.df$ID - gsub(\\..*$,, ZoneCover.df$ID) -- David Winsemius, MD West Hartford, CT __ 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 import simple java/mathematica expression to R
On Jul 23, 2010, at 10:29 PM, jim holtman wrote:
Well, I took you equation and put the following at the start:
Power - function(x,y) x^y
EE - function(x) x
alp - 2
x - 900*Power(-0.2030178326474623 + 0.23024073983368956*(1 -
alp) +
0.2807352820970084*(1 - alp)*(1 - alp*(1 + EE(1))) +
0.2145643524071315*(1 - alp)*
Power(1 - alp*(1 + EE(1)),2) + 0.11519022530097237*(1 -
alp)*Power(1 - alp*(1 + EE(1)),3) +
0.046127977611990736*(1 - alp)*Power(1 - alp*(1 + EE(1)),4) +
0.014279410543117517*(1 - alp)*
Power(1 - alp*(1 + EE(1)),5) + 0.2145643524071315*(Power(1 -
alp,2)*alp*(1 + EE(1)) +
(1 - alp)*alp*(1 + EE(1))*(1 - alp*(1 + EE(1))) + Power(1 -
alp,2)*alp*(1 + EE(2))) +
0.11519022530097237*(Power(1 - alp,2)*alp*(1 + EE(1))*(1 - alp*(1
+ EE(1))) +
2*(1 - alp)*alp*(1 + EE(1))*Power(1 - alp*(1 + EE(1)),2) +
Power(1 - alp,2)*alp*(1 - alp*(1 + EE(1)))*(1 + EE(2))) +
.
since there appeared to be functions Power and EE, and the
variable 'alp'. This evaluated as-is to: 32157617213
so what you have appears to be a legal equation that R can parse as
is. There you can probably put it in a function and use one of the R
routines to minimize it, It does not look like you will have any
problem importing it to R; just have to make sure you have the
appropriate functions defined.
It has some interesting properties with that identity definition of
EE(x). Local maxima at 1 and 0, blows up beyond -1 and 2, and
several local minima nearby:
Math.Fn - function(alp) { big-long-expression }
plot( seq(-.3,1.5,by=0.01), Math.Fn(seq(-.3,1.5,by=0.01) ), cex=0.2)
--
David.
On Fri, Jul 23, 2010 at 1:29 PM, Andrey Siver
andrey.si...@gmail.com wrote:
Hello,
2010/7/23 jim holtman jholt...@gmail.com:
It would be nice if you could post what the data looks like that you
want to import. R can import any text file and then you have string
manipulation that you can do to parse it. So the basic answer is
probably yes, but we do need to understand the format of the data to
give a more precise answer.
I put the target expression to minimize (with some constrains) here:
http://analytic-products4you.com/target.txt
Is it possible to import it as a function to minimize?
What is the problem that you are trying to solve?
We solve a problem for parameters estimation with ties.
--
David Winsemius, MD
West Hartford, CT
__
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] Trouble retrieving the second largest value from each row of a data.frame
On Jul 23, 2010, at 9:20 PM, mpw...@illinois.edu wrote: I have a data frame with a couple million lines and want to retrieve the largest and second largest values in each row, along with the label of the column these values are in. For example row 1 strongest=-11072 secondstrongest=-11707 strongestantenna=value120 secondstrongantenna=value60 Below is the code I am using and a truncated data.frame. Retrieving the largest value was easy, but I have been getting errors every way I have tried to retrieve the second largest value. I have not even tried to retrieve the labels for the value yet. Any help would be appreciated Mike Using Holtman's extract of your data with x as the name and the order function to generate an index to names of the dataframe: t(apply(x, 1, sort, decreasing=TRUE)[1:3, ]) [,1] [,2] [,3] 1 -11072 -11707 -12471 2 -11176 -11799 -12838 3 -3 -11778 -12439 4 -11071 -11561 -11653 5 -11067 -11638 -12834 6 -11068 -11698 -12430 7 -11092 -11607 -11709 8 -11061 -11426 -11665 9 -11137 -11736 -12570 10 -11146 -11779 -12537 Putting it all together: matrix( paste( t(apply(x, 1, sort, decreasing=TRUE)[1:3, ]), names(x)[ t(apply(x, 1, order, decreasing=TRUE) [1:3, ]) ]), ncol=3) [,1] [,2] [,3] [1,] -11072 value120 -11707 value60 -12471 value180 [2,] -11176 value120 -11799 value180 -12838 value0 [3,] -3 value120 -11778 value60 -12439 value180 [4,] -11071 value120 -11561 value240 -11653 value60 [5,] -11067 value120 -11638 value180 -12834 value0 [6,] -11068 value0 -11698 value60 -12430 value120 [7,] -11092 value120 -11607 value240 -11709 value180 [8,] -11061 value120 -11426 value240 -11665 value60 [9,] -11137 value120 -11736 value60 -12570 value180 [10,] -11146 value300 -11779 value0 -12537 value180 -- David. data-data.frame(value0,value60,value120,value180,value240,value300) data value0 value60 value120 value180 value240 value300 1 -13007 -11707 -11072 -12471 -12838 -13357 2 -12838 -13210 -11176 -11799 -13210 -13845 3 -12880 -11778 -3 -12439 -13089 -13880 4 -12805 -11653 -11071 -12385 -11561 -13317 5 -12834 -13527 -11067 -11638 -13527 -13873 6 -11068 -11698 -12430 -12430 -12430 -12814 7 -12807 -14068 -11092 -11709 -11607 -13025 8 -12770 -11665 -11061 -12373 -11426 -12805 9 -12988 -11736 -11137 -12570 -13467 -13739 10 -11779 -12873 -12973 -12537 -12973 -11146 #largest value in the row strongest-apply(data,1,max) #second largest value in the row n-function(data)(1/(min(1/(data[1,]-max(data[1,]+ (max(data[1,]))) secondstrongest-apply(data,1,n) Error in data[1, ] : incorrect number of dimensions __ 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] Trouble retrieving the second largest value from each row of a data.frame
On Jul 24, 2010, at 4:54 PM, mpw...@illinois.edu wrote:
THANKS, but I have one issue and one question.
For some reason the secondstrongest value for row 3 and 6 are
incorrect (they are the strongest) the remaining 10 are correct??
In my run of Wiley's code I instead get identical values for rows
2,5,6. Holtman's and my solutions did not suffer from that defect,
although mine suffered from my misreading of your request, thinking
that you wanted the top 3. The fix is trivial
These data are being used to track radio-tagged birds, they are from
automated radio telemetry receivers. I will applying the following
formula
diff - ((strongest- secondstrongest)/100)
bearingdiff -30-(-0.0624*(diff**2))-(2.8346*diff)
vals - c(value0, value60, value120, value180, value240,
value300)
value.str2 - (match(yourdata$secondstrongestantenna, vals)-1)*60
value.str1 - (match(yourdata$strongestantenna, vals)-1)*60
change.ind - abs(match(yourdata, vals) - which(match(yourdata, vals) )
A) Then the bearing diff is added to strongestantenna (value0 =
0degrees) if the secondstrongestatenna is greater (eg value0 and
value60),
B) or if the secondstrongestantenna is smaller than the
strongestantenna,
then the bearingdiff is substracted from the strongestantenna.
C) The only exception is that if value0 (0degrees) is strongest and
value300(360degrees) is the secondstrongestantenna then the bearing
is 360-bearingdiff.
D) Also the strongestantenna and secondstrongestantenna have to be
next to each other (e.g. value0 with value60, value240 with
value300, value0 with value300) or the results should be NA.
After setting finalbearing with A, B, and C then:
yourdata$finalbearing - with(yourdata, ifelse (
change.ind 5 change.ind 1 ,
NA, finalbearing) )
I have been trying to use a series of if,else statements to produce
these bearing, but all I am producing is errors. Any suggestion
would be appreciated.
Again THANKS for you efforts.
Mike
Original message
Date: Fri, 23 Jul 2010 23:01:56 -0700
From: Joshua Wiley jwiley.ps...@gmail.com
Subject: Re: [R] Trouble retrieving the second largest value from
each row of a data.frame
To: mpw...@illinois.edu
Cc: r-help@r-project.org
Hi,
Here is a little function that will do what you want and return a
nice output:
#Function To calculate top two values and return
my.finder - function(mydata) {
my.fun - function(data) {
strongest - which.max(data)
secondstrongest - which.max(data[-strongest])
strongestantenna - names(data)[strongest]
secondstrongantenna - names(data[-strongest])[secondstrongest]
value - matrix(c(data[strongest], data[secondstrongest],
strongestantenna, secondstrongantenna), ncol =4)
return(value)
}
dat - apply(mydata, 1, my.fun)
dat - t(dat)
dat - as.data.frame(dat, stringsAsFactors = FALSE)
colnames(dat) - c(strongest, secondstrongest,
strongestantenna, secondstrongantenna)
dat[ , strongest] - as.numeric(dat[ , strongest])
dat[ , secondstrongest] - as.numeric(dat[ , secondstrongest])
return(dat)
}
#Using your example data:
yourdata - structure(list(value0 = c(-13007L, -12838L, -12880L,
-12805L,
-12834L, -11068L, -12807L, -12770L, -12988L, -11779L), value60 =
c(-11707L,
-13210L, -11778L, -11653L, -13527L, -11698L, -14068L, -11665L,
-11736L, -12873L), value120 = c(-11072L, -11176L, -3L, -11071L,
-11067L, -12430L, -11092L, -11061L, -11137L, -12973L), value180 =
c(-12471L,
-11799L, -12439L, -12385L, -11638L, -12430L, -11709L, -12373L,
-12570L, -12537L), value240 = c(-12838L, -13210L, -13089L, -11561L,
-13527L, -12430L, -11607L, -11426L, -13467L, -12973L), value300 =
c(-13357L,
-13845L, -13880L, -13317L, -13873L, -12814L, -13025L, -12805L,
-13739L, -11146L)), .Names = c(value0, value60, value120,
value180, value240, value300), class = data.frame,
row.names = c(1,
2, 3, 4, 5, 6, 7, 8, 9, 10))
my.finder(yourdata) #and what you want is in a nicely labeled data
frame
#A potential problem is that it is not very efficient
#Here is a test using a matrix of 100,000 rows
#sampled from the same range as your data
#with the same number of columns
data.test - matrix(
sample(seq(min(yourdata),max(yourdata)), size = 50, replace =
TRUE),
ncol = 5)
system.time(my.finder(data.test))
#On my system I get
system.time(my.finder(data.test))
user system elapsed
2.890.002.89
Hope that helps,
Josh
On Fri, Jul 23, 2010 at 6:20 PM, mpw...@illinois.edu wrote:
I have a data frame with a couple million lines and want to
retrieve the largest and second largest values in each row, along
with the label of the column these values are in. For example
row 1
strongest=-11072
secondstrongest=-11707
strongestantenna=value120
secondstrongantenna=value60
Below is the code I am using and a truncated data.frame.
Retrieving the largest value was easy, but
Re: [R] Trouble retrieving the second largest value from each row of a data.frame
On Jul 24, 2010, at 8:09 PM, David Winsemius wrote:
On Jul 24, 2010, at 4:54 PM, mpw...@illinois.edu wrote:
THANKS, but I have one issue and one question.
For some reason the secondstrongest value for row 3 and 6 are
incorrect (they are the strongest) the remaining 10 are correct??
In my run of Wiley's code I instead get identical values for rows
2,5,6. Holtman's and my solutions did not suffer from that defect,
although mine suffered from my misreading of your request, thinking
that you wanted the top 3. The fix is trivial
These data are being used to track radio-tagged birds, they are
from automated radio telemetry receivers. I will applying the
following formula
diff - ((strongest- secondstrongest)/100)
bearingdiff -30-(-0.0624*(diff**2))-(2.8346*diff)
vals - c(value0, value60, value120, value180, value240,
value300)
value.str2 - (match(yourdata$secondstrongestantenna, vals)-1)*60
value.str1 - (match(yourdata$strongestantenna, vals)-1)*60
change.ind - abs(match(yourdata, vals) - which(match(yourdata,
vals) )
OOOPs should have been
change.ind - abs(match(yourdata, vals) - match(yourdata, vals) )
A) Then the bearing diff is added to strongestantenna (value0 =
0degrees) if the secondstrongestatenna is greater (eg value0 and
value60),
B) or if the secondstrongestantenna is smaller than the
strongestantenna,
then the bearingdiff is substracted from the strongestantenna.
yourdata$finalbearing - with(yourdata, ifelse (value.str2value.str1,
bearingdiff+value.str1, value.str1-bearingdiff) )
C) The only exception is that if value0 (0degrees) is strongest and
value300(360degrees) is the secondstrongestantenna then the bearing
is 360-bearingdiff.
yourdata$finalbearing - with(yourdata, ifelse (strongestantenna ==
value0 secondstrongantenna == value300, 360- bearingdiff,
finalbearing) );
D) Also the strongestantenna and secondstrongestantenna have to be
next to each other (e.g. value0 with value60, value240 with
value300, value0 with value300) or the results should be NA.
After setting finalbearing with A, B, and C then:
yourdata$finalbearing - with(yourdata, ifelse (
change.ind 5 change.ind 1 ,
NA, finalbearing) )
yourdata
strongest secondstrongest strongestantenna secondstrongantenna
finalbearing
1 -11072 -11707 value120 value60
-11086.52
2 -11176 -11799 value120value180
-11190.76
3 -3 -11778 value120 value60
-11126.91
4 -11071 -11561 value120
value240 NA
5 -11067 -11638 value120value180
-11082.85
6 -11068 -11698 value0 value60
-11082.62
7 -11092 -11607 value120
value240 NA
8 -11061 -11426 value120
value240 NA
9 -11137 -11736 value120 value60
-11152.26
10-11146 -11779 value300 value0
-11160.56
I have been trying to use a series of if,else statements to produce
these bearing,
ifelse is the correct construct for processing vectors
--
David.
but all I am producing is errors. Any suggestion would be
appreciated.
Again THANKS for you efforts.
Mike
Original message
Date: Fri, 23 Jul 2010 23:01:56 -0700
From: Joshua Wiley jwiley.ps...@gmail.com
Subject: Re: [R] Trouble retrieving the second largest value from
each row of a data.frame
To: mpw...@illinois.edu
Cc: r-help@r-project.org
Hi,
Here is a little function that will do what you want and return a
nice output:
#Function To calculate top two values and return
my.finder - function(mydata) {
my.fun - function(data) {
strongest - which.max(data)
secondstrongest - which.max(data[-strongest])
strongestantenna - names(data)[strongest]
secondstrongantenna - names(data[-strongest])[secondstrongest]
value - matrix(c(data[strongest], data[secondstrongest],
strongestantenna, secondstrongantenna), ncol =4)
return(value)
}
dat - apply(mydata, 1, my.fun)
dat - t(dat)
dat - as.data.frame(dat, stringsAsFactors = FALSE)
colnames(dat) - c(strongest, secondstrongest,
strongestantenna, secondstrongantenna)
dat[ , strongest] - as.numeric(dat[ , strongest])
dat[ , secondstrongest] - as.numeric(dat[ , secondstrongest])
return(dat)
}
#Using your example data:
yourdata - structure(list(value0 = c(-13007L, -12838L, -12880L,
-12805L,
-12834L, -11068L, -12807L, -12770L, -12988L, -11779L), value60 =
c(-11707L,
-13210L, -11778L, -11653L, -13527L, -11698L, -14068L, -11665L,
-11736L, -12873L), value120 = c(-11072L, -11176L, -3L, -11071L,
-11067L, -12430L, -11092L, -11061L, -11137L, -12973L), value180 =
c(-12471L,
-11799L, -12439L
Re: [R] Trouble retrieving the second largest value from each row of a data.frame
On Jul 24, 2010, at 9:27 PM, David Winsemius wrote:
On Jul 24, 2010, at 8:09 PM, David Winsemius wrote:
On Jul 24, 2010, at 4:54 PM, mpw...@illinois.edu wrote:
THANKS, but I have one issue and one question.
For some reason the secondstrongest value for row 3 and 6 are
incorrect (they are the strongest) the remaining 10 are correct??
In my run of Wiley's code I instead get identical values for rows
2,5,6. Holtman's and my solutions did not suffer from that defect,
although mine suffered from my misreading of your request, thinking
that you wanted the top 3. The fix is trivial
These data are being used to track radio-tagged birds, they are
from automated radio telemetry receivers. I will applying the
following formula
diff - ((strongest- secondstrongest)/100)
bearingdiff -30-(-0.0624*(diff**2))-(2.8346*diff)
vals - c(value0, value60, value120, value180, value240,
value300)
value.str2 - (match(yourdata$secondstrongestantenna, vals)-1)*60
Had a misspelling ... rather:
match(yourdata$secondstrongantenna, vals)
value.str1 - (match(yourdata$strongestantenna, vals)-1)*60
change.ind - abs(match(yourdata, vals) - which(match(yourdata,
vals) )
OOOPs should have been
change.ind - abs(match(yourdata, vals) - match(yourdata, vals) )
A) Then the bearing diff is added to strongestantenna (value0 =
0degrees) if the secondstrongestatenna is greater (eg value0 and
value60),
B) or if the secondstrongestantenna is smaller than the
strongestantenna,
then the bearingdiff is substracted from the strongestantenna.
yourdata$finalbearing - with(yourdata, ifelse
(value.str2value.str1, bearingdiff+value.str1, value.str1-
bearingdiff) )
C) The only exception is that if value0 (0degrees) is strongest
and value300(360degrees) is the secondstrongestantenna then the
bearing is 360-bearingdiff.
yourdata$finalbearing - with(yourdata, ifelse (strongestantenna ==
value0 secondstrongantenna == value300, 360- bearingdiff,
finalbearing) );
D) Also the strongestantenna and secondstrongestantenna have to be
next to each other (e.g. value0 with value60, value240 with
value300, value0 with value300) or the results should be NA.
After setting finalbearing with A, B, and C then:
yourdata$finalbearing - with(yourdata, ifelse (
change.ind 5 change.ind 1 ,
NA, finalbearing) )
Better result with proper creation of value.str2:
yourdata
strongest secondstrongest strongestantenna secondstrongantenna
finalbearing
1 -11072 -11707 value120 value60
105.48359
2 -11176 -11799 value120value180
134.76237
3 -3 -11778 value120 value60
106.09061
4 -11071 -11561 value120
value240 NA
5 -11067 -11638 value120value180
135.84893
6 -11068 -11698 value0 value60
14.61868
7 -11092 -11607 value120
value240 NA
8 -11061 -11426 value120
value240 NA
9 -11137 -11736 value120 value60
104.74034
10-11146 -11779 value300 value0
285.44272
I have been trying to use a series of if,else statements to
produce these bearing,
ifelse is the correct construct for processing vectors
--
David.
but all I am producing is errors. Any suggestion would be
appreciated.
Again THANKS for you efforts.
Mike
Original message
Date: Fri, 23 Jul 2010 23:01:56 -0700
From: Joshua Wiley jwiley.ps...@gmail.com
Subject: Re: [R] Trouble retrieving the second largest value from
each row of a data.frame
To: mpw...@illinois.edu
Cc: r-help@r-project.org
Hi,
Here is a little function that will do what you want and return a
nice output:
#Function To calculate top two values and return
my.finder - function(mydata) {
my.fun - function(data) {
strongest - which.max(data)
secondstrongest - which.max(data[-strongest])
strongestantenna - names(data)[strongest]
secondstrongantenna - names(data[-strongest])[secondstrongest]
value - matrix(c(data[strongest], data[secondstrongest],
strongestantenna, secondstrongantenna), ncol =4)
return(value)
}
dat - apply(mydata, 1, my.fun)
dat - t(dat)
dat - as.data.frame(dat, stringsAsFactors = FALSE)
colnames(dat) - c(strongest, secondstrongest,
strongestantenna, secondstrongantenna)
dat[ , strongest] - as.numeric(dat[ , strongest])
dat[ , secondstrongest] - as.numeric(dat[ , secondstrongest])
return(dat)
}
#Using your example data:
yourdata - structure(list(value0 = c(-13007L, -12838L, -12880L,
-12805L,
-12834L, -11068L, -12807L, -12770L, -12988L, -11779L), value60 =
c(-11707L,
-13210L, -11778L, -11653L, -13527L, -11698L, -14068L, -11665L
Re: [R] R equivalent of SAS proc freq
On Jul 25, 2010, at 9:32 AM, Sébastien Bihorel wrote: Dear R-users, I am looking for a R function that would be the equivalent of the SAS proc freq ( http://support.sas.com/documentation/cdl/en/procstat/63104/HTML/default/viewer.htm#/ documentation/cdl/en/procstat/63104/HTML/default/ procstat_freq_sect006.htm). The table, ftable, xtabs functions are close but do not quite offer the same capabilities (e.g. they just return counts and no %ages as far as I could tell). I wanted to check with the group if there would already be such a function in base or in a contribution package, before I start coding my own function. There are a couple of efforts to imitate the Proc FREQ output in the archives. If you use one of the several search facilities (e.g. RSeek, the RSiteSearch function or the ??? function in the sos library) you will get links to those: (A test of this theory with RSeek pointed to the CrossTable function in gmodels.) Try also with... ?RsiteSearch RSiteSearch() Or with ... require(sos) ???proc freq sas -- David Winsemius, MD West Hartford, CT __ 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] 3d topographic map
On Jul 25, 2010, at 6:31 PM, sh...@ucar.edu wrote: Hi All- I would like to create a 3d topographic map using lat/lon and z(height). I have been scouring the R help pages and have not located the package I am looking for. Does anyone have a suggestion of package that will work for this? thanks- I suspect most viewers of this message are going to be puzzled. If wireframe and contourplot are not doing it for you, what are the problems? Surely you found references to Lattice and wireframe. The Lattice system has a worked example in Figure 6.11 from this page: http://lmdvr.r-forge.r-project.org/figures/figures.html A search for topographic in one of the r searching facilities is sure to bring multiple hits, . including the usual starting point: https://svn.r-project.org/R/trunk/src/library/graphics/demo/graphics.R So there must be details and issues that you have not chosen to disclose. -- 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.
Re: [R] Marginal effects from interaction regression model
On Jul 25, 2010, at 10:24 PM, Guillem R. wrote: Dear all, I'd like to plot the marginal effect of a variable in a multiplicative interaction regression, that is, the effect of a variable conditional on the values of another variable. As an illustration, given model lm1 lm1 - lm(y ~ x*z) ? predict Perhaps: predict(lm1, newdata=data.frame(x=1:10, z=5), interval=confidence) I'd like to get the effects of x on y conditional on the values of z, with the corresponding confidence intervals if possible. Does anyone know of any package or simple way to do this? Thanks -- David Winsemius, MD West Hartford, CT __ 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] Marginal effects from interaction regression model
On Jul 25, 2010, at 11:47 PM, Guillem R. wrote: As far as I know, the predict command gives the predicted values (and intervals) of y, but what I'm looking for is the conditional effects (betas) of x on y conditional on values of z. If you want the betas, then simply print the object. Don't for get the default parameterization is for treatment effects. I'm trying to produce a plot similar to the first shown in this link: http://homepages.nyu.edu/~mrg217/interaction.html#code If you want to provide an example ...and a more complete description of the desire output, I sure someone here can finish the job of showing you how setting up a call to predict will get you to that goal. -- David. Thanks again David Winsemius wrote: On Jul 25, 2010, at 10:24 PM, Guillem R. wrote: Dear all, I'd like to plot the marginal effect of a variable in a multiplicative interaction regression, that is, the effect of a variable conditional on the values of another variable. As an illustration, given model lm1 lm1 - lm(y ~ x*z) ? predict Perhaps: predict(lm1, newdata=data.frame(x=1:10, z=5), interval=confidence) I'd like to get the effects of x on y conditional on the values of z, with the corresponding confidence intervals if possible. Does anyone know of any package or simple way to do this? Thanks -- David Winsemius, MD West Hartford, CT __ 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. -- View this message in context: http://r.789695.n4.nabble.com/Marginal-effects-from-interaction-regression-model-tp2301858p2301884.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, MD West Hartford, CT __ 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] Plot of a subset of a data.frame()
On Jul 26, 2010, at 7:38 AM, Steffen Uhlig wrote: Hello, my data.frame is sort of a collection of process values, i.e. huge run-chart. It consists of a time-stamp in the first column (date as string), factors in the following columns (used for subset- filtering), and some process-data columns. Hereafter, two examples are listed, showing the problems that occour during print: At first the example, that works fine: ~~ a = c(1:10) # create a vector of integers b = rep(c(a,b),5) # create a vector of chars, used # as factor-levels d = rnorm(10) # some random numbers e = data.frame(a,b,d) # connect to a data.frame You've gotten several answers, but none have addressed an aspect of R behavior that took me longer to appreciate than it perhaps should have. The b column inside the e data.frame is now a factor column. I mention that because you later referred to it as a string which it is not. It is an integer with an associated indexed level character vector. Many of the functions that you might think would work on strings will give either errors or unexpected results when applied to factors. e.1 = subset(e, b==a) # create two subsets e.2 = subset(e, b==b) plot(d~a, e.1, pch=3, col=2) # plot first data-subset points(d~a, e.2, pch=4, col=3) # plot the 2nd one ~~ all looks fine in theses plots. However, changing the content of vector a to a set of strings the following happens: ~~ a = c(a,b,c,d,e,f,g,h,i,j) e = data.frame(a,b,d) # re-build data.frame e.1 = subset(e, b==a) # create two subsets e.2 = subset(e, b==b) plot(d~a, e.1, pch=3, col=2) points(d~a, e.2, pch=4, col=3) ~~ The plot-command produces horizontal lines instead of dots. This seems to happen when the x-axis contains strings rather than numbers. is there a way out? Best regards, /Steffen -- David Winsemius, MD Heritage Laboratories West Hartford, CT __ 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.
