Re: [R] zero line

2013-03-18 Thread PIKAL Petr
Hi

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
 project.org] On Behalf Of bgnumis
 Sent: Sunday, March 17, 2013 4:12 PM
 To: r-help@r-project.org
 Subject: [R] zero line
 
 Hi all,
 
 I´m plotting cf (with two axis) and addind a shaded color up and down
 on the 0 line x axis (tfr1 is the time). The thing is that when I plot
 this graph adds a line up on the first plot.
 
 I hope you can understand what I mean.

Not much

 plot(cf[,1],type=l,col=black, xlab=Time, ylab=Correlation, axes=F)
Error in plot(cf[, 1], type = l, col = black, xlab = Time, ylab = 
Correlation,  : 
  object 'cf' not found

No object cf here.
It would be good if you make your code reproducible by putting a result of

dput(cf) 

and make available other objects used for plotting. Or you can try to present 
some fake data illustrating your points.

 
 How should I erase this sencond line, it is suposed they have the same
 data? Why zero line doesn´t complete the shaded area?
 
 plot(cf[,1],type=l,col=black, xlab=Time, ylab=Correlation,
 axes=F)
 grid()
 
 cord.x - c(0,1:length(cf[,1]), length(cf[,1])) cord.y - c(0,cf[,1],0)
 polygon(cord.x,cord.y,col='skyblue')
 par(new=T)

From below code you have some object named T. Hence your above code probably 
do not do what you expect as your T is put in par not reserved word TRUE.

Regards
Petr

 variable-matrix(0,length(cf[,1]),1)
 plot(tfr1[(T+Tcor):(nrow(tfr1)),1],variable,type=l,ylab=,col=black
 ,xlab=)
 mtext(Correlation,side=4,line=2,col=4)
 
 title(Dinamic Correlation,font=4)
 
   [[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] how to plot u-v wind by R?

2013-03-18 Thread Jie Tang
hi R users:
   I have a dataset including u wind in x-axis and v wind in y-axis.
 How can I plot the u,v wind data in vector or  barb figure?
 which command ?
thank you .

-- 
TANG Jie

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] how to plot u-v wind by R?

2013-03-18 Thread Jim Lemon

On 03/18/2013 08:50 PM, Jie Tang wrote:

hi R users:
I have a dataset including u wind in x-axis and v wind in y-axis.
  How can I plot the u,v wind data in vector or  barb figure?
  which command ?
thank you .


Hi Jie,
Maybe the vectorField function (plotrix)?

Jim

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] plotting CIE chromaticity diagram?

2013-03-18 Thread Bryan Hanson
As it turns out, I am working on this right now, and it's a bit of a challenge 
to get the color gradient correct.  I don't think this is really the right kind 
of question for this list, so let's see if anyone has attempted it or has 
ideas, and then take the conversation private.  I can send you my not quite 
perfect attempt a little later today.

Bryan


Prof. Bryan Hanson
Dept of Chemistry  Biochemistry
DePauw University
Greencastle IN 46135 USA
academic.depauw.edu/~hanson/deadpezsociety.html
github.com/bryanhanson
academic.depauw.edu/~hanson/UMP/Index.html

On Mar 18, 2013, at 6:12 AM, ishi soichi soichi...@gmail.com wrote:

 Has anyone plotted or is it possible to plot
 
 CIE *xy* chromaticity diagram
 
 http://en.wikipedia.org/wiki/File:CIE1931xy_blank.svg
 
 
 I need this plot in color.
 
 ishida
 
   [[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] plotting CIE chromaticity diagram?

2013-03-18 Thread Ken Knoblauch
ishi soichi soichi777 at gmail.com writes:
 Has anyone plotted or is it possible to plot
 
 CIE *xy* chromaticity diagram
 
 http://en.wikipedia.org/wiki/File:CIE1931xy_blank.svg
 
 I need this plot in color.
 
 ishida
 

I think that plotting the spectral locus and the line of purples is trivial, 
but that's not actually what you are asking.  Though filling in the interior 
with a color gradient would be a challenging and interesting problem, 
there are those who would consider this a misguided attempt at 
representation (even though it is frequently done) as the CIE coordinates 
provide no information on color appearance, per se.  
They specify what lights look alike, not what they look like.

-- 
Kenneth Knoblauch
Inserm U846
Stem-cell and Brain Research Institute
Department of Integrative Neurosciences
18 avenue du Doyen Lépine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: +33 (0)6 84 10 64 10
http://www.sbri.fr/members/kenneth-knoblauch.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.


Re: [R] plotting CIE chromaticity diagram?

2013-03-18 Thread Ken Knoblauch
ishi soichi soichi777 at gmail.com writes:
 
 Has anyone plotted or is it possible to plot
 
 CIE *xy* chromaticity diagram
 
 http://en.wikipedia.org/wiki/File:CIE1931xy_blank.svg
 
 I need this plot in color.
 
 ishida
 
sermon
And following up on my previous mail (diatribe), after
having contemplated the image at the link that you
provide, there is another significant distortion to take
into account in the representation that you propose.
Suppose that you are considering the appearance of
isolated, centrally viewed lights in isolation, in neutral
adaptation by an observer whose vision corresponds to
the average observer of CIE1931.  Then, the set of primaries 
that you will be using to generate those colors (on a screen
or on paper) will have a gamut that excludes a significant
portion of the diagram, so that if you take this limited gamut
and stretch it out to fill the diagram, then the coordinates of
the colors indicated will deviate significantly from the positions
in which they are represented (not to even speak of their
appearance).  So, in general, while this leads to a pretty
diagram, there is a lot of potential for misunderstanding
in such a (mis)representation.  /sermon

-- 
Kenneth Knoblauch
Inserm U846
Stem-cell and Brain Research Institute
Department of Integrative Neurosciences
18 avenue du Doyen Lépine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: +33 (0)6 84 10 64 10
http://www.sbri.fr/members/kenneth-knoblauch.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.


Re: [R] plotting CIE chromaticity diagram?

2013-03-18 Thread Bryan Hanson
I am, unfortunately, well-aware of the limitations that Ken points out (and I 
do appreciate him making these points).  One can readily demonstrate the gamut 
limitations by printing the diagram Ishida links to on different devices.  My 
hope is to get something close and include a disclaimer.  Bryan

On Mar 18, 2013, at 7:08 AM, Ken Knoblauch ken.knobla...@inserm.fr wrote:

 ishi soichi soichi777 at gmail.com writes:
 
 Has anyone plotted or is it possible to plot
 
 CIE *xy* chromaticity diagram
 
 http://en.wikipedia.org/wiki/File:CIE1931xy_blank.svg
 
 I need this plot in color.
 
 ishida
 
 sermon
 And following up on my previous mail (diatribe), after
 having contemplated the image at the link that you
 provide, there is another significant distortion to take
 into account in the representation that you propose.
 Suppose that you are considering the appearance of
 isolated, centrally viewed lights in isolation, in neutral
 adaptation by an observer whose vision corresponds to
 the average observer of CIE1931.  Then, the set of primaries 
 that you will be using to generate those colors (on a screen
 or on paper) will have a gamut that excludes a significant
 portion of the diagram, so that if you take this limited gamut
 and stretch it out to fill the diagram, then the coordinates of
 the colors indicated will deviate significantly from the positions
 in which they are represented (not to even speak of their
 appearance).  So, in general, while this leads to a pretty
 diagram, there is a lot of potential for misunderstanding
 in such a (mis)representation.  /sermon
 
 -- 
 Kenneth Knoblauch
 Inserm U846
 Stem-cell and Brain Research Institute
 Department of Integrative Neurosciences
 18 avenue du Doyen Lépine
 69500 Bron
 France
 tel: +33 (0)4 72 91 34 77
 fax: +33 (0)4 72 91 34 61
 portable: +33 (0)6 84 10 64 10
 http://www.sbri.fr/members/kenneth-knoblauch.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.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] plotting CIE chromaticity diagram?

2013-03-18 Thread Ken Knoblauch

Unfortunately, a bad idea that looks good never dies.

That said, I still think that it is an interesting challenge
to create this graph by programming.  I can imagine that if
you specified the coordinates at a certain number of points
in the interior and along the spectrum locus, you could use
an interpolation algorithm, to reproduce what is on the
wikipedia page.  Of course, I'm kicking myself for encouraging
this be done given the opportunity it provides for misguidance
and misrepresentation.

Ken


Quoting Bryan Hanson han...@depauw.edu:

I am, unfortunately, well-aware of the limitations that Ken points   
out (and I do appreciate him making these points).  One can readily   
demonstrate the gamut limitations by printing the diagram Ishida   
links to on different devices.  My hope is to get something close   
and include a disclaimer.  Bryan


On Mar 18, 2013, at 7:08 AM, Ken Knoblauch ken.knobla...@inserm.fr wrote:


ishi soichi soichi777 at gmail.com writes:


Has anyone plotted or is it possible to plot

CIE *xy* chromaticity diagram

http://en.wikipedia.org/wiki/File:CIE1931xy_blank.svg

I need this plot in color.

ishida


sermon
And following up on my previous mail (diatribe), after
having contemplated the image at the link that you
provide, there is another significant distortion to take
into account in the representation that you propose.
Suppose that you are considering the appearance of
isolated, centrally viewed lights in isolation, in neutral
adaptation by an observer whose vision corresponds to
the average observer of CIE1931.  Then, the set of primaries
that you will be using to generate those colors (on a screen
or on paper) will have a gamut that excludes a significant
portion of the diagram, so that if you take this limited gamut
and stretch it out to fill the diagram, then the coordinates of
the colors indicated will deviate significantly from the positions
in which they are represented (not to even speak of their
appearance).  So, in general, while this leads to a pretty
diagram, there is a lot of potential for misunderstanding
in such a (mis)representation.  /sermon

--
Kenneth Knoblauch
Inserm U846
Stem-cell and Brain Research Institute
Department of Integrative Neurosciences
18 avenue du Doyen Lépine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: +33 (0)6 84 10 64 10
http://www.sbri.fr/members/kenneth-knoblauch.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.







--
Kenneth Knoblauch
Inserm U846
Stem-cell and Brain Research Institute
Department of Integrative Neurosciences
18 avenue du Doyen Lépine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: +33 (0)6 84 10 64 10
http://www.sbri.fr/members/kenneth-knoblauch.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.


Re: [R] Q-Q plot for chi-quadrat w/o knowing the contingency table

2013-03-18 Thread Michael Dewey

At 14:56 17/03/2013, Wim Kreinen wrote:

Hello,

I have chi2-values from 2x2 contingency tables, each degree of freedom =1.
But I dont know the values of the contingency table. All I have in addtion:
The allelic odds ratio.
Now I would like to do a Q-Q-plot. Is that possible?


If I understand you correctly
?qqplot
should help you here.


If yes, please let me know. If no, please let me know as well.

Thanks Wim

[[alternative HTML version deleted]]


Michael Dewey
i...@aghmed.fsnet.co.uk
http://www.aghmed.fsnet.co.uk/home.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.


[R] How many samples ACTUALLY used in regression?

2013-03-18 Thread Federico Calboli
Dear All,

is there a simple way that covers all regression models to extract the number 
of samples from a data frame/matrix actually used in a regression model?

For instance I might have a data of 100 rows and 4 colums (1 response + 3 
explanatory variables).  If 3 samples have one or more NAs in the explanatory 
variable columns these samples will be dropped in any model:

my.model = lm(y ~ x + w + z, my.data)
my.model = glm(y ~ x + w + z, my.data, family = binomial)
my.model = polr(y ~ x + w + z, my.data)
…

I don't seem to be able to find one single method that works in the exact same 
way -- irrespective of the model type -- to interrogate my.model to see how 
many samples of my.data were actually used.  Is there such function or do I 
need to hack something together?

Best wishes

Federico

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Counting confidence intervals

2013-03-18 Thread Jim Silverton
Hi,
I have a 2 x 1 matrix of confidence intervals. The first column is the
lower and the next column is the upper. I want to cont how many times a
number say 12 lies in the interval. Can anyone assist?

-- 
Thanks,
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.


Re: [R] !0 + !0 == !0 - !0

2013-03-18 Thread Ben Bolker
Bert Gunter gunter.berton at gene.com writes:

 
 But this has nothing to do with 7.31 and everything to do with operator
 precedence and automatic casting  from integers to logical and vice-versa.
 
 I also think it fair to say that all (??) languages have these sorts of
 malapropisms due to operator precedence.
 

  Maybe FAQ 7.31 was referred to not for its direct relevance but as
a measure of the old-hand-ness of the people who will get the joke.
I have to admit it took me a minute ...

  Ben Bolker

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 many samples ACTUALLY used in regression?

2013-03-18 Thread Ben Bolker
Federico Calboli f.calboli at imperial.ac.uk writes:

 is there a simple way that covers all regression models to extract 
 the number of samples from a data
 frame/matrix actually used in a regression model?
 
 For instance I might have a data of 100 rows and 4 colums 
 (1 response + 3 explanatory variables).  If 3 samples
 have one or more NAs in the explanatory variable columns 
 these samples will be dropped in any model:

my.model = lm(y ~ x + w + z, my.data)
my.model = glm(y ~ x + w + z, my.data, family = binomial)
my.model = polr(y ~ x + w + z, my.data)

 I don't seem to be able to find one single method that works 
 in the exact same way -- irrespective of the model
 type -- to interrogate my.model to see how many samples of 
 my.data were actually used.  Is there such
 function or do I need to hack something together?

  I haven't tested it (don't want to bother to put together the
test data), but does nrow(model.frame(my.model)) work ?

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Counting confidence intervals

2013-03-18 Thread Jorge I Velez
Hi Jim,

Try either of the following (untested):

sum( x[1, ]  12  x[2, ]  12)
sum(apply(x, 2, function(x)  x[1]  12  x[2]  12))

where x is your 2x1000 matrix.

HTH,
Jorge.-


On Tue, Mar 19, 2013 at 12:03 AM, Jim Silverton  wrote:

 Hi,
 I have a 2 x 1 matrix of confidence intervals. The first column is the
 lower and the next column is the upper. I want to cont how many times a
 number say 12 lies in the interval. Can anyone assist?

 --
 Thanks,
 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.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Counting confidence intervals

2013-03-18 Thread Jeff Newmiller
sum(M[1]12  12=M[2]) untested, no data
---
Jeff NewmillerThe .   .  Go Live...
DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live Go...
  Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
--- 
Sent from my phone. Please excuse my brevity.

Jim Silverton jim.silver...@gmail.com wrote:

Hi,
I have a 2 x 1 matrix of confidence intervals. The first column is
the
lower and the next column is the upper. I want to cont how many times a
number say 12 lies in the interval. Can anyone assist?

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] header intact

2013-03-18 Thread Kamal Kishore
-- 
Thanks  Regards,

Kamal Kishore
Research  Scholar
Department of Biostatistics
National Institute of Mental Health and Neuro Sciences (NIMHANS)
Banglore-560029
Mob-+91-9591349768
Email:kamalkishorest...@gmail.com

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] try/tryCatch

2013-03-18 Thread lisa
Hi All,

I have tried every fix on my try or tryCatch that I have found on the
internet, but so far have not been able to get my R code to continue with
the for loop after the lmer model results in an error.

Here is two attemps of my code, the input is a 3D array file, but really
any function would do

metatrialstry-function(mydata){

a-matrix(data=NA, nrow=dim(mydata)[3], ncol=5)
#colnames(a)=c(sens, spec, corr, sens_se, spec_se,
counter)#colnames(a)=c(sens, spec, corr, sens_se, spec_se)
k=1
for(ii in 1:dim(mydata)[3]){
tmp-mydata[,,ii]
tmp1-as.data.frame(tmp)
names(tmp1)=c(persons, d1, tp, fn, fp, fn, detect, d0,
outcome)
lm1-try(lmer(outcome~0+d1+d0+(0+d1+d0 | persons), family=binomial,
data=tmp1, nAGQ=3), silent=T)
if(class(lm1)[1]!='try-error'){
a[ii,1]=lm1@fixef[1]
a[ii,2]=lm1@fixef[2]
a[ii,3]=vcov(lm1)[1,2]/prod(sqrt(diag(vcov(lm1
a[ii,4:5]=sqrt(diag(vcov(lm1)))
}
}
#k=k+1
#a[ii,6]=k

return(a)
}

#
# try / try catch ###
#



metatrialstry-function(mydata){

a-matrix(data=NA, nrow=dim(mydata)[3], ncol=5)
#colnames(a)=c(sens, spec, corr, sens_se, spec_se, counter)
colnames(a)=c(sens, spec, corr, sens_se, spec_se)
#a[,6]=rep(0, length(a[,6]))
for(ii in 1:dim(mydata)[3]){
tmp-mydata[,,ii]
tmp1-as.data.frame(tmp)
names(tmp1)=c(persons, d1, tp, fn, fp, fn, detect, d0,
outcome)
lm1-tryCatch(lmer(outcome~0+d1+d0+(0+d1+d0 | persons),
family=binomial, data=tmp1, nAGQ=3), error=function(e) e)
a[ii,1]=lm1@fixef[1]
a[ii,2]=lm1@fixef[2]
a[ii,3]=vcov(lm1)[1,2]/prod(sqrt(diag(vcov(lm1
a[ii,4:5]=sqrt(diag(vcov(lm1)))
}
return(a)
}


Any guidance would be greatly appreciated...

thanks!
Lisa

-- 
Lisa Wang
email: wang.li...@gmail.com
cell: +49 -0176-87786557
Tübingen, Germany, 72070

[[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] biclustering error in drawHeatmap

2013-03-18 Thread Bart Cuypers
Hi Everyone,

 

I have a problem that keeps occurring using the biclust function (info
below).

 

Does anyone have an idea what is going wrong?

 

Best Wishes and thank you,

 

Bart

 

 library(biclust)

 data-read.table(file.choose(), header = T, sep = \t)

 data2-as.matrix(data[,-1])

 data3-discretize(data2)

 dim(data3)

[1] 266  30

 clust-biclust(data3,method = BCCC(), number = 10)

 clust

 

An object of class Biclust 

 

call:

biclust(x = data3, method = BCCC(), number = 10)

 

Number of Clusters found:  10 

 

First  5  Cluster sizes:

   BC 1 BC 2 BC 3 BC 4 BC 5

Number of Rows:  82   58   39   19   15

Number of Columns:   20   15   19   16   17

 

 

 drawHeatmap(x = data3, bicResult = clust, bicluster = 10)

Error in drawHeatmap(x = data3, bicResult = clust, bicluster = 10) : 

  Error: the bicluster does not exist in the result set

 


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Problem with generated parameter estimates

2013-03-18 Thread Kamal Kishore
Dear All,

I would be very grateful for your help concerning the following question:

 Below mentioned programme is available on net to generate longitudinal
data. Usually we get almost same parameter estimates as used to generate
the data. The problem here is I am not able to get it for data used here,
despite increasing sample size and number of simulations. Is it normal to
expect this type of behavior from mixed effect models? The problem is this
code running fine for data(#data) provided by the author. I am not good
with programming, can anyone comments on the code or
possible explanation for this behavior?
Any help would be greatly appreciated.


sim.lmer - function(n,p,error,B0,varb0,B1,varb1,cor)
 {
#start function
###
### INPUTS 
###
#n is number of subjects
#p is number of time points
#error is Residuals of measurement
#B0 is fixed intercept effect (average group intercept)
#B1 is fixed slope effect (average group slope)
#varb0 is the variance of individual intercepts
#varb1 is the variance of individual slopes
#cor is the correlation between ind. intercepts and ind. slopes
# ASSUMPTIONS###
# Random effects have a bivariate normal distribution with zero mean
# Random errors have a normal distribution with zero mean
# The responses are generated based on a linear mixed effect regression
model
# Y = B0 + B1*Weeks + b0 + b1 + e
#
#
#
# Start function
cov - cor*sqrt(varb0*varb1)  #correlation between random effects
d - matrix(c(varb0,cov,cov,varb1),nrow=2,ncol=2) #var-cov matrix of random
effects
ind - mvrnorm(n,c(0,0),d)   #generate bivariate normal random
effects
b0 - ind[,1] # individual intercepts' deviation from fixed intercept
b1 - ind[,2] # individual slopes' deviation from fixed slope
d2 - (error^2)*diag(p)   # var-cov matrix of error terms at time points
err - mvrnorm(n,rep(0,p),d2) # generate multivariate normal error terms
with zero mean
ind.slo - B1+b1
ind.int - B0+b0
rand.eff - cbind(ind.slo,ind.int)
##
data - matrix(nrow=n,ncol=p)
for(i in 1:n) {
for(k in 1:p) {
data[i,k]= B0 + B1*(k-1) + (b0[i] + b1[i]*(k-1)) }}
data2 - data + err
##
data2 - as.data.frame(data2)
names - c()
for(i in 1:p) names[i]=paste(Score,i,sep=)
colnames(data2) - names ### Add column names to data
mynames - colnames(data2)
data2$ID - 1:n
d - reshape(data2,varying=mynames,idvar=ID,
v.names=Score,timevar=Weeks,times=1:p,direction=long)
list(data=d,rand.eff=rand.eff)
} ### End of Function

require(MASS)
require(plyr)
require(lme4)
data -
sim.lmer(n=1000,p=6,error=10,B0=38.93,varb0=30.92,B1=-2.31,varb1=0.56,cor=-0.7)
#data -
sim.lmer(n=1000,p=10,error=5,B0=40,varb0=1050,B1=.2,varb1=.4,cor=.7)
d - data$data
LMER.1 - lmer(Score ~ Weeks + (1 + Weeks | ID),  data = d, REML=F)
summary(LMER.1)
-- 
Thanks  Regards,

Kamal Kishore
PhD Student

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] date time manipulation- R 2.15.1 windows 7

2013-03-18 Thread yash kajaria
Hi Arun,
Thank u so much for your help. It helped me solve the problem
THANK YOU

-Yashvardhan Kajaria
On Thu, Mar 14, 2013 at 11:28 PM, arun smartpink...@yahoo.com wrote:

 HI,

 1.

 date1-c(5 jan 2013, 1 jan 2013)
  date1-as.Date(date1,format=%d %b %Y)
 date1[1]-date1[2]
 #Time difference of 4 days


 2.


 If you only have the week number of year without any other information, it
 would be difficult to predict which day that would be.
 You could get the week number from the date:
 library(lubridate)
  date2-2013-01-26
 wday(ymd(date2),label=TRUE)
 # 1 parsed with %Y-%m-%d
 #[1] Sat


 week(ymd(date2))
 # 1 parsed with %Y-%m-%d
 #[1] 4
 A.K.

 - Original Message -
 From: yash kajaria yash.kaja...@gmail.com
 To: r-help@r-project.org
 Cc:
 Sent: Thursday, March 14, 2013 9:45 AM
 Subject: [R] date  time manipulation- R 2.15.1 windows 7

 Hi,
 I wanted to learn how to solve a date and time manipulation where i can
 do the following two
 1. difference of two dates eg (differnce between 5th jan 2013 and 1st
 jan 2013)

 2.Suppose i have week number of the year, i want to know if i can find
 out the day it refers to eg( say week 2 of 2013 would be  6th jan 2013 and
 the day is sunday)
i need my result to tell me that its the 6th of jan 2013 as well as
 the day (sunday)

 Can u please help me out?

 Thanks,
 Yashvardhan Kajaria

 [[alternative HTML version deleted]]

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.



[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] save scores from sem

2013-03-18 Thread marKo
I'm not aware of any routine that those the job, although I think that 
it could be relatively easily done by multiplication the manifest 
variable vector with the estimates for the specific effect.

To make an example:

v1; v2; v3; v4 are manifest variables that loads on one y latent 
variablein a data frame called A

the code for the model should be like:
model -specifymodel(
y - v1, lam1, NA
y - v2, lam2, NA
y - v3, lam3, NA
y - v4, lam4, NA

After fitting the model with sem

model.sem - sem(model, data=A)

you should be able to compute the y variable like:
attach(data)
data$y-v1*lam1+v2*lam2+v3*lam3+v4*lam4 #change the loading name with 
the actual loading (number) or extract them from the objectiveML object 
(they are located in model.sem[[15]])


Note that those loadings are unstandardized and that the resulting 
variable will not be standardized.


Hope it helps

Regards,

Marko



--
Marko Tončić
Assistant Researcher
University of Rijeka
Faculty of Humanities and Social Sciences
Department of Psychology
Sveučilišna Avenija 4, 51000 Rijeka, Croatia

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Assign Observations to Tree Nodes

2013-03-18 Thread Jeff Coughlin
Hi, is there a way to assign observations from a dataset to nodes of a tree?  

For example, if I have a data frame like the following:

ID Var1 Var2 Var3 … VarX 
1           
.           
.           
.           
n           

And then run a classification tree with this form:
tree - ctree(Var1 ~ Var2 + Var3 + ... + VarX, data=set)

I would like to produce a vector for each node containing the observations 
(IDs) that are contained within that node.  Or really, any output that would 
tell me what observations are contained within nodes.  My goal is to run 
logits/other analyses on certain nodes and be able to easily pinpoint subsets 
of the original data frame who meet the criteria of certain nodes.

Thanks in advance for any help/suggestions.

Jeff
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Counting confidence intervals

2013-03-18 Thread arun
Hi,
Try this:
set.seed(25)
mat1- 
matrix(cbind(sample(1:15,20,replace=TRUE),sample(16:30,20,replace=TRUE)),ncol=2)
 nrow(mat1[sapply(seq_len(nrow(mat1)),function(i) 
any(seq(mat1[i,1],mat1[i,2])==12)),])
#[1] 17


set.seed(25)
mat2- 
matrix(cbind(sample(1:15,1e5,replace=TRUE),sample(16:30,1e5,replace=TRUE)),ncol=2)

system.time(res-nrow(mat2[sapply(seq_len(nrow(mat2)),function(i) 
any(seq(mat2[i,1],mat2[i,2])==12)),]))
 #  user  system elapsed 
 # 1.552   0.000   1.549 
res
#[1] 80070
 head(mat2[sapply(seq_len(nrow(mat2)),function(i) 
any(seq(mat2[i,1],mat2[i,2])==12)),])
# [,1] [,2]
#[1,]    7   29
#[2,]   11   30
#[3,]    3   30
#[4,]    2   26
#[5,]   10   22
#[6,]    6   22
A.K.







From: Jim Silverton jim.silver...@gmail.com
To: r-help@r-project.org 
Sent: Monday, March 18, 2013 9:03 AM
Subject: Re: [R] Counting confidence intervals

Hi,
I have a 2 x 1 matrix of confidence intervals. The first column is the
lower and the next column is the upper. I want to cont how many times a
number say 12 lies in the interval. Can anyone assist?

-- 
Thanks,
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.


[R] Problem with write.table

2013-03-18 Thread Sahana Srinivasan
Hi everyone,
I'm trying to create unique filenames and then write a data frame into
these files. This is the code I am using.


str1-fname
str2-.ext.txt;
FILENAME-paste0(str1,str2);

write.table(df,
file=FILENAME,col.names=FALSE,row.names=FALSE,quote=FALSE,sep=\t);



Ideally, a file called fname.ext.txt should be created, containing the data
frame df. However, the file is not getting created at all. Am I going wrong
somewhere?

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Counting confidence intervals

2013-03-18 Thread Jim Silverton
Thanks.

Jeff

On Mon, Mar 18, 2013 at 9:30 AM, Jeff Newmiller jdnew...@dcn.davis.ca.uswrote:

 sum(M[1]12  12=M[2]) untested, no data
 ---
 Jeff NewmillerThe .   .  Go Live...
 DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live
 Go...
   Live:   OO#.. Dead: OO#..  Playing
 Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
 /Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
 ---
 Sent from my phone. Please excuse my brevity.

 Jim Silverton jim.silver...@gmail.com wrote:

 Hi,
 I have a 2 x 1 matrix of confidence intervals. The first column is
 the
 lower and the next column is the upper. I want to cont how many times a
 number say 12 lies in the interval. Can anyone assist?




-- 
Thanks,
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.


Re: [R] Counting confidence intervals

2013-03-18 Thread Jorge I Velez
Thats cumbersome, Arun.

 sum(mat1[,1]  12  mat1[,2]  12)
[1] 17

will do the job and even faster:

 system.time(replicate(1, sum(mat1[,1]  12  mat1[,2]  12)))
#   user  system elapsed
#  0.067   0.001   0.078

HTH,
Jorge.-


On Tue, Mar 19, 2013 at 1:06 AM, arun  wrote:

 Hi,
 Try this:
 set.seed(25)
 mat1-
 matrix(cbind(sample(1:15,20,replace=TRUE),sample(16:30,20,replace=TRUE)),ncol=2)
  nrow(mat1[sapply(seq_len(nrow(mat1)),function(i)
 any(seq(mat1[i,1],mat1[i,2])==12)),])
 #[1] 17


 set.seed(25)
 mat2-
 matrix(cbind(sample(1:15,1e5,replace=TRUE),sample(16:30,1e5,replace=TRUE)),ncol=2)

 system.time(res-nrow(mat2[sapply(seq_len(nrow(mat2)),function(i)
 any(seq(mat2[i,1],mat2[i,2])==12)),]))
  #  user  system elapsed
  # 1.552   0.000   1.549
 res
 #[1] 80070
  head(mat2[sapply(seq_len(nrow(mat2)),function(i)
 any(seq(mat2[i,1],mat2[i,2])==12)),])
 # [,1] [,2]
 #[1,]7   29
 #[2,]   11   30
 #[3,]3   30
 #[4,]2   26
 #[5,]   10   22
 #[6,]6   22
 A.K.






 
 From: Jim Silverton 
 To: r-help@r-project.org
 Sent: Monday, March 18, 2013 9:03 AM
 Subject: Re: [R] Counting confidence intervals

 Hi,
 I have a 2 x 1 matrix of confidence intervals. The first column is the
 lower and the next column is the upper. I want to cont how many times a
 number say 12 lies in the interval. Can anyone assist?

 --
 Thanks,
 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.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Counting confidence intervals

2013-03-18 Thread arun
Hi,
Jorge's method will be faster.
#system.time(res1-sum(apply(mat2,1,function(x) x[1]12  x[2]12))) #instead 
of 2, it should be 1
 #  user  system elapsed 
 # 0.440   0.000   0.445 

 system.time(res1-sum(apply(mat2,1,function(x) x[1]=12  x[2]12))) #
 #  user  system elapsed 
 # 0.500   0.000   0.502 
 res1
#[1] 80070


A.K.







From: Jim Silverton jim.silver...@gmail.com
To: arun smartpink...@yahoo.com 
Sent: Monday, March 18, 2013 10:08 AM
Subject: Re: [R] Counting confidence intervals


thanks arun!!


On Mon, Mar 18, 2013 at 10:06 AM, arun smartpink...@yahoo.com wrote:

Hi,
Try this:
set.seed(25)
mat1- 
matrix(cbind(sample(1:15,20,replace=TRUE),sample(16:30,20,replace=TRUE)),ncol=2)
 nrow(mat1[sapply(seq_len(nrow(mat1)),function(i) 
any(seq(mat1[i,1],mat1[i,2])==12)),])
#[1] 17


set.seed(25)
mat2- 
matrix(cbind(sample(1:15,1e5,replace=TRUE),sample(16:30,1e5,replace=TRUE)),ncol=2)

system.time(res-nrow(mat2[sapply(seq_len(nrow(mat2)),function(i) 
any(seq(mat2[i,1],mat2[i,2])==12)),]))
 #  user  system elapsed
 # 1.552   0.000   1.549
res
#[1] 80070
 head(mat2[sapply(seq_len(nrow(mat2)),function(i) 
any(seq(mat2[i,1],mat2[i,2])==12)),])
# [,1] [,2]
#[1,]    7   29
#[2,]   11   30
#[3,]    3   30
#[4,]    2   26
#[5,]   10   22
#[6,]    6   22
A.K.







From: Jim Silverton jim.silver...@gmail.com
To: r-help@r-project.org
Sent: Monday, March 18, 2013 9:03 AM
Subject: Re: [R] Counting confidence intervals


Hi,
I have a 2 x 1 matrix of confidence intervals. The first column is the
lower and the next column is the upper. I want to cont how many times a
number say 12 lies in the interval. Can anyone assist?

--
Thanks,
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.



-- 
Thanks,
Jim.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Counting confidence intervals

2013-03-18 Thread Jorge I Velez
If you don't use apply() it would be even faster:

 system.time(sum(mat2[,1]  12  mat2[,2]  12))
   user  system elapsed
  0.004   0.000   0.003

Regards,
Jorge.-


On Tue, Mar 19, 2013 at 1:21 AM, arun  wrote:

 Hi,
 Jorge's method will be faster.
 #system.time(res1-sum(apply(mat2,1,function(x) x[1]12  x[2]12)))
 #instead of 2, it should be 1
  #  user  system elapsed
  # 0.440   0.000   0.445

  system.time(res1-sum(apply(mat2,1,function(x) x[1]=12  x[2]12))) #
  #  user  system elapsed
  # 0.500   0.000   0.502
  res1
 #[1] 80070


 A.K.






 
 From: Jim Silverton jim.silver...@gmail.com
 To: arun smartpink...@yahoo.com
 Sent: Monday, March 18, 2013 10:08 AM
 Subject: Re: [R] Counting confidence intervals


 thanks arun!!


 On Mon, Mar 18, 2013 at 10:06 AM, arun smartpink...@yahoo.com wrote:

 Hi,
 Try this:
 set.seed(25)
 mat1-
 matrix(cbind(sample(1:15,20,replace=TRUE),sample(16:30,20,replace=TRUE)),ncol=2)
  nrow(mat1[sapply(seq_len(nrow(mat1)),function(i)
 any(seq(mat1[i,1],mat1[i,2])==12)),])
 #[1] 17
 
 
 set.seed(25)
 mat2-
 matrix(cbind(sample(1:15,1e5,replace=TRUE),sample(16:30,1e5,replace=TRUE)),ncol=2)
 
 system.time(res-nrow(mat2[sapply(seq_len(nrow(mat2)),function(i)
 any(seq(mat2[i,1],mat2[i,2])==12)),]))
  #  user  system elapsed
  # 1.552   0.000   1.549
 res
 #[1] 80070
  head(mat2[sapply(seq_len(nrow(mat2)),function(i)
 any(seq(mat2[i,1],mat2[i,2])==12)),])
 # [,1] [,2]
 #[1,]7   29
 #[2,]   11   30
 #[3,]3   30
 #[4,]2   26
 #[5,]   10   22
 #[6,]6   22
 A.K.
 
 
 
 
 
 
 
 From: Jim Silverton jim.silver...@gmail.com
 To: r-help@r-project.org
 Sent: Monday, March 18, 2013 9:03 AM
 Subject: Re: [R] Counting confidence intervals
 
 
 Hi,
 I have a 2 x 1 matrix of confidence intervals. The first column is the
 lower and the next column is the upper. I want to cont how many times a
 number say 12 lies in the interval. Can anyone assist?
 
 --
 Thanks,
 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.
 


 --
 Thanks,
 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.


Re: [R] How many samples ACTUALLY used in regression?

2013-03-18 Thread Marc Schwartz

On Mar 18, 2013, at 7:36 AM, Federico Calboli f.calb...@imperial.ac.uk wrote:

 Dear All,
 
 is there a simple way that covers all regression models to extract the number 
 of samples from a data frame/matrix actually used in a regression model?
 
 For instance I might have a data of 100 rows and 4 colums (1 response + 3 
 explanatory variables).  If 3 samples have one or more NAs in the explanatory 
 variable columns these samples will be dropped in any model:
 
 my.model = lm(y ~ x + w + z, my.data)
 my.model = glm(y ~ x + w + z, my.data, family = binomial)
 my.model = polr(y ~ x + w + z, my.data)
 …
 
 I don't seem to be able to find one single method that works in the exact 
 same way -- irrespective of the model type -- to interrogate my.model to see 
 how many samples of my.data were actually used.  Is there such function or do 
 I need to hack something together?
 
 Best wishes
 
 Federico


I don't know that this would be universal to all possible R model 
implementations, but should work for those that at least abide by certain 
standards[1] relative to the internal use of ?model.frame.

In the case where model functions use 'model = TRUE' as the default in their 
call (eg. lm(),  glm() and MASS::polr()), the returned model object will have a 
component called 'model', such that:

  nrow(my.model$model)

returns the number of rows in the internally created data frame.

Note that 'model = TRUE' is not the default for many functions, for example 
Terry's coxph() in survival or Frank's lrm() in rms. 

Note also that the value of 'na.action' in the modeling function call may have 
a potential effect on whether the number of rows in the retained 'model' data 
frame is really the correct value.

You can also use model.frame(), independently matching arguments passed to the 
model function, to replicate what takes place internally in many modeling 
functions. The result of model.frame() will be a data frame, again, subject to 
similar limitations as above.

Regards,

Marc Schwartz

[1]: http://developer.r-project.org/model-fitting-functions.txt

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] save scores from sem

2013-03-18 Thread John Fox
Dear Marko,

I can't quite tell whether this is an original question (as suggested by no 
Re: in the subject field and no text from an earlier message) or an answer to 
a question already posed (as suggested by the phrasing of the message); if the 
latter, I apologize for missing the original question.

In any event, see the fscore() function in the sem package (?fscore) for 
computation of factor scores.

I hope this helps,
 John


John Fox
Sen. William McMaster Prof. of Social Statistics
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

On Mon, 18 Mar 2013 13:21:48 +0100
 marKo mton...@ffri.hr wrote:
 I'm not aware of any routine that those the job, although I think that it 
 could be relatively easily done by multiplication the manifest variable 
 vector with the estimates for the specific effect.
 To make an example:
 
 v1; v2; v3; v4 are manifest variables that loads on one y latent variablein a 
 data frame called A
 the code for the model should be like:
 model -specifymodel(
 y - v1, lam1, NA
 y - v2, lam2, NA
 y - v3, lam3, NA
 y - v4, lam4, NA
 
 After fitting the model with sem
 
 model.sem - sem(model, data=A)
 
 you should be able to compute the y variable like:
 attach(data)
 data$y-v1*lam1+v2*lam2+v3*lam3+v4*lam4 #change the loading name with the 
 actual loading (number) or extract them from the objectiveML object (they are 
 located in model.sem[[15]])
 
 Note that those loadings are unstandardized and that the resulting variable 
 will not be standardized.
 
 Hope it helps
 
 Regards,
 
 Marko
 
 
 
 --
 Marko Ton?i?
 Assistant Researcher
 University of Rijeka
 Faculty of Humanities and Social Sciences
 Department of Psychology
 Sveu?ilišna Avenija 4, 51000 Rijeka, Croatia
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Problem with write.table

2013-03-18 Thread PIKAL Petr
Hi

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
 project.org] On Behalf Of Sahana Srinivasan
 Sent: Monday, March 18, 2013 3:04 PM
 To: r-help@r-project.org
 Subject: [R] Problem with write.table
 
 Hi everyone,
 I'm trying to create unique filenames and then write a data frame into
 these files. This is the code I am using.
 
 
 str1-fname
 str2-.ext.txt;
 FILENAME-paste0(str1,str2);
 
 write.table(df,
 file=FILENAME,col.names=FALSE,row.names=FALSE,quote=FALSE,sep=\t);
 
 
 
 Ideally, a file called fname.ext.txt should be created, containing the
 data frame df. However, the file is not getting created at all. Am I
 going wrong somewhere?

Maybe OS and write permissions? I tried your code with some data.frame and I 
got a file fname.ext.txt without problems.

Check your objecd df and FILENAME.

Regards
Petr

 
   [[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] How many samples ACTUALLY used in regression?

2013-03-18 Thread Cade, Brian
Perhaps a crude but reliable way is to check the number of residuals, e.g.,
length(my.model$resid).

Brian

Brian S. Cade, PhD

U. S. Geological Survey
Fort Collins Science Center
2150 Centre Ave., Bldg. C
Fort Collins, CO  80526-8818

email:  ca...@usgs.gov brian_c...@usgs.gov
tel:  970 226-9326



On Mon, Mar 18, 2013 at 8:39 AM, Marc Schwartz marc_schwa...@me.com wrote:


 On Mar 18, 2013, at 7:36 AM, Federico Calboli f.calb...@imperial.ac.uk
 wrote:

  Dear All,
 
  is there a simple way that covers all regression models to extract the
 number of samples from a data frame/matrix actually used in a regression
 model?
 
  For instance I might have a data of 100 rows and 4 colums (1 response +
 3 explanatory variables).  If 3 samples have one or more NAs in the
 explanatory variable columns these samples will be dropped in any model:
 
  my.model = lm(y ~ x + w + z, my.data)
  my.model = glm(y ~ x + w + z, my.data, family = binomial)
  my.model = polr(y ~ x + w + z, my.data)
  …
 
  I don't seem to be able to find one single method that works in the
 exact same way -- irrespective of the model type -- to interrogate my.model
 to see how many samples of my.data were actually used.  Is there such
 function or do I need to hack something together?
 
  Best wishes
 
  Federico


 I don't know that this would be universal to all possible R model
 implementations, but should work for those that at least abide by certain
 standards[1] relative to the internal use of ?model.frame.

 In the case where model functions use 'model = TRUE' as the default in
 their call (eg. lm(),  glm() and MASS::polr()), the returned model object
 will have a component called 'model', such that:

   nrow(my.model$model)

 returns the number of rows in the internally created data frame.

 Note that 'model = TRUE' is not the default for many functions, for
 example Terry's coxph() in survival or Frank's lrm() in rms.

 Note also that the value of 'na.action' in the modeling function call may
 have a potential effect on whether the number of rows in the retained
 'model' data frame is really the correct value.

 You can also use model.frame(), independently matching arguments passed to
 the model function, to replicate what takes place internally in many
 modeling functions. The result of model.frame() will be a data frame,
 again, subject to similar limitations as above.

 Regards,

 Marc Schwartz

 [1]: http://developer.r-project.org/model-fitting-functions.txt

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Confirmatory factor analysis using the sem package. TLI CFI and RMSEA absent from model summary.

2013-03-18 Thread Kevin Cheung
Hi R-help,

I am using the sem package to run confirmatory factor analysis (cfa) on some 
questionnaire data collected from 307 participants. I have been running 
R-2.15.3 in Windows in conjunction with R studio. The model I am using was 
developed from exploratory factor analysis of a separate dataset (n=439); it 
includes 18 items that load onto 3 factors. I have used the sem package 
documentation and this video (http://vimeo.com/38941937) to run the cfa and 
obtain a chi-square statistic for the model. However, when I use the summary() 
function, the model does not provide indices of fit; I need these as part of my 
analysis output. In particular, I am looking for the Tucker Lewis Index (TLI), 
Comparative Fit Index (CFI),  the Root Mean Square of Approximation (RMSEA). I 
have checked the documentation and cannot seem to find any reason for this; 
none of the arguments listed with the sem command indicate that I have to 
specify these as part of the output. In addition, the analysis example f!
 rom the video includes these indices as part of the output, but my analysis 
does not provide them. I have included my code with comments below:



library(sem)

validation.data -
structure(list(V1 = c(5L, 4L, 2L, 4L, 5L, 6L, 6L, 4L, 5L, 3L,
6L, 5L, 4L, 5L, 3L, 5L, 5L, 5L, 5L, 5L, 5L, 4L, 5L, 5L, 5L, 4L,
5L, 4L, 4L, 5L, 4L, 4L, 5L, 5L, 5L, 5L, 2L, 6L, 5L, 6L, 4L, 5L,
6L, 5L, 5L, 4L, 5L, 5L, 3L, 5L, 5L, 5L, 5L, 4L, 2L, 5L, 5L, 5L,
4L, 6L, 4L, 6L, 5L, 5L, 5L, 4L, 5L, 5L, 4L, 5L, 6L, 4L, 5L, 4L,
5L, 5L, 5L, 3L, 5L, 5L, 3L, 5L, 4L, 5L, 2L, 6L, 4L, 4L, 4L, 5L,
5L, 4L, 6L, 2L, 4L, 5L, 4L, 5L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 4L, 2L, 4L, 4L, 5L, 5L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 4L,
5L, 5L, 5L, 5L, 5L, 6L, 5L, 5L, 4L, 4L, 5L, 4L, 4L, 5L, 4L, 4L,
4L, 5L, 5L, 6L, 5L, 5L, 5L, 4L, 4L, 3L, 4L, 5L, 5L, 5L, 2L, 5L,
5L, 5L, 4L, 5L, 5L, 5L, 5L, 6L, 4L, 5L, 5L, 6L, 5L, 5L, 5L, 5L,
4L, 5L, 4L, 5L, 5L, 4L, 4L, 4L, 4L, 5L, 2L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 5L, 4L, 4L, 5L, 4L, 5L, 5L, 4L, 5L, 5L, 4L, 4L, 5L,
5L, 5L, 5L, 4L, 4L, 2L, 6L, 5L, 5L, 5L, 6L, 4L, 3L, 5L, 5L, 5L,
5L, 4L, 4L, 6L, 6L, 5L, 5L, 5L, 3L, 6L, 5L, 5L, 5L, 4L, 5L, 5L,
4L, 5L, 4L, 6L, 5L, 4L, 4L, 5L, 4L, 3L, 4L, 5L, 4L, 5L, 6L, 2L,
4L, 4L, 5L, 4L, 4L, 4L, 5L, 6L, 5L, 4L, 4L, 4L, 6L, 6L, 4L, 5L,
5L, 5L, 2L, 4L, 4L, 5L, 4L, 5L, 5L, 4L, 5L, 3L, 5L, 6L, 5L, 4L,
4L, 3L, 3L, 4L, 5L, 5L, 1L, 4L, 5L, 3L, 5L, 1L, 6L, 5L, 4L, 4L,
5L, 5L, 4L, 5L, 5L, 6L, 5L, 5L, 5L), V2 = c(5L, 5L, 6L, 6L, 6L,
6L, 6L, 6L, 5L, 5L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 6L, 5L, 6L, 6L,
4L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 5L,
5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 5L, 5L, 6L, 5L, 6L, 5L,
5L, 5L, 6L, 6L, 6L, 6L, 6L, 4L, 6L, 5L, 6L, 5L, 5L, 6L, 6L, 6L,
6L, 5L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 5L,
6L, 6L, 5L, 6L, 6L, 5L, 4L, 6L, 5L, 6L, 5L, 5L, 6L, 5L, 6L, 6L,
5L, 6L, 6L, 6L, 5L, 5L, 6L, 6L, 5L, 4L, 6L, 4L, 6L, 6L, 6L, 6L,
6L, 6L, 5L, 5L, 5L, 6L, 6L, 6L, 5L, 5L, 6L, 6L, 5L, 6L, 5L, 6L,
6L, 5L, 6L, 6L, 5L, 6L, 6L, 6L, 6L, 5L, 6L, 5L, 5L, 5L, 6L, 5L,
6L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 4L, 2L, 5L, 6L, 4L,
5L, 5L, 6L, 6L, 5L, 6L, 4L, 6L, 5L, 5L, 6L, 5L, 6L, 6L, 5L, 6L,
5L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 6L, 6L, 5L, 6L, 5L, 6L, 6L,
5L, 6L, 5L, 6L, 5L, 6L, 6L, 5L, 5L, 6L, 6L, 6L, 5L, 5L, 6L, 6L,
6L, 6L, 5L, 5L, 6L, 6L, 6L, 6L, 4L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 5L, 6L, 6L, 3L, 5L, 6L, 6L, 6L, 5L, 6L, 5L, 5L, 6L, 6L,
6L, 5L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 5L, 6L, 6L, 6L, 6L, 5L, 5L,
6L, 6L, 6L, 6L, 6L, 5L, 5L, 6L, 4L, 5L, 6L, 6L, 5L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 5L, 6L, 6L, 5L, 5L, 6L, 5L, 5L, 6L, 5L, 5L, 6L,
5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 5L, 6L, 6L), V3 = c(5L,
5L, 3L, 6L, 5L, 2L, 4L, 4L, 4L, 4L, 6L, 5L, 5L, 4L, 4L, 4L, 4L,
5L, 4L, 4L, 1L, 3L, 4L, 5L, 5L, 4L, 4L, 5L, 5L, 5L, 4L, 5L, 4L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 4L, 1L, 5L, 5L, 4L, 5L, 4L, 5L,
5L, 4L, 5L, 4L, 3L, 2L, 5L, 6L, 6L, 5L, 5L, 5L, 6L, 5L, 6L, 5L,
4L, 5L, 5L, 5L, 4L, 5L, 5L, 5L, 5L, 3L, 5L, 4L, 5L, 5L, 5L, 3L,
5L, 5L, 5L, 3L, 5L, 4L, 5L, 4L, 5L, 4L, 3L, 5L, 3L, 5L, 3L, 4L,
4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 2L,
5L, 4L, 5L, 4L, 6L, 4L, 5L, 5L, 4L, 5L, 4L, 5L, 5L, 4L, 6L, 5L,
5L, 5L, 5L, 4L, 5L, 4L, 5L, 4L, 4L, 5L, 5L, 4L, 6L, 4L, 6L, 5L,
4L, 4L, 5L, 5L, 5L, 4L, 4L, 1L, 5L, 4L, 4L, 5L, 5L, 4L, 6L, 3L,
4L, 4L, 5L, 4L, 4L, 5L, 3L, 5L, 6L, 5L, 4L, 4L, 5L, 5L, 5L, 5L,
3L, 5L, 4L, 6L, 5L, 4L, 6L, 3L, 4L, 2L, 4L, 4L, 5L, 4L, 4L, 5L,
3L, 4L, 5L, 5L, 4L, 3L, 3L, 5L, 5L, 5L, 5L, 4L, 3L, 4L, 2L, 5L,
5L, 6L, 4L, 5L, 4L, 4L, 5L, 4L, 6L, 6L, 4L, 4L, 4L, 5L, 5L, 6L,
5L, 4L, 6L, 5L, 5L, 5L, 4L, 6L, 6L, 3L, 2L, 3L, 6L, 4L, 5L, 3L,
6L, 3L, 4L, 4L, 5L, 4L, 6L, 4L, 5L, 4L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 4L, 4L, 5L, 5L, 6L, 5L, 4L, 4L, 5L, 5L, 5L, 5L, 4L, 3L, 5L,
5L, 4L, 5L, 5L, 5L, 5L, 6L, 5L, 4L, 4L, 3L, 4L, 5L, 5L, 4L, 1L,
6L, 4L, 4L, 4L, 2L, 6L, 5L, 4L, 5L, 5L, 5L, 5L, 5L, 4L, 6L, 5L,
5L, 6L), V4 = c(5L, 3L, 4L, 6L, 5L, 4L, 6L, 4L, 

Re: [R] How many samples ACTUALLY used in regression?

2013-03-18 Thread Prof Brian Ripley

On 18/03/2013 14:51, Cade, Brian wrote:

Perhaps a crude but reliable way is to check the number of residuals, e.g.,
length(my.model$resid).


Not very reliable (what about zero weights, for example?), and the 
component is usually 'residuals'.


No one has so far mentioned nobs(), which seems to me to be the closest.


Brian

Brian S. Cade, PhD

U. S. Geological Survey
Fort Collins Science Center
2150 Centre Ave., Bldg. C
Fort Collins, CO  80526-8818

email:  ca...@usgs.gov brian_c...@usgs.gov
tel:  970 226-9326



On Mon, Mar 18, 2013 at 8:39 AM, Marc Schwartz marc_schwa...@me.com wrote:



On Mar 18, 2013, at 7:36 AM, Federico Calboli f.calb...@imperial.ac.uk
wrote:


Dear All,

is there a simple way that covers all regression models to extract the

number of samples from a data frame/matrix actually used in a regression
model?


For instance I might have a data of 100 rows and 4 colums (1 response +

3 explanatory variables).  If 3 samples have one or more NAs in the
explanatory variable columns these samples will be dropped in any model:


my.model = lm(y ~ x + w + z, my.data)
my.model = glm(y ~ x + w + z, my.data, family = binomial)
my.model = polr(y ~ x + w + z, my.data)
…

I don't seem to be able to find one single method that works in the

exact same way -- irrespective of the model type -- to interrogate my.model
to see how many samples of my.data were actually used.  Is there such
function or do I need to hack something together?


Best wishes

Federico



I don't know that this would be universal to all possible R model
implementations, but should work for those that at least abide by certain
standards[1] relative to the internal use of ?model.frame.

In the case where model functions use 'model = TRUE' as the default in
their call (eg. lm(),  glm() and MASS::polr()), the returned model object
will have a component called 'model', such that:

   nrow(my.model$model)

returns the number of rows in the internally created data frame.

Note that 'model = TRUE' is not the default for many functions, for
example Terry's coxph() in survival or Frank's lrm() in rms.

Note also that the value of 'na.action' in the modeling function call may
have a potential effect on whether the number of rows in the retained
'model' data frame is really the correct value.

You can also use model.frame(), independently matching arguments passed to
the model function, to replicate what takes place internally in many
modeling functions. The result of model.frame() will be a data frame,
again, subject to similar limitations as above.

Regards,

Marc Schwartz

[1]: http://developer.r-project.org/model-fitting-functions.txt

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



[[alternative HTML version deleted]]



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.




--
Brian D. Ripley,  rip...@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] How many samples ACTUALLY used in regression?

2013-03-18 Thread Federico Calboli
On 18 Mar 2013, at 15:07, Prof Brian Ripley rip...@stats.ox.ac.uk wrote:

 On 18/03/2013 14:51, Cade, Brian wrote:
 Perhaps a crude but reliable way is to check the number of residuals, e.g.,
 length(my.model$resid).
 
 Not very reliable (what about zero weights, for example?), and the component 
 is usually 'residuals'.
 
 No one has so far mentioned nobs(), which seems to me to be the closest.

Given a my.data where 3 out of 100 rows will be discarded due to NAs

test = lm(formula = y ~ x + w, my.data, model = T)
nobs(test) 
[1] 97 # as expected

But if I substitute 1 NA in one of the row of the housing data:

house.plr = polr(formula = Sat ~ Infl + Type + Cont, data = housing, weights = 
Freq)
nobs(house.plr)
[1] 1661

because of weights (which would not be take into account in a glm() fit).

Because I only care about the raw number of observations, is there a 
(hopefully) trivial way of getting nobs(poor.fit) to behave like a 
nobs(vlm.fit)?

BW

Federico






 
 Brian
 
 Brian S. Cade, PhD
 
 U. S. Geological Survey
 Fort Collins Science Center
 2150 Centre Ave., Bldg. C
 Fort Collins, CO  80526-8818
 
 email:  ca...@usgs.gov brian_c...@usgs.gov
 tel:  970 226-9326
 
 
 
 On Mon, Mar 18, 2013 at 8:39 AM, Marc Schwartz marc_schwa...@me.com wrote:
 
 
 On Mar 18, 2013, at 7:36 AM, Federico Calboli f.calb...@imperial.ac.uk
 wrote:
 
 Dear All,
 
 is there a simple way that covers all regression models to extract the
 number of samples from a data frame/matrix actually used in a regression
 model?
 
 For instance I might have a data of 100 rows and 4 colums (1 response +
 3 explanatory variables).  If 3 samples have one or more NAs in the
 explanatory variable columns these samples will be dropped in any model:
 
 my.model = lm(y ~ x + w + z, my.data)
 my.model = glm(y ~ x + w + z, my.data, family = binomial)
 my.model = polr(y ~ x + w + z, my.data)
 …
 
 I don't seem to be able to find one single method that works in the
 exact same way -- irrespective of the model type -- to interrogate my.model
 to see how many samples of my.data were actually used.  Is there such
 function or do I need to hack something together?
 
 Best wishes
 
 Federico
 
 
 I don't know that this would be universal to all possible R model
 implementations, but should work for those that at least abide by certain
 standards[1] relative to the internal use of ?model.frame.
 
 In the case where model functions use 'model = TRUE' as the default in
 their call (eg. lm(),  glm() and MASS::polr()), the returned model object
 will have a component called 'model', such that:
 
   nrow(my.model$model)
 
 returns the number of rows in the internally created data frame.
 
 Note that 'model = TRUE' is not the default for many functions, for
 example Terry's coxph() in survival or Frank's lrm() in rms.
 
 Note also that the value of 'na.action' in the modeling function call may
 have a potential effect on whether the number of rows in the retained
 'model' data frame is really the correct value.
 
 You can also use model.frame(), independently matching arguments passed to
 the model function, to replicate what takes place internally in many
 modeling functions. The result of model.frame() will be a data frame,
 again, subject to similar limitations as above.
 
 Regards,
 
 Marc Schwartz
 
 [1]: http://developer.r-project.org/model-fitting-functions.txt
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 
  [[alternative HTML version deleted]]
 
 
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 
 
 -- 
 Brian D. Ripley,  rip...@stats.ox.ac.uk
 Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
 University of Oxford, Tel:  +44 1865 272861 (self)
 1 South Parks Road, +44 1865 272866 (PA)
 Oxford OX1 3TG, UKFax:  +44 1865 272595

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] try/tryCatch

2013-03-18 Thread jim holtman
It would help if you told us what type of error you are getting and to also
provide sample data so that we could run it to see what happens.  I use
'try' a lot to catch errors and have not had any problems with it.

On Mon, Mar 18, 2013 at 6:11 AM, lisa wang.li...@gmail.com wrote:

 Hi All,

 I have tried every fix on my try or tryCatch that I have found on the
 internet, but so far have not been able to get my R code to continue with
 the for loop after the lmer model results in an error.

 Here is two attemps of my code, the input is a 3D array file, but really
 any function would do

 metatrialstry-function(mydata){

 a-matrix(data=NA, nrow=dim(mydata)[3], ncol=5)
 #colnames(a)=c(sens, spec, corr, sens_se, spec_se,
 counter)#colnames(a)=c(sens, spec, corr, sens_se, spec_se)
 k=1
 for(ii in 1:dim(mydata)[3]){
 tmp-mydata[,,ii]
 tmp1-as.data.frame(tmp)
 names(tmp1)=c(persons, d1, tp, fn, fp, fn, detect, d0,
 outcome)
 lm1-try(lmer(outcome~0+d1+d0+(0+d1+d0 | persons), family=binomial,
 data=tmp1, nAGQ=3), silent=T)
 if(class(lm1)[1]!='try-error'){
 a[ii,1]=lm1@fixef[1]
 a[ii,2]=lm1@fixef[2]
 a[ii,3]=vcov(lm1)[1,2]/prod(sqrt(diag(vcov(lm1
 a[ii,4:5]=sqrt(diag(vcov(lm1)))
 }
 }
 #k=k+1
 #a[ii,6]=k

 return(a)
 }

 #
 # try / try catch ###
 #



 metatrialstry-function(mydata){

 a-matrix(data=NA, nrow=dim(mydata)[3], ncol=5)
 #colnames(a)=c(sens, spec, corr, sens_se, spec_se, counter)
 colnames(a)=c(sens, spec, corr, sens_se, spec_se)
 #a[,6]=rep(0, length(a[,6]))
 for(ii in 1:dim(mydata)[3]){
 tmp-mydata[,,ii]
 tmp1-as.data.frame(tmp)
 names(tmp1)=c(persons, d1, tp, fn, fp, fn, detect, d0,
 outcome)
 lm1-tryCatch(lmer(outcome~0+d1+d0+(0+d1+d0 | persons),
 family=binomial, data=tmp1, nAGQ=3), error=function(e) e)
 a[ii,1]=lm1@fixef[1]
 a[ii,2]=lm1@fixef[2]
 a[ii,3]=vcov(lm1)[1,2]/prod(sqrt(diag(vcov(lm1
 a[ii,4:5]=sqrt(diag(vcov(lm1)))
 }
 return(a)
 }


 Any guidance would be greatly appreciated...

 thanks!
 Lisa

 --
 Lisa Wang
 email: wang.li...@gmail.com
 cell: +49 -0176-87786557
 Tübingen, Germany, 72070

 [[alternative HTML version deleted]]


 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.




-- 
Jim Holtman
Data Munger Guru

What is the problem that you are trying to solve?
Tell me what you want to do, not how you want to do it.

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Confirmatory factor analysis using the sem package. TLI CFI and RMSEA absent from model summary.

2013-03-18 Thread John Fox
Dear Kevin,

See ?summary.objectiveML, and in particular the description of the fit.indices 
argument. By default, the summary() method doesn't print many fit indices, but 
many are available optionally.

I hope this helps,
 John


John Fox
Sen. William McMaster Prof. of Social Statistics
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/

On Mon, 18 Mar 2013 15:00:06 +
 Kevin Cheung k.che...@derby.ac.uk wrote:
 Hi R-help,
 
 I am using the sem package to run confirmatory factor analysis (cfa) on some 
 questionnaire data collected from 307 participants. I have been running 
 R-2.15.3 in Windows in conjunction with R studio. The model I am using was 
 developed from exploratory factor analysis of a separate dataset (n=439); it 
 includes 18 items that load onto 3 factors. I have used the sem package 
 documentation and this video (http://vimeo.com/38941937) to run the cfa and 
 obtain a chi-square statistic for the model. However, when I use the 
 summary() function, the model does not provide indices of fit; I need these 
 as part of my analysis output. In particular, I am looking for the Tucker 
 Lewis Index (TLI), Comparative Fit Index (CFI),  the Root Mean Square of 
 Approximation (RMSEA). I have checked the documentation and cannot seem to 
 find any reason for this; none of the arguments listed with the sem command 
 indicate that I have to specify these as part of the output. In addition, the 
 analysis example!
  f!
  rom the video includes these indices as part of the output, but my analysis 
 does not provide them. I have included my code with comments below:
 
 
 
 library(sem)
 
 validation.data -
 structure(list(V1 = c(5L, 4L, 2L, 4L, 5L, 6L, 6L, 4L, 5L, 3L,
 6L, 5L, 4L, 5L, 3L, 5L, 5L, 5L, 5L, 5L, 5L, 4L, 5L, 5L, 5L, 4L,
 5L, 4L, 4L, 5L, 4L, 4L, 5L, 5L, 5L, 5L, 2L, 6L, 5L, 6L, 4L, 5L,
 6L, 5L, 5L, 4L, 5L, 5L, 3L, 5L, 5L, 5L, 5L, 4L, 2L, 5L, 5L, 5L,
 4L, 6L, 4L, 6L, 5L, 5L, 5L, 4L, 5L, 5L, 4L, 5L, 6L, 4L, 5L, 4L,
 5L, 5L, 5L, 3L, 5L, 5L, 3L, 5L, 4L, 5L, 2L, 6L, 4L, 4L, 4L, 5L,
 5L, 4L, 6L, 2L, 4L, 5L, 4L, 5L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
 5L, 4L, 2L, 4L, 4L, 5L, 5L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 4L,
 5L, 5L, 5L, 5L, 5L, 6L, 5L, 5L, 4L, 4L, 5L, 4L, 4L, 5L, 4L, 4L,
 4L, 5L, 5L, 6L, 5L, 5L, 5L, 4L, 4L, 3L, 4L, 5L, 5L, 5L, 2L, 5L,
 5L, 5L, 4L, 5L, 5L, 5L, 5L, 6L, 4L, 5L, 5L, 6L, 5L, 5L, 5L, 5L,
 4L, 5L, 4L, 5L, 5L, 4L, 4L, 4L, 4L, 5L, 2L, 4L, 4L, 4L, 4L, 4L,
 4L, 4L, 4L, 5L, 4L, 4L, 5L, 4L, 5L, 5L, 4L, 5L, 5L, 4L, 4L, 5L,
 5L, 5L, 5L, 4L, 4L, 2L, 6L, 5L, 5L, 5L, 6L, 4L, 3L, 5L, 5L, 5L,
 5L, 4L, 4L, 6L, 6L, 5L, 5L, 5L, 3L, 6L, 5L, 5L, 5L, 4L, 5L, 5L,
 4L, 5L, 4L, 6L, 5L, 4L, 4L, 5L, 4L, 3L, 4L, 5L, 4L, 5L, 6L, 2L,
 4L, 4L, 5L, 4L, 4L, 4L, 5L, 6L, 5L, 4L, 4L, 4L, 6L, 6L, 4L, 5L,
 5L, 5L, 2L, 4L, 4L, 5L, 4L, 5L, 5L, 4L, 5L, 3L, 5L, 6L, 5L, 4L,
 4L, 3L, 3L, 4L, 5L, 5L, 1L, 4L, 5L, 3L, 5L, 1L, 6L, 5L, 4L, 4L,
 5L, 5L, 4L, 5L, 5L, 6L, 5L, 5L, 5L), V2 = c(5L, 5L, 6L, 6L, 6L,
 6L, 6L, 6L, 5L, 5L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 6L, 5L, 6L, 6L,
 4L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 5L,
 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 5L, 5L, 6L, 5L, 6L, 5L,
 5L, 5L, 6L, 6L, 6L, 6L, 6L, 4L, 6L, 5L, 6L, 5L, 5L, 6L, 6L, 6L,
 6L, 5L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 5L,
 6L, 6L, 5L, 6L, 6L, 5L, 4L, 6L, 5L, 6L, 5L, 5L, 6L, 5L, 6L, 6L,
 5L, 6L, 6L, 6L, 5L, 5L, 6L, 6L, 5L, 4L, 6L, 4L, 6L, 6L, 6L, 6L,
 6L, 6L, 5L, 5L, 5L, 6L, 6L, 6L, 5L, 5L, 6L, 6L, 5L, 6L, 5L, 6L,
 6L, 5L, 6L, 6L, 5L, 6L, 6L, 6L, 6L, 5L, 6L, 5L, 5L, 5L, 6L, 5L,
 6L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 4L, 2L, 5L, 6L, 4L,
 5L, 5L, 6L, 6L, 5L, 6L, 4L, 6L, 5L, 5L, 6L, 5L, 6L, 6L, 5L, 6L,
 5L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 6L, 6L, 5L, 6L, 5L, 6L, 6L,
 5L, 6L, 5L, 6L, 5L, 6L, 6L, 5L, 5L, 6L, 6L, 6L, 5L, 5L, 6L, 6L,
 6L, 6L, 5L, 5L, 6L, 6L, 6L, 6L, 4L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
 6L, 6L, 5L, 6L, 6L, 3L, 5L, 6L, 6L, 6L, 5L, 6L, 5L, 5L, 6L, 6L,
 6L, 5L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 5L, 6L, 6L, 6L, 6L, 5L, 5L,
 6L, 6L, 6L, 6L, 6L, 5L, 5L, 6L, 4L, 5L, 6L, 6L, 5L, 6L, 6L, 6L,
 6L, 6L, 6L, 6L, 5L, 6L, 6L, 5L, 5L, 6L, 5L, 5L, 6L, 5L, 5L, 6L,
 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 5L, 6L, 6L), V3 = c(5L,
 5L, 3L, 6L, 5L, 2L, 4L, 4L, 4L, 4L, 6L, 5L, 5L, 4L, 4L, 4L, 4L,
 5L, 4L, 4L, 1L, 3L, 4L, 5L, 5L, 4L, 4L, 5L, 5L, 5L, 4L, 5L, 4L,
 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 4L, 1L, 5L, 5L, 4L, 5L, 4L, 5L,
 5L, 4L, 5L, 4L, 3L, 2L, 5L, 6L, 6L, 5L, 5L, 5L, 6L, 5L, 6L, 5L,
 4L, 5L, 5L, 5L, 4L, 5L, 5L, 5L, 5L, 3L, 5L, 4L, 5L, 5L, 5L, 3L,
 5L, 5L, 5L, 3L, 5L, 4L, 5L, 4L, 5L, 4L, 3L, 5L, 3L, 5L, 3L, 4L,
 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 2L,
 5L, 4L, 5L, 4L, 6L, 4L, 5L, 5L, 4L, 5L, 4L, 5L, 5L, 4L, 6L, 5L,
 5L, 5L, 5L, 4L, 5L, 4L, 5L, 4L, 4L, 5L, 5L, 4L, 6L, 4L, 6L, 5L,
 4L, 4L, 5L, 5L, 5L, 4L, 4L, 1L, 5L, 4L, 4L, 5L, 5L, 4L, 6L, 3L,
 4L, 4L, 5L, 4L, 4L, 5L, 3L, 

Re: [R] Counting confidence intervals

2013-03-18 Thread S Ellison
  I want to cont how many 
  times a number say 12 lies in the interval. Can anyone assist?

Has anyone else ever wished there was a moderately general 'inside' or 'within' 
function in R for this problem?

For example, something that behaves more or less like 

within - function(x, interval=NULL, closed=c(TRUE, TRUE), lower=min(interval), 
upper=max(interval)) {
#interval must be a length 2 vector
#closed is taken in the order (lower, upper)
#lower and upper may be vectors and will be recycled (by  etc) if 
not of length length(x)

low.comp - if(closed[1]) = else  
high.comp - if(closed[2]) = else 

do.call(low.comp, list(lower, x))  do.call(high.comp, list(upper, x))
}


#Examples
within(1:5, c(2,4))

within(1:5, c(2,4), closed=c(FALSE, TRUE))

within(1:5, lower=5:1, upper=10:14)


S Ellison
LGC

***
This email and any attachments are confidential. Any use...{{dropped:8}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] !0 + !0 == !0 - !0

2013-03-18 Thread Charles Berry
Ben Bolker bbolker at gmail.com writes:


   Maybe FAQ 7.31 was referred to not for its direct relevance but as
 a measure of the old-hand-ness of the people who will get the joke.


   !1i|!0

Chuck

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Confirmatory factor analysis using the sem package. TLI CFI and RMSEA absent from model summary.

2013-03-18 Thread Kevin Cheung
Dear John,

Thank you for taking the time to help me with this, I have been able to obtain 
fit indices using the information that you provided.

Note to users searching archived R-help posts about this issue: The 
instructional video I looked at (http://vimeo.com/38941937) gave fit indices 
using the default summary() command without any additional arguments. This may 
have been due to a different version of R (I noticed that the instructor was 
using a mac based OS).

With regards,
Kevin

Kevin Yet Fong Cheung, Bsc., MRes., MBPsS.
Postgraduate Researcher
Centre for Psychological Research
University of Derby
Kedleston Road
Derby DE22 1GB
k.che...@derby.ac.ukmailto:k.che...@derby.ac.uk
01332592081

http://derby.academia.edu/KevinCheung


-Original Message-
From: John Fox [mailto:j...@mcmaster.ca]
Sent: 18 March 2013 15:55
To: Kevin Cheung
Cc: r-help@r-project.org
Subject: Re: [R] Confirmatory factor analysis using the sem package. TLI CFI 
and RMSEA absent from model summary.

Dear Kevin,

See ?summary.objectiveML, and in particular the description of the fit.indices 
argument. By default, the summary() method doesn't print many fit indices, but 
many are available optionally.

I hope this helps,
 John


John Fox
Sen. William McMaster Prof. of Social Statistics Department of Sociology 
McMaster University Hamilton, Ontario, Canada http://socserv.mcmaster.ca/jfox/

On Mon, 18 Mar 2013 15:00:06 +
 Kevin Cheung k.che...@derby.ac.uk wrote:
 Hi R-help,

 I am using the sem package to run confirmatory factor analysis (cfa) on some 
 questionnaire data collected from 307 participants. I have been running 
 R-2.15.3 in Windows in conjunction with R studio. The model I am using was 
 developed from exploratory factor analysis of a separate dataset (n=439); it 
 includes 18 items that load onto 3 factors. I have used the sem package 
 documentation and this video (http://vimeo.com/38941937) to run the cfa and 
 obtain a chi-square statistic for the model. However, when I use the 
 summary() function, the model does not provide indices of fit; I need these 
 as part of my analysis output. In particular, I am looking for the Tucker 
 Lewis Index (TLI), Comparative Fit Index (CFI),  the Root Mean Square of 
 Approximation (RMSEA). I have checked the documentation and cannot seem to 
 find any reason for this; none of the arguments listed with the sem command 
 indicate that I have to specify these as part of the output. In addition, the 
 analysis example!
  f!
  rom the video includes these indices as part of the output, but my analysis 
 does not provide them. I have included my code with comments below:

 

 library(sem)

 validation.data -
 structure(list(V1 = c(5L, 4L, 2L, 4L, 5L, 6L, 6L, 4L, 5L, 3L, 6L, 5L,
 4L, 5L, 3L, 5L, 5L, 5L, 5L, 5L, 5L, 4L, 5L, 5L, 5L, 4L, 5L, 4L, 4L,
 5L, 4L, 4L, 5L, 5L, 5L, 5L, 2L, 6L, 5L, 6L, 4L, 5L, 6L, 5L, 5L, 4L,
 5L, 5L, 3L, 5L, 5L, 5L, 5L, 4L, 2L, 5L, 5L, 5L, 4L, 6L, 4L, 6L, 5L,
 5L, 5L, 4L, 5L, 5L, 4L, 5L, 6L, 4L, 5L, 4L, 5L, 5L, 5L, 3L, 5L, 5L,
 3L, 5L, 4L, 5L, 2L, 6L, 4L, 4L, 4L, 5L, 5L, 4L, 6L, 2L, 4L, 5L, 4L,
 5L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 4L, 2L, 4L, 4L, 5L, 5L, 4L,
 4L, 5L, 5L, 5L, 5L, 5L, 5L, 4L, 5L, 5L, 5L, 5L, 5L, 6L, 5L, 5L, 4L,
 4L, 5L, 4L, 4L, 5L, 4L, 4L, 4L, 5L, 5L, 6L, 5L, 5L, 5L, 4L, 4L, 3L,
 4L, 5L, 5L, 5L, 2L, 5L, 5L, 5L, 4L, 5L, 5L, 5L, 5L, 6L, 4L, 5L, 5L,
 6L, 5L, 5L, 5L, 5L, 4L, 5L, 4L, 5L, 5L, 4L, 4L, 4L, 4L, 5L, 2L, 4L,
 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 4L, 4L, 5L, 4L, 5L, 5L, 4L, 5L, 5L,
 4L, 4L, 5L, 5L, 5L, 5L, 4L, 4L, 2L, 6L, 5L, 5L, 5L, 6L, 4L, 3L, 5L,
 5L, 5L, 5L, 4L, 4L, 6L, 6L, 5L, 5L, 5L, 3L, 6L, 5L, 5L, 5L, 4L, 5L,
 5L, 4L, 5L, 4L, 6L, 5L, 4L, 4L, 5L, 4L, 3L, 4L, 5L, 4L, 5L, 6L, 2L,
 4L, 4L, 5L, 4L, 4L, 4L, 5L, 6L, 5L, 4L, 4L, 4L, 6L, 6L, 4L, 5L, 5L,
 5L, 2L, 4L, 4L, 5L, 4L, 5L, 5L, 4L, 5L, 3L, 5L, 6L, 5L, 4L, 4L, 3L,
 3L, 4L, 5L, 5L, 1L, 4L, 5L, 3L, 5L, 1L, 6L, 5L, 4L, 4L, 5L, 5L, 4L,
 5L, 5L, 6L, 5L, 5L, 5L), V2 = c(5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 5L,
 5L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 6L, 5L, 6L, 6L, 4L, 5L, 6L, 6L, 6L,
 6L, 6L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 5L, 5L, 5L, 6L, 6L, 6L, 6L,
 6L, 6L, 6L, 6L, 5L, 5L, 6L, 5L, 6L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L,
 4L, 6L, 5L, 6L, 5L, 5L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 6L, 6L, 6L, 5L,
 6L, 6L, 6L, 6L, 6L, 5L, 6L, 5L, 6L, 6L, 5L, 6L, 6L, 5L, 4L, 6L, 5L,
 6L, 5L, 5L, 6L, 5L, 6L, 6L, 5L, 6L, 6L, 6L, 5L, 5L, 6L, 6L, 5L, 4L,
 6L, 4L, 6L, 6L, 6L, 6L, 6L, 6L, 5L, 5L, 5L, 6L, 6L, 6L, 5L, 5L, 6L,
 6L, 5L, 6L, 5L, 6L, 6L, 5L, 6L, 6L, 5L, 6L, 6L, 6L, 6L, 5L, 6L, 5L,
 5L, 5L, 6L, 5L, 6L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 4L, 2L,
 5L, 6L, 4L, 5L, 5L, 6L, 6L, 5L, 6L, 4L, 6L, 5L, 5L, 6L, 5L, 6L, 6L,
 5L, 6L, 5L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 6L, 6L, 5L, 6L, 5L, 6L,
 6L, 5L, 6L, 5L, 6L, 5L, 6L, 6L, 5L, 5L, 6L, 6L, 6L, 5L, 5L, 6L, 6L,
 6L, 6L, 5L, 5L, 6L, 6L, 6L, 6L, 4L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
 6L, 5L, 

Re: [R] !0 + !0 == !0 - !0

2013-03-18 Thread Sam Steingold
 * Bert Gunter thagre.ore...@trar.pbz [2013-03-17 20:30:56 -0700]:

 I also think it fair to say that all (??) languages have these sorts
 of malapropisms due to operator precedence.

Except for those languages which do _not_ have operator precedence.
Like, e.g., Lisp.

-- 
Sam Steingold (http://sds.podval.org/) on Ubuntu 12.04 (precise) X 11.0.11103000
http://www.childpsy.net/ http://openvotingconsortium.org http://jihadwatch.org
http://palestinefacts.org http://mideasttruth.com http://camera.org
DRM access management == prison freedom management.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] !0 + !0 == !0 - !0

2013-03-18 Thread Bert Gunter
Sam:

Yes.  Good point. (which is why my ?? was necessary).

-- Bert

On Mon, Mar 18, 2013 at 10:05 AM, Sam Steingold s...@gnu.org wrote:

  * Bert Gunter thagre.ore...@trar.pbz [2013-03-17 20:30:56 -0700]:
 
  I also think it fair to say that all (??) languages have these sorts
  of malapropisms due to operator precedence.

 Except for those languages which do _not_ have operator precedence.
 Like, e.g., Lisp.

 --
 Sam Steingold (http://sds.podval.org/) on Ubuntu 12.04 (precise) X
 11.0.11103000
 http://www.childpsy.net/ http://openvotingconsortium.org
 http://jihadwatch.org
 http://palestinefacts.org http://mideasttruth.com http://camera.org
 DRM access management == prison freedom management.




-- 

Bert Gunter
Genentech Nonclinical Biostatistics

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Writing a hyperlink to a csv file

2013-03-18 Thread MacQueen, Don
Besides what others have suggested, there are at some packages for writing
Excel files directly, and either of those might have a way to write
something that Excel will recognize as a hyperlink.

xlsx, XLConnect are both Java-based; I would start with these.

-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 3/16/13 6:28 AM, Brian Smith bsmith030...@gmail.com wrote:

Hi Marc,

Thanks for the reply.

The question is whether it is possible to write some text (with a
hyperlink) to a csv file, such that when you open the file in excel, it
shows the text as hyperlinked. I guess it boils down to whether there are
any 'tags' that you can put in the csv/txt file so that excel recognizes
it
as hyperlinked text.

Does that make sense?

thanks!

On Sat, Mar 16, 2013 at 4:16 AM, Marc Girondot marc_...@yahoo.fr wrote:

 Le 15/03/13 12:53, Brian Smith a écrit :

  Hi,

 I was wondering if it is possible to create a hyperlink in a csv file
 using
 R code and some package. For example, in the following code:

 links - cbind(rep('Click for Google',3),google search address goes
 here)
 ## R Mailing list blocks if I put the actual web address here
 write.table(links,'test.csv',
 sep=',',row.names=F,col.names=**F)


 the web address should be linked to 'Click for Google'.

 The browseURL() function open your internet browser with the url
indicated
 as parameter:
 browseURL(http://www.r-**project.org http://www.r-project.org)

 But I am not sure how you want call it. You should be more precise about
 the context you want to use it.
 For example:

 links - data.frame(c(' for Google', ' for Bing'),
c(http://www.google.com;, http://www.bing.com;),
 stringsAsFactors = FALSE)
 cat(Choose an option:\n, paste(1:2, links[,1],\n))
 f-scan(nmax=1, quiet=TRUE)

 browseURL(links[f,2])

 Sincerely
 Marc

 --
 __**
 Marc Girondot, Pr

 Laboratoire Ecologie, Systématique et Evolution
 Equipe de Conservation des Populations et des Communautés
 CNRS, AgroParisTech et Université Paris-Sud 11 , UMR 8079
 Bâtiment 362
 91405 Orsay Cedex, France

 Tel:  33 1 (0)1.69.15.72.30   Fax: 33 1 (0)1.69.15.73.53
 e-mail: marc.giron...@u-psud.fr
 Web: 
http://www.ese.u-psud.fr/epc/**conservation/Marc.htmlhttp://www.ese.u-ps
ud.fr/epc/conservation/Marc.html
 Skype: girondot



   [[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] melt with complications

2013-03-18 Thread Richard M. Heiberger
## Can someone suggest a simpler expression than either of these, with the
goal
## of taking a long matrix into a wide one with exactly one of the factors
converted to
## columns and all the rest retained as factors.  I want something that
generalizes beyond
## the three factors illustrated here.

## Rich

meltTest - data.frame(A=rep(c(B,C), each=12),
   D=rep(c(E,F,G), each=4, times=2),
   H=rep(c(I,J,K,L), times=6),
   M=1:24)

meltTest

result.melt - do.call(rbind, {
  tmp - cast(D ~ H | A, value=M, data=meltTest)
  lapply(names(tmp), function(x) cbind(A=x, tmp[[x]])) ## explicit use of
name A
})
result.melt

result.reshape - reshape(meltTest, direction=wide, timevar=H,
idvar=c(A,D))
names(result.reshape)[3:6] - unique(as.character(meltTest$H))  ## explicit
use of name H
result.reshape

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] try/tryCatch

2013-03-18 Thread lisa
here is the error:

 aa-metatrialstry(beta_5_50)

Error in asMethod(object) : matrix is not symmetric [1,2]



metatrials, the function that i am attempting to convert with try/tryCatch,
gives me back a matrix with as many rows are there are simulations  (z) in
the aray with dim(x,y,z). with the data i attached, x is 500(number of
patients), y is 9 (these are covariates), and z is 500.


metatrials-function(mydata){

a-matrix(data=NA, nrow=dim(mydata)[3], ncol=5)
colnames(a)=c(sens, spec, corr, sens_se, spec_se)

for(ii in 1:dim(mydata)[3]){
tmp-mydata[,,ii]
tmp1-as.data.frame(tmp)
names(tmp1)=c(persons, d1, tp, fn, fp, fn, detect, d0,
outcome)
lm1-lmer(outcome~0+d1+d0+(0+d1+d0 | persons), family=binomial,
data=tmp1, nAGQ=3)
a[ii,1]=lm1@fixef[1]
a[ii,2]=lm1@fixef[2]
a[ii,3]=vcov(lm1)[1,2]/prod(sqrt(diag(vcov(lm1
a[ii,4:5]=sqrt(diag(vcov(lm1)))

}
return(a)
}


##


what i want is for the function to go on to the next data set in the array
and simply return an NA for that line in the metatrials results. so
basically, just keep going.


thanks so much for your help!


-Lisa


On Mon, Mar 18, 2013 at 4:24 PM, jim holtman jholt...@gmail.com wrote:

 It would help if you told us what type of error you are getting and to
 also provide sample data so that we could run it to see what happens.  I
 use 'try' a lot to catch errors and have not had any problems with it.

 On Mon, Mar 18, 2013 at 6:11 AM, lisa wang.li...@gmail.com wrote:

 Hi All,

 I have tried every fix on my try or tryCatch that I have found on the
 internet, but so far have not been able to get my R code to continue with
 the for loop after the lmer model results in an error.

 Here is two attemps of my code, the input is a 3D array file, but really
 any function would do

 metatrialstry-function(mydata){

 a-matrix(data=NA, nrow=dim(mydata)[3], ncol=5)
 #colnames(a)=c(sens, spec, corr, sens_se, spec_se,
 counter)#colnames(a)=c(sens, spec, corr, sens_se, spec_se)
 k=1
 for(ii in 1:dim(mydata)[3]){
 tmp-mydata[,,ii]
 tmp1-as.data.frame(tmp)
 names(tmp1)=c(persons, d1, tp, fn, fp, fn, detect, d0,
 outcome)
 lm1-try(lmer(outcome~0+d1+d0+(0+d1+d0 | persons), family=binomial,
 data=tmp1, nAGQ=3), silent=T)
 if(class(lm1)[1]!='try-error'){
 a[ii,1]=lm1@fixef[1]
 a[ii,2]=lm1@fixef[2]
 a[ii,3]=vcov(lm1)[1,2]/prod(sqrt(diag(vcov(lm1
 a[ii,4:5]=sqrt(diag(vcov(lm1)))
 }
 }
 #k=k+1
 #a[ii,6]=k

 return(a)
 }

 #
 # try / try catch ###
 #



 metatrialstry-function(mydata){

 a-matrix(data=NA, nrow=dim(mydata)[3], ncol=5)
 #colnames(a)=c(sens, spec, corr, sens_se, spec_se, counter)
 colnames(a)=c(sens, spec, corr, sens_se, spec_se)
 #a[,6]=rep(0, length(a[,6]))
 for(ii in 1:dim(mydata)[3]){
 tmp-mydata[,,ii]
 tmp1-as.data.frame(tmp)
 names(tmp1)=c(persons, d1, tp, fn, fp, fn, detect, d0,
 outcome)
 lm1-tryCatch(lmer(outcome~0+d1+d0+(0+d1+d0 | persons),
 family=binomial, data=tmp1, nAGQ=3), error=function(e) e)
 a[ii,1]=lm1@fixef[1]
 a[ii,2]=lm1@fixef[2]
 a[ii,3]=vcov(lm1)[1,2]/prod(sqrt(diag(vcov(lm1
 a[ii,4:5]=sqrt(diag(vcov(lm1)))
 }
 return(a)
 }


 Any guidance would be greatly appreciated...

 thanks!
 Lisa

 --
 Lisa Wang
 email: wang.li...@gmail.com
 cell: +49 -0176-87786557
 Tübingen, Germany, 72070

 [[alternative HTML version deleted]]


 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.




 --
 Jim Holtman
 Data Munger Guru

 What is the problem that you are trying to solve?
 Tell me what you want to do, not how you want to do it.




-- 
Lisa Wang
email: wang.li...@gmail.com
cell: +49 -0176-87786557
Tübingen, Germany, 72070
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Fit a mixture of lognormal and normal distributions

2013-03-18 Thread To . .

Hello
I am trying to find an automated way of fitting a mixture of normal and 
log-normal distributions to data which is clearly bimodal.
Here's a simulated example:
x.1-rnorm(6000, 2.4, 0.6)x.2-rlnorm(1, 1.3,0.1)X-c(x.1, x.2)
hist(X,100,freq=FALSE, ylim=c(0,1.5))lines(density(x.1), lty=2, 
lwd=2)lines(density(x.2), lty=2, lwd=2)lines(density(X), lty=4)

Currently i am using mixtools and just run:
library(mixtools)mixmdl = normalmixEM(X, k=2, epsilon = 1e-08, maxit = 1000, 
maxrestarts=20, verb = TRUE, fast=FALSE, ECM = FALSE, arbmean = TRUE, arbvar = 
TRUE) plot(mixmdl,which=2)lines(density(X), lty=2, lwd=2)
This is obviously not the best way of doing this. The estimates it gives me 
are:mu 3.6595737(x.2 log()=1.29) 2.3113135(x.1)
These are not too far off but I was wondering if someone knows of a better R 
package/way of doing this or has any hints that would help me coding it from 
scratch (EM+writing down the formulae for ML? but... would this really be 
better than using a more-advanced-already-available-package-made-by-pros?).
The objective is to estimate threshold values at specific FDRs. (some help with 
that would also be most helpful.)
Thanks to all in advance!To'
 


  
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] data.frame with NA

2013-03-18 Thread Pete

I have this little data.frame

http://dl.dropbox.com/u/102669/nanotna.rdata

Two column contains NA, so the best thing to do is use na.locf function (with
fromLast = T)

But locf function doesn't work because NA in my data.frame are not recognized as
real NA.

Is there a way to substitute fake NA with real NA? In this case na.locf function
should work

Thank you

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Confirmatory factor analysis using the sem package. TLI CFI and RMSEA absent from model summary.

2013-03-18 Thread Hervé Guyon

Hi Kevin,

With sem package use : summary(X,fit.indices=c(RMSEA,...)) to get 
RMSEA or anothers fit indices.

See the section ML.methods in sem.pdf

Hervé

Hervé
Le 18/03/2013 16:00, Kevin Cheung a écrit :

Hi R-help,

I am using the sem package to run confirmatory factor analysis (cfa) on some 
questionnaire data collected from 307 participants. I have been running R-2.15.3 in 
Windows in conjunction with R studio. The model I am using was developed from 
exploratory factor analysis of a separate dataset (n=439); it includes 18 items 
that load onto 3 factors. I have used the sem package documentation and this video 
(http://vimeo.com/38941937) to run the cfa and obtain a chi-square statistic for 
the model. However, when I use the summary() function, the model does not provide 
indices of fit; I need these as part of my analysis output. In particular, I am 
looking for the Tucker Lewis Index (TLI), Comparative Fit Index (CFI),  the 
Root Mean Square of Approximation (RMSEA). I have checked the documentation and 
cannot seem to find any reason for this; none of the arguments listed with the sem 
command indicate that I have to specify these as part of the output. In addition, 
the analysis example f!
  rom the video includes these indices as part of the output, but my analysis 
does not provide them. I have included my code with comments below:



library(sem)

validation.data -
structure(list(V1 = c(5L, 4L, 2L, 4L, 5L, 6L, 6L, 4L, 5L, 3L,
6L, 5L, 4L, 5L, 3L, 5L, 5L, 5L, 5L, 5L, 5L, 4L, 5L, 5L, 5L, 4L,
5L, 4L, 4L, 5L, 4L, 4L, 5L, 5L, 5L, 5L, 2L, 6L, 5L, 6L, 4L, 5L,
6L, 5L, 5L, 4L, 5L, 5L, 3L, 5L, 5L, 5L, 5L, 4L, 2L, 5L, 5L, 5L,
4L, 6L, 4L, 6L, 5L, 5L, 5L, 4L, 5L, 5L, 4L, 5L, 6L, 4L, 5L, 4L,
5L, 5L, 5L, 3L, 5L, 5L, 3L, 5L, 4L, 5L, 2L, 6L, 4L, 4L, 4L, 5L,
5L, 4L, 6L, 2L, 4L, 5L, 4L, 5L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 4L, 2L, 4L, 4L, 5L, 5L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 4L,
5L, 5L, 5L, 5L, 5L, 6L, 5L, 5L, 4L, 4L, 5L, 4L, 4L, 5L, 4L, 4L,
4L, 5L, 5L, 6L, 5L, 5L, 5L, 4L, 4L, 3L, 4L, 5L, 5L, 5L, 2L, 5L,
5L, 5L, 4L, 5L, 5L, 5L, 5L, 6L, 4L, 5L, 5L, 6L, 5L, 5L, 5L, 5L,
4L, 5L, 4L, 5L, 5L, 4L, 4L, 4L, 4L, 5L, 2L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 5L, 4L, 4L, 5L, 4L, 5L, 5L, 4L, 5L, 5L, 4L, 4L, 5L,
5L, 5L, 5L, 4L, 4L, 2L, 6L, 5L, 5L, 5L, 6L, 4L, 3L, 5L, 5L, 5L,
5L, 4L, 4L, 6L, 6L, 5L, 5L, 5L, 3L, 6L, 5L, 5L, 5L, 4L, 5L, 5L,
4L, 5L, 4L, 6L, 5L, 4L, 4L, 5L, 4L, 3L, 4L, 5L, 4L, 5L, 6L, 2L,
4L, 4L, 5L, 4L, 4L, 4L, 5L, 6L, 5L, 4L, 4L, 4L, 6L, 6L, 4L, 5L,
5L, 5L, 2L, 4L, 4L, 5L, 4L, 5L, 5L, 4L, 5L, 3L, 5L, 6L, 5L, 4L,
4L, 3L, 3L, 4L, 5L, 5L, 1L, 4L, 5L, 3L, 5L, 1L, 6L, 5L, 4L, 4L,
5L, 5L, 4L, 5L, 5L, 6L, 5L, 5L, 5L), V2 = c(5L, 5L, 6L, 6L, 6L,
6L, 6L, 6L, 5L, 5L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 6L, 5L, 6L, 6L,
4L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 5L,
5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 5L, 5L, 6L, 5L, 6L, 5L,
5L, 5L, 6L, 6L, 6L, 6L, 6L, 4L, 6L, 5L, 6L, 5L, 5L, 6L, 6L, 6L,
6L, 5L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 5L,
6L, 6L, 5L, 6L, 6L, 5L, 4L, 6L, 5L, 6L, 5L, 5L, 6L, 5L, 6L, 6L,
5L, 6L, 6L, 6L, 5L, 5L, 6L, 6L, 5L, 4L, 6L, 4L, 6L, 6L, 6L, 6L,
6L, 6L, 5L, 5L, 5L, 6L, 6L, 6L, 5L, 5L, 6L, 6L, 5L, 6L, 5L, 6L,
6L, 5L, 6L, 6L, 5L, 6L, 6L, 6L, 6L, 5L, 6L, 5L, 5L, 5L, 6L, 5L,
6L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 4L, 2L, 5L, 6L, 4L,
5L, 5L, 6L, 6L, 5L, 6L, 4L, 6L, 5L, 5L, 6L, 5L, 6L, 6L, 5L, 6L,
5L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 6L, 6L, 5L, 6L, 5L, 6L, 6L,
5L, 6L, 5L, 6L, 5L, 6L, 6L, 5L, 5L, 6L, 6L, 6L, 5L, 5L, 6L, 6L,
6L, 6L, 5L, 5L, 6L, 6L, 6L, 6L, 4L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 5L, 6L, 6L, 3L, 5L, 6L, 6L, 6L, 5L, 6L, 5L, 5L, 6L, 6L,
6L, 5L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 5L, 6L, 6L, 6L, 6L, 5L, 5L,
6L, 6L, 6L, 6L, 6L, 5L, 5L, 6L, 4L, 5L, 6L, 6L, 5L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 5L, 6L, 6L, 5L, 5L, 6L, 5L, 5L, 6L, 5L, 5L, 6L,
5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 5L, 6L, 6L), V3 = c(5L,
5L, 3L, 6L, 5L, 2L, 4L, 4L, 4L, 4L, 6L, 5L, 5L, 4L, 4L, 4L, 4L,
5L, 4L, 4L, 1L, 3L, 4L, 5L, 5L, 4L, 4L, 5L, 5L, 5L, 4L, 5L, 4L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 4L, 1L, 5L, 5L, 4L, 5L, 4L, 5L,
5L, 4L, 5L, 4L, 3L, 2L, 5L, 6L, 6L, 5L, 5L, 5L, 6L, 5L, 6L, 5L,
4L, 5L, 5L, 5L, 4L, 5L, 5L, 5L, 5L, 3L, 5L, 4L, 5L, 5L, 5L, 3L,
5L, 5L, 5L, 3L, 5L, 4L, 5L, 4L, 5L, 4L, 3L, 5L, 3L, 5L, 3L, 4L,
4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 2L,
5L, 4L, 5L, 4L, 6L, 4L, 5L, 5L, 4L, 5L, 4L, 5L, 5L, 4L, 6L, 5L,
5L, 5L, 5L, 4L, 5L, 4L, 5L, 4L, 4L, 5L, 5L, 4L, 6L, 4L, 6L, 5L,
4L, 4L, 5L, 5L, 5L, 4L, 4L, 1L, 5L, 4L, 4L, 5L, 5L, 4L, 6L, 3L,
4L, 4L, 5L, 4L, 4L, 5L, 3L, 5L, 6L, 5L, 4L, 4L, 5L, 5L, 5L, 5L,
3L, 5L, 4L, 6L, 5L, 4L, 6L, 3L, 4L, 2L, 4L, 4L, 5L, 4L, 4L, 5L,
3L, 4L, 5L, 5L, 4L, 3L, 3L, 5L, 5L, 5L, 5L, 4L, 3L, 4L, 2L, 5L,
5L, 6L, 4L, 5L, 4L, 4L, 5L, 4L, 6L, 6L, 4L, 4L, 4L, 5L, 5L, 6L,
5L, 4L, 6L, 5L, 5L, 5L, 4L, 6L, 6L, 3L, 2L, 3L, 6L, 4L, 5L, 3L,
6L, 3L, 4L, 4L, 5L, 4L, 6L, 4L, 5L, 4L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 4L, 4L, 5L, 5L, 6L, 5L, 4L, 

[R] Why stacking rasters return NAs?

2013-03-18 Thread Jonsson
I have several rasters that I want to do some calculations ,basically
calculating the moving average. 

dir2 - list.files(D:\\2010+2011, *.bin, full.names = TRUE)
   saf=stack(dir2)
  movi -  overlay(stack(saf),fun=function(x) movingFun(x, fun=mean,
n=3,  na.rm=TRUE))
 Error in .overlayList(x, fun = fun, filename = filename, ...) : 
 cannot use this formula, probably because it is not vectorized
 I then checked the data but found that all values were returnd as NA and
this may explain why i am getting the error. 

saf
   class   : RasterStack 
   dimensions  : 720, 1440, 1036800, 601  (nrow, ncol, ncell, nlayers)
   resolution  : 0.25, 0.25  (x, y)
   extent  : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
  coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0
+no_defs 
 names   : Vol_025_H//00_1_wgs84, Vol_025_H//00_1_wgs84,
Vol_025_H//00_1_wgs84,Vol_025_H//00_1_wgs84, Vol_025_H//00_1_wgs84,
Vol_025_H//00_1_wgs84, , ... 
min values  :NA,NA, 
  
NA,   NA,NA,NA, 
  
NA, NA,NA,NA,   

NA,NA,NA,NA,
   
NA, ... 
max values  :NA,NA, 
  
NA,NA,NA,NA,
   
NA,NA,NA,NA,
   
NA,   NA,NA,NA, 
  
NA, ... 


I wonder why this is happening, I checked the files separably(summary) and
everything was right!as you can see bellow:

 ol_025_H14_2011092000_1_wgs84 Vol_025_H14_2011092100_1_wgs84  
Vol_025_H14_2011092200_1_wgs84 Vol_025_H14_2011092300_1_wgs84
Vol_025_H14_2011092400_1_wgs84
  Min.   0.0  0.000 
  
0.000  0.000  0.000
 1st Qu.0.31883  0.3163167  
   
0.3146436  0.3113111  0.3064551
 Median  .0   .000  

.000   .000   .000
 3rd Qu. .0   .000  

.000   .000   .000
  Max..0   .000 
 
.000   .000   .000
 NA's   0.0  0.000  
 
0.000  0.0

I am gratful to anyhelp



--
View this message in context: 
http://r.789695.n4.nabble.com/Why-stacking-rasters-return-NAs-tp4661706.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] Loop or some other way to parse by data generated values when it is not linear

2013-03-18 Thread plessthanpointohfive
I'm sorry for the really vague subject line but I am not sure how to 
succinctly describe what I am doing and what the problem is.


But, here goes:

1.  I have data with two-way data with frequencies.  Below is an 
example, though in reality I am looking at about 10 different variables 
that I am crossing so the values of X1 and X2 change.  X1 and X2 are 
place holders.


Here's the dataset (though using this first part does not happen in 
reality):


X1 - matrix(c(0, 1, 2, 3, 4, 99), nrow=18, ncol=1, byrow=T)
X2 - sort(matrix(c(0, 2, 4), nrow=18, ncol=1, byrow=T), decreasing=F)
Y - matrix(c(83, 107, 47, 27, 38, 1, 12 ,25, 14, 4, 9, 0, 14, 27, 28, 
13, 18, 0), nrow=18, ncol=1, byrow=T)

tmp.n - data.frame(X1, X2, Y)

The final data frame is what I actually get:


   X1 X2   Y
1   0  0  83
2   1  0 107
3   2  0  47
4   3  0  27
5   4  0  38
6  99  0   1
7   0  2  12
8   1  2  25
9   2  2  14
10  3  2   4
11  4  2   9
12 99  2   0
13  0  4  14
14  1  4  27
15  2  4  28
16  3  4  13
17  4  4  18
18 99  4   0


2.  What I want is:


 0   2   4
0 83 12 14
1   107 25 27
2 47 14 28
3 27   4 13
4 38   9 18
99 1   0   0


3.  I've been trying to do it using this (which is inside a function so 
I can vary what variables X1 and X2 are):



X1 - table(tmp.n[,1])
X2 - table(tmp.n[,2])

# Create the tmp.n.# datasets that contain the Y's.  Do this in a loop 
to automate

dta - NULL
for (i in 0:length(X1)) {
assign(tmp.n_, tmp.n[tmp.n[,1] == i, c(1,3)])
tmp.n_ - data.frame(tmp.n_[,2])
dta[i] - assign(paste(tmp.n., i, sep=), tmp.n_)
dta
}
dta2 - (data.frame(matrix(unlist(dta), nrow=n2[1], byrow=T)))
colnames(dta2) - names(X2)
dta2


And that works so long as X1 and X2 are linear.  In other words, if X1 
- seq(0, 4, 1).  But that 99 throws the whole thing off and it gives me 
this:


   X1 X2
1 107 25
2  27 47
3  14 28
4  27  4
5  13 38
6   9 18

It's basically breaks the whole thing.

I've not been able to figure this out and I've been like a dog with a 
bone trying to make it work with modifications to the for loop.  I know 
there is an easier way to do this, but my brain is no longer capable of 
thinking outside the box I've put it in.  So, I am turning to you for help.


Best,

Jen

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 with NA

2013-03-18 Thread Berend Hasselman

On 18-03-2013, at 16:49, Pete freeri...@gmail.com wrote:

 
 I have this little data.frame
 
 http://dl.dropbox.com/u/102669/nanotna.rdata
 
 Two column contains NA, so the best thing to do is use na.locf function (with
 fromLast = T)
 
 But locf function doesn't work because NA in my data.frame are not recognized 
 as
 real NA.
 
 Is there a way to substitute fake NA with real NA? In this case na.locf 
 function
 should work
 

Your data are all characters. Do

str(db)

to see that. What is probably supposed to be numeric is also character,
Somehow you have managed to read in data that R thinks is all chr.
Your NA are NA in reality: a character string NA.

You will have to review the method you used to get the data into R.
And make sure that what you want to be numeric is indeed numeric.
Then you can start to think about doing something about the NA's.

Berend

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 or some other way to parse by data generated values when it is not linear

2013-03-18 Thread arun
Hi,


library(reshape2)
 dcast(tmp.n,X1~X2,value.var=Y)
 # X1   0  2  4
#1  0  83 12 14
#2  1 107 25 27
#3  2  47 14 28
#4  3  27  4 13
#5  4  38  9 18
#6 99   1  0  0
A.K.


- Original Message -
From: plessthanpointohf...@gmail.com plessthanpointohf...@gmail.com
To: r-help@r-project.org
Cc: 
Sent: Monday, March 18, 2013 1:44 PM
Subject: [R] Loop or some other way to parse by data generated values when it 
is not linear

I'm sorry for the really vague subject line but I am not sure how to succinctly 
describe what I am doing and what the problem is.

But, here goes:

1.  I have data with two-way data with frequencies.  Below is an example, 
though in reality I am looking at about 10 different variables that I am 
crossing so the values of X1 and X2 change.  X1 and X2 are place holders.

Here's the dataset (though using this first part does not happen in reality):

X1 - matrix(c(0, 1, 2, 3, 4, 99), nrow=18, ncol=1, byrow=T)
X2 - sort(matrix(c(0, 2, 4), nrow=18, ncol=1, byrow=T), decreasing=F)
Y - matrix(c(83, 107, 47, 27, 38, 1, 12 ,25, 14, 4, 9, 0, 14, 27, 28, 13, 18, 
0), nrow=18, ncol=1, byrow=T)
tmp.n - data.frame(X1, X2, Y)

The final data frame is what I actually get:


   X1 X2   Y
1   0  0  83
2   1  0 107
3   2  0  47
4   3  0  27
5   4  0  38
6  99  0   1
7   0  2  12
8   1  2  25
9   2  2  14
10  3  2   4
11  4  2   9
12 99  2   0
13  0  4  14
14  1  4  27
15  2  4  28
16  3  4  13
17  4  4  18
18 99  4   0


2.  What I want is:


         0   2   4
0     83 12 14
1   107 25 27
2     47 14 28
3     27   4 13
4     38   9 18
99     1   0   0


3.  I've been trying to do it using this (which is inside a function so I can 
vary what variables X1 and X2 are):


X1 - table(tmp.n[,1])
X2 - table(tmp.n[,2])

# Create the tmp.n.# datasets that contain the Y's.  Do this in a loop to 
automate
dta - NULL
for (i in 0:length(X1)) {
assign(tmp.n_, tmp.n[tmp.n[,1] == i, c(1,3)])
tmp.n_ - data.frame(tmp.n_[,2])
dta[i] - assign(paste(tmp.n., i, sep=), tmp.n_)
dta
}
dta2 - (data.frame(matrix(unlist(dta), nrow=n2[1], byrow=T)))
colnames(dta2) - names(X2)
dta2


And that works so long as X1 and X2 are linear.  In other words, if X1 - 
seq(0, 4, 1).  But that 99 throws the whole thing off and it gives me this:

   X1 X2
1 107 25
2  27 47
3  14 28
4  27  4
5  13 38
6   9 18

It's basically breaks the whole thing.

I've not been able to figure this out and I've been like a dog with a bone 
trying to make it work with modifications to the for loop.  I know there is an 
easier way to do this, but my brain is no longer capable of thinking outside 
the box I've put it in.  So, I am turning to you for help.

Best,

Jen

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Loop or some other way to parse by data generated values when it is not linear

2013-03-18 Thread Marc Schwartz
On Mar 18, 2013, at 12:44 PM, plessthanpointohf...@gmail.com wrote:

 I'm sorry for the really vague subject line but I am not sure how to 
 succinctly describe what I am doing and what the problem is.
 
 But, here goes:
 
 1.  I have data with two-way data with frequencies.  Below is an example, 
 though in reality I am looking at about 10 different variables that I am 
 crossing so the values of X1 and X2 change.  X1 and X2 are place holders.
 
 Here's the dataset (though using this first part does not happen in reality):
 
 X1 - matrix(c(0, 1, 2, 3, 4, 99), nrow=18, ncol=1, byrow=T)
 X2 - sort(matrix(c(0, 2, 4), nrow=18, ncol=1, byrow=T), decreasing=F)
 Y - matrix(c(83, 107, 47, 27, 38, 1, 12 ,25, 14, 4, 9, 0, 14, 27, 28, 13, 
 18, 0), nrow=18, ncol=1, byrow=T)
 tmp.n - data.frame(X1, X2, Y)
 
 The final data frame is what I actually get:
 
 
   X1 X2   Y
 1   0  0  83
 2   1  0 107
 3   2  0  47
 4   3  0  27
 5   4  0  38
 6  99  0   1
 7   0  2  12
 8   1  2  25
 9   2  2  14
 10  3  2   4
 11  4  2   9
 12 99  2   0
 13  0  4  14
 14  1  4  27
 15  2  4  28
 16  3  4  13
 17  4  4  18
 18 99  4   0
 
 
 2.  What I want is:
 
 
 0   2   4
 0 83 12 14
 1   107 25 27
 2 47 14 28
 3 27   4 13
 4 38   9 18
 99 1   0   0
 
 
 3.  I've been trying to do it using this (which is inside a function so I can 
 vary what variables X1 and X2 are):
 
 
 X1 - table(tmp.n[,1])
 X2 - table(tmp.n[,2])
 
 # Create the tmp.n.# datasets that contain the Y's.  Do this in a loop to 
 automate
 dta - NULL
 for (i in 0:length(X1)) {
 assign(tmp.n_, tmp.n[tmp.n[,1] == i, c(1,3)])
 tmp.n_ - data.frame(tmp.n_[,2])
 dta[i] - assign(paste(tmp.n., i, sep=), tmp.n_)
 dta
 }
 dta2 - (data.frame(matrix(unlist(dta), nrow=n2[1], byrow=T)))
 colnames(dta2) - names(X2)
 dta2
 
 
 And that works so long as X1 and X2 are linear.  In other words, if X1 - 
 seq(0, 4, 1).  But that 99 throws the whole thing off and it gives me this:
 
   X1 X2
 1 107 25
 2  27 47
 3  14 28
 4  27  4
 5  13 38
 6   9 18
 
 It's basically breaks the whole thing.
 
 I've not been able to figure this out and I've been like a dog with a bone 
 trying to make it work with modifications to the for loop.  I know there is 
 an easier way to do this, but my brain is no longer capable of thinking 
 outside the box I've put it in.  So, I am turning to you for help.
 
 Best,
 
 Jen


Something like this?

 xtabs(Y ~ X1 + X2, data = tmp.n)
X2
X1 0   2   4
  0   83  12  14
  1  107  25  27
  2   47  14  28
  3   27   4  13
  4   38   9  18
  99   1   0   0

See ?xtabs

Regards,

Marc Schwartz

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Confirmatory factor analysis using the sem package. TLI CFI and RMSEA absent from model summary.

2013-03-18 Thread John Fox
Dear Kevin,

In earlier versions of the sem package, the summary() method for model objects 
produced many (and from version to version, an increasing number of) fit 
indices by default. For various reasons, I decided to provide only the LR 
chisquare test for the model, the BIC, and the AIC by default, and to leave it 
up the user to decide what fit indices to show. This can be controlled with the 
fit.indices option as well as the fit.indices argument to the summary() method.

Best,
 John

On Mon, 18 Mar 2013 16:46:03 +
 Kevin Cheung k.che...@derby.ac.uk wrote:
 Dear John,
 
 Thank you for taking the time to help me with this, I have been able to 
 obtain fit indices using the information that you provided.
 
 Note to users searching archived R-help posts about this issue: The 
 instructional video I looked at (http://vimeo.com/38941937) gave fit indices 
 using the default summary() command without any additional arguments. This 
 may have been due to a different version of R (I noticed that the instructor 
 was using a mac based OS).
 
 With regards,
 Kevin
 
 Kevin Yet Fong Cheung, Bsc., MRes., MBPsS.
 Postgraduate Researcher
 Centre for Psychological Research
 University of Derby
 Kedleston Road
 Derby DE22 1GB
 k.che...@derby.ac.ukmailto:k.che...@derby.ac.uk
 01332592081
 
 http://derby.academia.edu/KevinCheung
 
 
 -Original Message-
 From: John Fox [mailto:j...@mcmaster.ca]
 Sent: 18 March 2013 15:55
 To: Kevin Cheung
 Cc: r-help@r-project.org
 Subject: Re: [R] Confirmatory factor analysis using the sem package. TLI CFI 
 and RMSEA absent from model summary.
 
 Dear Kevin,
 
 See ?summary.objectiveML, and in particular the description of the 
 fit.indices argument. By default, the summary() method doesn't print many fit 
 indices, but many are available optionally.
 
 I hope this helps,
  John
 
 
 John Fox
 Sen. William McMaster Prof. of Social Statistics Department of Sociology 
 McMaster University Hamilton, Ontario, Canada http://socserv.mcmaster.ca/jfox/
 
 On Mon, 18 Mar 2013 15:00:06 +
  Kevin Cheung k.che...@derby.ac.uk wrote:
  Hi R-help,
 
  I am using the sem package to run confirmatory factor analysis (cfa) on 
  some questionnaire data collected from 307 participants. I have been 
  running R-2.15.3 in Windows in conjunction with R studio. The model I am 
  using was developed from exploratory factor analysis of a separate dataset 
  (n=439); it includes 18 items that load onto 3 factors. I have used the sem 
  package documentation and this video (http://vimeo.com/38941937) to run the 
  cfa and obtain a chi-square statistic for the model. However, when I use 
  the summary() function, the model does not provide indices of fit; I need 
  these as part of my analysis output. In particular, I am looking for the 
  Tucker Lewis Index (TLI), Comparative Fit Index (CFI),  the Root Mean 
  Square of Approximation (RMSEA). I have checked the documentation and 
  cannot seem to find any reason for this; none of the arguments listed with 
  the sem command indicate that I have to specify these as part of the 
  output. In addition, the analysis examp!
 le!
   f!
   rom the video includes these indices as part of the output, but my 
  analysis does not provide them. I have included my code with comments below:
 
  
 
  library(sem)
 
  validation.data -
  structure(list(V1 = c(5L, 4L, 2L, 4L, 5L, 6L, 6L, 4L, 5L, 3L, 6L, 5L,
  4L, 5L, 3L, 5L, 5L, 5L, 5L, 5L, 5L, 4L, 5L, 5L, 5L, 4L, 5L, 4L, 4L,
  5L, 4L, 4L, 5L, 5L, 5L, 5L, 2L, 6L, 5L, 6L, 4L, 5L, 6L, 5L, 5L, 4L,
  5L, 5L, 3L, 5L, 5L, 5L, 5L, 4L, 2L, 5L, 5L, 5L, 4L, 6L, 4L, 6L, 5L,
  5L, 5L, 4L, 5L, 5L, 4L, 5L, 6L, 4L, 5L, 4L, 5L, 5L, 5L, 3L, 5L, 5L,
  3L, 5L, 4L, 5L, 2L, 6L, 4L, 4L, 4L, 5L, 5L, 4L, 6L, 2L, 4L, 5L, 4L,
  5L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 4L, 2L, 4L, 4L, 5L, 5L, 4L,
  4L, 5L, 5L, 5L, 5L, 5L, 5L, 4L, 5L, 5L, 5L, 5L, 5L, 6L, 5L, 5L, 4L,
  4L, 5L, 4L, 4L, 5L, 4L, 4L, 4L, 5L, 5L, 6L, 5L, 5L, 5L, 4L, 4L, 3L,
  4L, 5L, 5L, 5L, 2L, 5L, 5L, 5L, 4L, 5L, 5L, 5L, 5L, 6L, 4L, 5L, 5L,
  6L, 5L, 5L, 5L, 5L, 4L, 5L, 4L, 5L, 5L, 4L, 4L, 4L, 4L, 5L, 2L, 4L,
  4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 4L, 4L, 5L, 4L, 5L, 5L, 4L, 5L, 5L,
  4L, 4L, 5L, 5L, 5L, 5L, 4L, 4L, 2L, 6L, 5L, 5L, 5L, 6L, 4L, 3L, 5L,
  5L, 5L, 5L, 4L, 4L, 6L, 6L, 5L, 5L, 5L, 3L, 6L, 5L, 5L, 5L, 4L, 5L,
  5L, 4L, 5L, 4L, 6L, 5L, 4L, 4L, 5L, 4L, 3L, 4L, 5L, 4L, 5L, 6L, 2L,
  4L, 4L, 5L, 4L, 4L, 4L, 5L, 6L, 5L, 4L, 4L, 4L, 6L, 6L, 4L, 5L, 5L,
  5L, 2L, 4L, 4L, 5L, 4L, 5L, 5L, 4L, 5L, 3L, 5L, 6L, 5L, 4L, 4L, 3L,
  3L, 4L, 5L, 5L, 1L, 4L, 5L, 3L, 5L, 1L, 6L, 5L, 4L, 4L, 5L, 5L, 4L,
  5L, 5L, 6L, 5L, 5L, 5L), V2 = c(5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 5L,
  5L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 6L, 5L, 6L, 6L, 4L, 5L, 6L, 6L, 6L,
  6L, 6L, 6L, 6L, 6L, 6L, 6L, 5L, 6L, 6L, 5L, 5L, 5L, 6L, 6L, 6L, 6L,
  6L, 6L, 6L, 6L, 5L, 5L, 6L, 5L, 6L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L,
  4L, 6L, 5L, 6L, 5L, 5L, 

[R] Lattice strip with custom shingle.interval

2013-03-18 Thread Alex van der Spek
Trying to make a strip in lattice xyplots to display shingle intervals
with a sensible number of digits:


First I make a 2x2 matrix with rounded values of the shingle itnervals:

lec - round(matrix(unlist(levels(equal.count(sd$Frother))), ncol = 2,
byrow = TRUE), 3)

Next I make a custom strip function to handle the difference between a
factor proper and a shingle:

my.strip - function(which.given, which.panel, ..., shingle.intervals) {
if (which.given == 1) {
strip.default(which.given, which.panel, ...,strip.name 
= FALSE,
strip.levels = TRUE, shingle.intervals = lec)
} else {
strip.default(which.given, which.panel, ..., strip.name 
= FALSE,
strip.levels = TRUE)
}
}

Then the call to xyplot.

xyplot(Grade ~ Recovery| equal.count(Frother) + Row,
   groups = interaction(Row1.Mode, Row3.Mode, drop = TRUE),
   data = sd, pch = 1, layout = c(6,2),
   main = 'By Row and Frother dosing',
 strip = my.strip,
   sub = 'Colour indicates control mode',
 auto.key = list(columns = 5))


This however still gives full precision levels.

Any help much appreciated.

Alex van der Spek

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 or some other way to parse by data generated values when it is not linear

2013-03-18 Thread plessthanpointohfive

That works great!  I didn't even know there was a reshape2 package.

Thanks!

Jen


On 03/18/2013 01:49 PM, arun wrote:

Hi,


library(reshape2)
  dcast(tmp.n,X1~X2,value.var=Y)
  # X1   0  2  4
#1  0  83 12 14
#2  1 107 25 27
#3  2  47 14 28
#4  3  27  4 13
#5  4  38  9 18
#6 99   1  0  0
A.K.


- Original Message -
From: plessthanpointohf...@gmail.com plessthanpointohf...@gmail.com
To: r-help@r-project.org
Cc:
Sent: Monday, March 18, 2013 1:44 PM
Subject: [R] Loop or some other way to parse by data generated values when it 
is not linear

I'm sorry for the really vague subject line but I am not sure how to succinctly 
describe what I am doing and what the problem is.

But, here goes:

1.  I have data with two-way data with frequencies.  Below is an example, 
though in reality I am looking at about 10 different variables that I am 
crossing so the values of X1 and X2 change.  X1 and X2 are place holders.

Here's the dataset (though using this first part does not happen in reality):

X1 - matrix(c(0, 1, 2, 3, 4, 99), nrow=18, ncol=1, byrow=T)
X2 - sort(matrix(c(0, 2, 4), nrow=18, ncol=1, byrow=T), decreasing=F)
Y - matrix(c(83, 107, 47, 27, 38, 1, 12 ,25, 14, 4, 9, 0, 14, 27, 28, 13, 18, 
0), nrow=18, ncol=1, byrow=T)
tmp.n - data.frame(X1, X2, Y)

The final data frame is what I actually get:


X1 X2   Y
1   0  0  83
2   1  0 107
3   2  0  47
4   3  0  27
5   4  0  38
6  99  0   1
7   0  2  12
8   1  2  25
9   2  2  14
10  3  2   4
11  4  2   9
12 99  2   0
13  0  4  14
14  1  4  27
15  2  4  28
16  3  4  13
17  4  4  18
18 99  4   0


2.  What I want is:


  0   2   4
0 83 12 14
1   107 25 27
2 47 14 28
3 27   4 13
4 38   9 18
99 1   0   0


3.  I've been trying to do it using this (which is inside a function so I can 
vary what variables X1 and X2 are):


X1 - table(tmp.n[,1])
X2 - table(tmp.n[,2])

# Create the tmp.n.# datasets that contain the Y's.  Do this in a loop to 
automate
dta - NULL
for (i in 0:length(X1)) {
assign(tmp.n_, tmp.n[tmp.n[,1] == i, c(1,3)])
tmp.n_ - data.frame(tmp.n_[,2])
dta[i] - assign(paste(tmp.n., i, sep=), tmp.n_)
dta
}
dta2 - (data.frame(matrix(unlist(dta), nrow=n2[1], byrow=T)))
colnames(dta2) - names(X2)
dta2


And that works so long as X1 and X2 are linear.  In other words, if X1 - 
seq(0, 4, 1).  But that 99 throws the whole thing off and it gives me this:

X1 X2
1 107 25
2  27 47
3  14 28
4  27  4
5  13 38
6   9 18

It's basically breaks the whole thing.

I've not been able to figure this out and I've been like a dog with a bone 
trying to make it work with modifications to the for loop.  I know there is an 
easier way to do this, but my brain is no longer capable of thinking outside 
the box I've put it in.  So, I am turning to you for help.

Best,

Jen

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] OrgMassSpecR peak area issue

2013-03-18 Thread Christopher Beaver
Hello!

I am having an issue with the OrgMassSpecR package.  I run my HPLC using a
DAD detector.  My raw data is exported form chemstation as a csv file.  I
then upload the csv into Rstudio no problem.  Using the DrawChromatogram
function, I get a nice chromatogram, and my retention time, peak area, and
apex intensity values are given as well.

The problem comes with the peak area value given. The peak area is much
smaller than a value that would make sense.  My peak area value is actually
less than my apex intensity value.  Is this because I am using a DAD
detector rather than an MS? If so, is there a simply way to edit the peak
area equation so that it will also work with absorbance values?

Any help is greatly appreciated.

Thanks for your time.

Chris Beaver

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Loop or some other way to parse by data generated values when it is not linear

2013-03-18 Thread plessthanpointohfive

This works really great, too.

Thanks for this option.

I am glad I have two ways to accomplish this.

I KNEW it was going to be that simple. I was convinced I was making more 
out of it than I needed to and I'm glad to see I was right.


Best,

Jen

(It's fun to be new at something when I'm this old).







On 03/18/2013 01:50 PM, Marc Schwartz wrote:

On Mar 18, 2013, at 12:44 PM, plessthanpointohf...@gmail.com wrote:


I'm sorry for the really vague subject line but I am not sure how to succinctly 
describe what I am doing and what the problem is.

But, here goes:

1.  I have data with two-way data with frequencies.  Below is an example, 
though in reality I am looking at about 10 different variables that I am 
crossing so the values of X1 and X2 change.  X1 and X2 are place holders.

Here's the dataset (though using this first part does not happen in reality):

X1 - matrix(c(0, 1, 2, 3, 4, 99), nrow=18, ncol=1, byrow=T)
X2 - sort(matrix(c(0, 2, 4), nrow=18, ncol=1, byrow=T), decreasing=F)
Y - matrix(c(83, 107, 47, 27, 38, 1, 12 ,25, 14, 4, 9, 0, 14, 27, 28, 13, 18, 
0), nrow=18, ncol=1, byrow=T)
tmp.n - data.frame(X1, X2, Y)

The final data frame is what I actually get:


   X1 X2   Y
1   0  0  83
2   1  0 107
3   2  0  47
4   3  0  27
5   4  0  38
6  99  0   1
7   0  2  12
8   1  2  25
9   2  2  14
10  3  2   4
11  4  2   9
12 99  2   0
13  0  4  14
14  1  4  27
15  2  4  28
16  3  4  13
17  4  4  18
18 99  4   0


2.  What I want is:


 0   2   4
0 83 12 14
1   107 25 27
2 47 14 28
3 27   4 13
4 38   9 18
99 1   0   0


3.  I've been trying to do it using this (which is inside a function so I can 
vary what variables X1 and X2 are):


X1 - table(tmp.n[,1])
X2 - table(tmp.n[,2])

# Create the tmp.n.# datasets that contain the Y's.  Do this in a loop to 
automate
dta - NULL
for (i in 0:length(X1)) {
assign(tmp.n_, tmp.n[tmp.n[,1] == i, c(1,3)])
tmp.n_ - data.frame(tmp.n_[,2])
dta[i] - assign(paste(tmp.n., i, sep=), tmp.n_)
dta
}
dta2 - (data.frame(matrix(unlist(dta), nrow=n2[1], byrow=T)))
colnames(dta2) - names(X2)
dta2


And that works so long as X1 and X2 are linear.  In other words, if X1 - 
seq(0, 4, 1).  But that 99 throws the whole thing off and it gives me this:

   X1 X2
1 107 25
2  27 47
3  14 28
4  27  4
5  13 38
6   9 18

It's basically breaks the whole thing.

I've not been able to figure this out and I've been like a dog with a bone 
trying to make it work with modifications to the for loop.  I know there is an 
easier way to do this, but my brain is no longer capable of thinking outside 
the box I've put it in.  So, I am turning to you for help.

Best,

Jen


Something like this?


xtabs(Y ~ X1 + X2, data = tmp.n)

 X2
X1 0   2   4
   0   83  12  14
   1  107  25  27
   2   47  14  28
   3   27   4  13
   4   38   9  18
   99   1   0   0

See ?xtabs

Regards,

Marc Schwartz



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] How to stop set.seed() besides exiting out of R?

2013-03-18 Thread C W
Hi list,

I am curious how to stop the set.seed(), I don't want the same repeated
random number.  I know I can set it to a different seed, but I don't want
to go through the trouble of setting different seed every time.

Thanks,
Mike

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] new question

2013-03-18 Thread arun



 z.boxplot- function(lst){
 new.list-  lapply(lst,function(x) x[x$FDR0.01,])
print(new.list)
  par(mfrow=c(2,2))
b1-lapply(names(new.list),function(x) lapply(new.list[x],function(y) 
boxplot(FDR~z,data=y,xlab=Charge,ylab=FDR,main=x)))

}
 z.boxplot(ListFacGroup) #prints new.list

If you want to turn off that:
 z.boxplot- function(lst){
 new.list-  lapply(lst,function(x) x[x$FDR0.01,])
#print(new.list)
  par(mfrow=c(2,2))
b1-lapply(names(new.list),function(x) lapply(new.list[x],function(y) 
boxplot(FDR~z,data=y,xlab=Charge,ylab=FDR,main=x)))

}
 z.boxplot(ListFacGroup)
A.K.






From: Vera Costa veracosta...@gmail.com
To: arun smartpink...@yahoo.com 
Sent: Monday, March 18, 2013 1:59 PM
Subject: Re: new question


For example, if I run you code without pdf and dev.off I have what I 
want 

directory- C:/Users/Vera Costa/Desktop/dados.lixo
 #modified the function
GetFileList - function(directory,number){
  setwd(directory)
  filelist1-dir()

    lista-dir(directory,pattern = paste(MSMS_,number,PepInfo.txt,sep=), 
full.names = TRUE, recursive = TRUE)
  output- list(filelist1,lista)
  return(output)
 }
file.list.names-GetFileList(directory,23)[[1]]
lista-GetFileList(directory,23)[[2]]
FacGroup-c(0,1,0,2,2,0,3)
ReadDir-function(FacGroup){
  list.new-lista[FacGroup!=0]
  read.list-lapply(list.new, function(x) read.table(x,header=TRUE, sep = \t))
  names(read.list)-file.list.names[FacGroup!=0]
  return (read.list)
 }
ListFacGroup-ReadDir(FacGroup)
ListFacGroup
 z.boxplot- function(lst){
 new.list-  lapply(lst,function(x) x[x$FDR0.01,])
 print(new.list)
 #pdf(VeraBP.pdf)
 par(mfrow=c(2,2))
 lapply(names(new.list),function(x) lapply(new.list[x],function(y) 
boxplot(FDR~z,data=y,xlab=Charge,ylab=FDR,main=x)))
 #dev.off()
 }
 z.boxplot(ListFacGroup)





But I have the results too (I don't need it)


[[1]]
[[1]]$a2
[[1]]$a2$stats
 [,1] [,2]
[1,] 0.00 0.00
[2,] 0.00 0.00
[3,] 0.0001355197 0.0002175095
[4,] 0.0010588777 0.0004350190
[5,] 0.0016571381 0.0004350190
[[1]]$a2$n
[1] 8 2
[[1]]$a2$conf
  [,1]  [,2]
[1,] -0.0004559846 -0.0002685062
[2,]  0.0007270240  0.0007035253
[[1]]$a2$out
[1] 0.00494012
[[1]]$a2$group
[1] 1
[[1]]$a2$names
[1] 2 3

[[2]]
[[2]]$c2
[[2]]$c2$stats
 [,1] [,2]
[1,] 0.00 0.00
[2,] 0.00 0.00
[3,] 0.0001355197 0.0002175095
[4,] 0.0010588777 0.0004350190
[5,] 0.0016571381 0.0004350190
[[2]]$c2$n
[1] 8 2
[[2]]$c2$conf
  [,1]  [,2]
[1,] -0.0004559846 -0.0002685062
[2,]  0.0007270240  0.0007035253
[[2]]$c2$out
[1] 0.00494012
[[2]]$c2$group
[1] 1
[[2]]$c2$names
[1] 2 3

[[3]]
[[3]]$c3
[[3]]$c3$stats
    [,1] [,2]
[1,] 0.0 0.00
[2,] 0.0 0.00
[3,] 0.0 0.00
[4,] 0.002226409 0.0002086594
[5,] 0.002226409 0.0004173187
[[3]]$c3$n
[1] 6 4
[[3]]$c3$conf
 [,1]  [,2]
[1,] -0.001436105 -0.0001648409
[2,]  0.001436105  0.0001648409
[[3]]$c3$out
[1] 0.00560348
[[3]]$c3$group
[1] 1
[[3]]$c3$names
[1] 2 3

[[4]]
[[4]]$t2
[[4]]$t2$stats
 [,1] [,2] [,3]
[1,]    0 0.00    0
[2,]    0 0.00    0
[3,]    0 0.0002908668    0
[4,]    0 0.0025929827    0
[5,]    0 0.005251    0
[[4]]$t2$n
[1]  1 10  5
[[4]]$t2$conf
 [,1] [,2] [,3]
[1,]    0 -0.001004691    0
[2,]    0  0.001586424    0
[[4]]$t2$out
[1] 0.0092051934 0.0007174888
[[4]]$t2$group
[1] 2 3
[[4]]$t2$names
[1] 1 2 3

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 stop set.seed() besides exiting out of R?

2013-03-18 Thread Nordlund, Dan (DSHS/RDA)
 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
 project.org] On Behalf Of C W
 Sent: Monday, March 18, 2013 11:27 AM
 To: r-help
 Subject: [R] How to stop set.seed() besides exiting out of R?
 
 Hi list,
 
 I am curious how to stop the set.seed(), I don't want the same repeated
 random number.  I know I can set it to a different seed, but I don't
 want
 to go through the trouble of setting different seed every time.
 
 Thanks,
 Mike
 

Can you show us how you are using set.seed() that results in getting the same 
sequence repeatedly?  If you are doing simulations in a loop, then set the seed 
once, outside the loop.  Otherwise, I am not sure what you are doing that 
causes problems.  A reproducible example would really help.

Dan

Daniel J. Nordlund
Washington State Department of Social and Health Services
Planning, Performance, and Accountability
Research and Data Analysis Division
Olympia, WA 98504-5204

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 stop set.seed() besides exiting out of R?

2013-03-18 Thread C W
set.seed(100)
for (i in 1:100){
a - rnorm(1000, mean=0, sd=1)
hist(a)
}

#Now say, I want to simulate without being confined to set.seed(100), I
just want to get a simulation like how R normally does it.

Mike

On Mon, Mar 18, 2013 at 2:40 PM, Nordlund, Dan (DSHS/RDA) 
nord...@dshs.wa.gov wrote:

  -Original Message-
  From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
  project.org] On Behalf Of C W
  Sent: Monday, March 18, 2013 11:27 AM
  To: r-help
  Subject: [R] How to stop set.seed() besides exiting out of R?
 
  Hi list,
 
  I am curious how to stop the set.seed(), I don't want the same repeated
  random number.  I know I can set it to a different seed, but I don't
  want
  to go through the trouble of setting different seed every time.
 
  Thanks,
  Mike
 

 Can you show us how you are using set.seed() that results in getting the
 same sequence repeatedly?  If you are doing simulations in a loop, then set
 the seed once, outside the loop.  Otherwise, I am not sure what you are
 doing that causes problems.  A reproducible example would really help.

 Dan

 Daniel J. Nordlund
 Washington State Department of Social and Health Services
 Planning, Performance, and Accountability
 Research and Data Analysis Division
 Olympia, WA 98504-5204

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] OrgMassSpecR peak area issue

2013-03-18 Thread Bryan Hanson
If you type DrawChromatogram you can see the method used to calculate the peak 
area.  Looks to me like you could easily hack it if you wanted.  The relevant 
part about peak areas is this:

for (j in 1:n) {
k - (j%%n) + 1
x[j] - peakTime[j] * peakIntensity[k] - peakTime[k] * 
peakIntensity[j]
}
peakArea[i] - abs(sum(x)/2)

which looks pretty standard to me, though I'm not clear right off the top of my 
head why they are dividing by 2.  You can always contact the maintainer.

Bryan

On Mar 18, 2013, at 1:34 PM, Christopher Beaver christopher.bea...@gmail.com 
wrote:

 Hello!
 
 I am having an issue with the OrgMassSpecR package.  I run my HPLC using a
 DAD detector.  My raw data is exported form chemstation as a csv file.  I
 then upload the csv into Rstudio no problem.  Using the DrawChromatogram
 function, I get a nice chromatogram, and my retention time, peak area, and
 apex intensity values are given as well.
 
 The problem comes with the peak area value given. The peak area is much
 smaller than a value that would make sense.  My peak area value is actually
 less than my apex intensity value.  Is this because I am using a DAD
 detector rather than an MS? If so, is there a simply way to edit the peak
 area equation so that it will also work with absorbance values?
 
 Any help is greatly appreciated.
 
 Thanks for your time.
 
 Chris Beaver
 
   [[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] OrgMassSpecR peak area issue

2013-03-18 Thread Claudia Beleites
Hi Chris,

 I am having an issue with the OrgMassSpecR package.  I run my HPLC
 using a DAD detector. 
You are on a statistics IDE related mailing list. Have mercy with
people from other fields and tell them that you are using a diode array
to measure UV/VIS absorption. (And possibly let them know that you
expect the absorbance A = lg I_0 - lg I ~ c.)

 My raw data is exported form chemstation as a
 csv file.  I then upload the csv into Rstudio no problem.  Using the
 DrawChromatogram function, I get a nice chromatogram, and my
 retention time, peak area, and apex intensity values are given as
 well.
 
 The problem comes with the peak area value given. The peak area is
 much smaller than a value that would make sense.
How do you know that (see next comment)?

 My peak area value is actually less than my apex intensity value. 
This is not a good criterion to determine what area value would
actually make sense: area and intensity have different units!

Possible solution: a glance on the code in DrawChromatogram reveals
that really the polygon area is calculated (as the manual specifies). 

Thus the area will be in counts*s or counts*min, and of course 1
count*min = 60 counts*s.  How long does your analyte
take to elute? Unless it is  2 min (if time is in min) or  2 s (for
time scale in s), the numeric value of the area should be  A_max
(approximating the peak as triangule).

Your apex (max) absorbance should ideally be a bit below 1, so a rough
guesstimate for peak area would be 1/2 A_max * Δt which will be quite
below 1 if you measure time in minutes.

If you detect by mass spec, you get ion counts which are large
numbers, so areas are likely to be  1 (regardless of min or s time
scale).


 Is this because I am using a DAD detector rather than an MS? If so,
 is there a simply way to edit the peak area equation so that it will
 also work with absorbance values?
Most probably you just want to get your units right!

Hope that helps, 

Claudia

PS: for future questions of this sort, you may want to consider asking
on stackoverflow.com (or chemistry.stackexchange.com) where you can
post nicely formatted code, calculation results and images with your
question.

-- 
Claudia Beleites
Spectroscopy/Imaging
Institute of Photonic Technology 
Albert-Einstein-Str. 9
07745 Jena
Germany

email: claudia.belei...@ipht-jena.de
phone: +49 3641 206-133
fax:   +49 2641 206-399

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 stop set.seed() besides exiting out of R?

2013-03-18 Thread Nordlund, Dan (DSHS/RDA)
As I understand it, how R “‘normally” does it is to use the system clock 
to set the seed once per session, unless you use set.seed() to set a new seed. 
You chose to set the seed to a different value.  But from that point on, the 
pseudo random number generation continues  in the same way it “normally” 
does.  In your code below, each of your 100 histograms will be different.  If 
you then execute the for loop again (but not the set.seed(100) statement), you 
will get a different set of histograms.  The only way you would be “confined 
to set.seed(100)” is if you keep resetting the seed to 100.

Dan

Daniel J. Nordlund
Washington State Department of Social and Health Services
Planning, Performance, and Accountability
Research and Data Analysis Division
Olympia, WA 98504-5204

From: C W [mailto:tmrs...@gmail.com]
Sent: Monday, March 18, 2013 11:50 AM
To: Nordlund, Dan (DSHS/RDA)
Cc: r-help
Subject: Re: [R] How to stop set.seed() besides exiting out of R?

set.seed(100)
for (i in 1:100){
a - rnorm(1000, mean=0, sd=1)
hist(a)
}

#Now say, I want to simulate without being confined to set.seed(100), I just 
want to get a simulation like how R normally does it.

Mike

On Mon, Mar 18, 2013 at 2:40 PM, Nordlund, Dan (DSHS/RDA) 
nord...@dshs.wa.govmailto:nord...@dshs.wa.gov wrote:
 -Original Message-
 From: r-help-boun...@r-project.orgmailto:r-help-boun...@r-project.org 
 [mailto:r-help-bounces@r-mailto:r-help-bounces@r-
 project.orghttp://project.org] On Behalf Of C W
 Sent: Monday, March 18, 2013 11:27 AM
 To: r-help
 Subject: [R] How to stop set.seed() besides exiting out of R?

 Hi list,

 I am curious how to stop the set.seed(), I don't want the same repeated
 random number.  I know I can set it to a different seed, but I don't
 want
 to go through the trouble of setting different seed every time.

 Thanks,
 Mike

Can you show us how you are using set.seed() that results in getting the same 
sequence repeatedly?  If you are doing simulations in a loop, then set the seed 
once, outside the loop.  Otherwise, I am not sure what you are doing that 
causes problems.  A reproducible example would really help.

Dan

Daniel J. Nordlund
Washington State Department of Social and Health Services
Planning, Performance, and Accountability
Research and Data Analysis Division
Olympia, WA 98504-5204

__
R-help@r-project.orgmailto:R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] data.frame with NA

2013-03-18 Thread Berend Hasselman

On 18-03-2013, at 19:57, Pietro freeri...@gmail.com wrote:

 Yes, it's true Berend!
 
 What i do is simply use read.xlsx  function
 
 db - read.xlsx2(c:/mydb.xlsx,1,as.data.frame=T)
 
 This is excel file i use:
 http://dl.dropbox.com/u/102669/mydb.xlsx
 

I don't have the required packages installed to read .xlsx files since I have 
no need.
But I can load your xlsx in LibreOffice without any issues.
I'm not sure what read.xslx2 actually does with columns. It seems that is 
regards the colClass of all columns as character.
So maybe you should try read.xlsx or export the file as a .csv
Or specify colClasses.

Berend


 I can't find  a way to import as numeric.
 My objective is to be able to work (in R) with my NA's
 
 
 At 18.46 18/03/2013, Berend Hasselman wrote:
 
 On 18-03-2013, at 16:49, Pete freeri...@gmail.com wrote:
 
 
  I have this little data.frame
 
  http://dl.dropbox.com/u/102669/nanotna.rdata
 
  Two column contains NA, so the best thing to do is use na.locf function 
  (with
  fromLast = T)
 
  But locf function doesn't work because NA in my data.frame are not 
  recognized as
  real NA.
 
  Is there a way to substitute fake NA with real NA? In this case na.locf 
  function
  should work
 
 
 Your data are all characters. Do
 
 str(db)
 
 to see that. What is probably supposed to be numeric is also character,
 Somehow you have managed to read in data that R thinks is all chr.
 Your NA are NA in reality: a character string NA.
 
 You will have to review the method you used to get the data into R.
 And make sure that what you want to be numeric is indeed numeric.
 Then you can start to think about doing something about the NA's.
 
 Berend
 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] ksvm

2013-03-18 Thread Li, Yan
Hi All,

I'm developing a ksvm application. I need to use the results obtained from ksvm 
in R to construct or restore the ksvm formula and further use this formula to 
predict new data.

The ksvm class contains xmatrix, ymatrix, coefficients, alpha, b and fitted 
value. Taking linear svm as an example ( vanilladot), the formula is like :

 F(x) = sign(sum(w*x*y*z)+b)

Where, x =  xmatrix, y = ymatrix, b is intercept and z is the input matrix? The 
w is alpha or coefficient from ksvm class? I did some calculation but the 
result is not the same as the fitted value.

Can you give me some idea on the prediction of ksvm by the formula?

Thanks a lot!!!

Regards,
Yan


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] How to stop set.seed() besides exiting out of R?

2013-03-18 Thread C W
Yes, I agree with you.  I guess what I was really looking for is a function
like UNset.seed()?

By having set.seed(), I can have reproducible code.  But what if I want to
check my work against what's produced from set.seed(100)?

I really want to escape from the shadow of set.seed(), can I unset it?

On Mon, Mar 18, 2013 at 3:07 PM, Nordlund, Dan (DSHS/RDA) 
nord...@dshs.wa.gov wrote:

 As I understand it, how R “‘normally” does it is to use the system clock
 to set the seed once per session, unless you use set.seed() to set a new
 seed. You chose to set the seed to a different value.  But from that point
 on, the pseudo random number generation continues  in the same way it
 “normally” does.  In your code below, each of your 100 histograms will be
 different.  If you then execute the for loop again (but not the
 set.seed(100) statement), you will get a different set of histograms.  The
 only way you would be “confined to set.seed(100)” is if you keep resetting
 the seed to 100.

 Dan

 Daniel J. Nordlund
 Washington State Department of Social and Health Services
 Planning, Performance, and Accountability
 Research and Data Analysis Division
 Olympia, WA 98504-5204

 From: C W [mailto:tmrs...@gmail.com]
 Sent: Monday, March 18, 2013 11:50 AM
 To: Nordlund, Dan (DSHS/RDA)
 Cc: r-help
 Subject: Re: [R] How to stop set.seed() besides exiting out of R?

 set.seed(100)
 for (i in 1:100){
 a - rnorm(1000, mean=0, sd=1)
 hist(a)
 }

 #Now say, I want to simulate without being confined to set.seed(100), I
 just want to get a simulation like how R normally does it.

 Mike

 On Mon, Mar 18, 2013 at 2:40 PM, Nordlund, Dan (DSHS/RDA) 
 nord...@dshs.wa.govmailto:nord...@dshs.wa.gov wrote:
  -Original Message-
  From: r-help-boun...@r-project.orgmailto:r-help-boun...@r-project.org
 [mailto:r-help-bounces@r-mailto:r-help-bounces@r-
  project.orghttp://project.org] On Behalf Of C W
  Sent: Monday, March 18, 2013 11:27 AM
  To: r-help
  Subject: [R] How to stop set.seed() besides exiting out of R?
 
  Hi list,
 
  I am curious how to stop the set.seed(), I don't want the same repeated
  random number.  I know I can set it to a different seed, but I don't
  want
  to go through the trouble of setting different seed every time.
 
  Thanks,
  Mike
 
 Can you show us how you are using set.seed() that results in getting the
 same sequence repeatedly?  If you are doing simulations in a loop, then set
 the seed once, outside the loop.  Otherwise, I am not sure what you are
 doing that causes problems.  A reproducible example would really help.

 Dan

 Daniel J. Nordlund
 Washington State Department of Social and Health Services
 Planning, Performance, and Accountability
 Research and Data Analysis Division
 Olympia, WA 98504-5204

 __
 R-help@r-project.orgmailto:R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


 [[alternative HTML version deleted]]


 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.



[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] OrgMassSpecR peak area issue

2013-03-18 Thread Claudia Beleites
Hi Bryan,

Division by 2 is correct, comes from trapezoid calculation. 
The modulo line is a funny way of producing c (2 : n, 1)

Best,

Claudia


Am Mon, 18 Mar 2013 15:00:06 -0400
schrieb Bryan Hanson han...@depauw.edu:

 If you type DrawChromatogram you can see the method used to calculate
 the peak area.  Looks to me like you could easily hack it if you
 wanted.  The relevant part about peak areas is this:
 
 for (j in 1:n) {
 k - (j%%n) + 1
 x[j] - peakTime[j] * peakIntensity[k] - peakTime[k] * 
 peakIntensity[j]
 }
 peakArea[i] - abs(sum(x)/2)
 
 which looks pretty standard to me, though I'm not clear right off the
 top of my head why they are dividing by 2.  You can always contact
 the maintainer.
 
 Bryan
 
 On Mar 18, 2013, at 1:34 PM, Christopher Beaver
 christopher.bea...@gmail.com wrote:
 
  Hello!
  
  I am having an issue with the OrgMassSpecR package.  I run my HPLC
  using a DAD detector.  My raw data is exported form chemstation as
  a csv file.  I then upload the csv into Rstudio no problem.  Using
  the DrawChromatogram function, I get a nice chromatogram, and my
  retention time, peak area, and apex intensity values are given as
  well.
  
  The problem comes with the peak area value given. The peak area is
  much smaller than a value that would make sense.  My peak area
  value is actually less than my apex intensity value.  Is this
  because I am using a DAD detector rather than an MS? If so, is
  there a simply way to edit the peak area equation so that it will
  also work with absorbance values?
  
  Any help is greatly appreciated.
  
  Thanks for your time.
  
  Chris Beaver
  
  [[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.



-- 
Claudia Beleites
Spectroscopy/Imaging
Institute of Photonic Technology 
Albert-Einstein-Str. 9
07745 Jena
Germany

email: claudia.belei...@ipht-jena.de
phone: +49 3641 206-133
fax:   +49 2641 206-399

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 with NA

2013-03-18 Thread Pietro

Yes, it's true Berend!

What i do is simply use read.xlsx  function

db - read.xlsx2(c:/mydb.xlsx,1,as.data.frame=T)

This is excel file i use:
http://dl.dropbox.com/u/102669/mydb.xlsx

I can't find  a way to import as numeric.
My objective is to be able to work (in R) with my NA's


At 18.46 18/03/2013, Berend Hasselman wrote:


On 18-03-2013, at 16:49, Pete freeri...@gmail.com wrote:


 I have this little data.frame

 http://dl.dropbox.com/u/102669/nanotna.rdata

 Two column contains NA, so the best thing to do is use na.locf 
function (with

 fromLast = T)

 But locf function doesn't work because NA in my data.frame are 
not recognized as

 real NA.

 Is there a way to substitute fake NA with real NA? In this case 
na.locf function

 should work


Your data are all characters. Do

str(db)

to see that. What is probably supposed to be numeric is also character,
Somehow you have managed to read in data that R thinks is all chr.
Your NA are NA in reality: a character string NA.

You will have to review the method you used to get the data into R.
And make sure that what you want to be numeric is indeed numeric.
Then you can start to think about doing something about the NA's.

Berend


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 stop set.seed() besides exiting out of R?

2013-03-18 Thread Nordlund, Dan (DSHS/RDA)
No, you cannot unset the seed.  You can set it to a different value, but a the 
random number generators always need a starting seed.  If you don’t set one, 
R will set one for you , you just won’t know what it is.  And as a practical 
matter, given a sequence of random numbers you can’t tell what the starting 
seed was.  That is the point of good random number generators.  Each sequence 
of random numbers for most intents and purposes can be considered independent 
from previous sets of numbers.

Hope this is helpful,

Dan

Daniel J. Nordlund
Washington State Department of Social and Health Services
Planning, Performance, and Accountability
Research and Data Analysis Division
Olympia, WA 98504-5204

From: C W [mailto:tmrs...@gmail.com]
Sent: Monday, March 18, 2013 12:19 PM
To: Nordlund, Dan (DSHS/RDA)
Cc: r-help
Subject: Re: [R] How to stop set.seed() besides exiting out of R?

Yes, I agree with you.  I guess what I was really looking for is a function 
like UNset.seed()?

By having set.seed(), I can have reproducible code.  But what if I want to 
check my work against what's produced from set.seed(100)?

I really want to escape from the shadow of set.seed(), can I unset it?
On Mon, Mar 18, 2013 at 3:07 PM, Nordlund, Dan (DSHS/RDA) 
nord...@dshs.wa.govmailto:nord...@dshs.wa.gov wrote:
As I understand it, how R “‘normally” does it is to use the system clock 
to set the seed once per session, unless you use set.seed() to set a new seed. 
You chose to set the seed to a different value.  But from that point on, the 
pseudo random number generation continues  in the same way it “normally” 
does.  In your code below, each of your 100 histograms will be different.  If 
you then execute the for loop again (but not the set.seed(100) statement), you 
will get a different set of histograms.  The only way you would be “confined 
to set.seed(100)” is if you keep resetting the seed to 100.

Dan

Daniel J. Nordlund
Washington State Department of Social and Health Services
Planning, Performance, and Accountability
Research and Data Analysis Division
Olympia, WA 98504-5204
From: C W [mailto:tmrs...@gmail.commailto:tmrs...@gmail.com]
Sent: Monday, March 18, 2013 11:50 AM
To: Nordlund, Dan (DSHS/RDA)
Cc: r-help
Subject: Re: [R] How to stop set.seed() besides exiting out of R?

set.seed(100)
for (i in 1:100){
a - rnorm(1000, mean=0, sd=1)
hist(a)
}

#Now say, I want to simulate without being confined to set.seed(100), I just 
want to get a simulation like how R normally does it.

Mike
On Mon, Mar 18, 2013 at 2:40 PM, Nordlund, Dan (DSHS/RDA) 
nord...@dshs.wa.govmailto:nord...@dshs.wa.govmailto:nord...@dshs.wa.govmailto:nord...@dshs.wa.gov
 wrote:
 -Original Message-
 From: 
 r-help-boun...@r-project.orgmailto:r-help-boun...@r-project.orgmailto:r-help-boun...@r-project.orgmailto:r-help-boun...@r-project.org
  
 [mailto:r-help-bounces@r-mailto:r-help-bounces@r-mailto:r-help-bounces@r-mailto:r-help-bounces@r-
 project.orghttp://project.orghttp://project.org] On Behalf Of C W
 Sent: Monday, March 18, 2013 11:27 AM
 To: r-help
 Subject: [R] How to stop set.seed() besides exiting out of R?

 Hi list,

 I am curious how to stop the set.seed(), I don't want the same repeated
 random number.  I know I can set it to a different seed, but I don't
 want
 to go through the trouble of setting different seed every time.

 Thanks,
 Mike

Can you show us how you are using set.seed() that results in getting the same 
sequence repeatedly?  If you are doing simulations in a loop, then set the seed 
once, outside the loop.  Otherwise, I am not sure what you are doing that 
causes problems.  A reproducible example would really help.

Dan

Daniel J. Nordlund
Washington State Department of Social and Health Services
Planning, Performance, and Accountability
Research and Data Analysis Division
Olympia, WA 98504-5204

__
R-help@r-project.orgmailto:R-help@r-project.orgmailto:R-help@r-project.orgmailto:R-help@r-project.org
 mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[[alternative HTML version deleted]]


__
R-help@r-project.orgmailto:R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Creating SubScale Scores

2013-03-18 Thread Amanda Charles
Hello all, 
I have just completed a principal components analysis and now which to create 
sub scales out of my factors. Does anyone know how to do this in R? 
Thanks!!
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Counting confidence intervals

2013-03-18 Thread Rui Barradas

Hello,

There _is_ a function ?within. Maybe your function can be named 'between'

Rui Barradas

Em 18-03-2013 16:16, S Ellison escreveu:

I want to cont how many
times a number say 12 lies in the interval. Can anyone assist?


Has anyone else ever wished there was a moderately general 'inside' or 'within' 
function in R for this problem?

For example, something that behaves more or less like

within - function(x, interval=NULL, closed=c(TRUE, TRUE), lower=min(interval), 
upper=max(interval)) {
#interval must be a length 2 vector
#closed is taken in the order (lower, upper)
#lower and upper may be vectors and will be recycled (by  etc) if 
not of length length(x)

low.comp - if(closed[1]) = else 
high.comp - if(closed[2]) = else 

do.call(low.comp, list(lower, x))  do.call(high.comp, list(upper, x))
}


#Examples
within(1:5, c(2,4))

within(1:5, c(2,4), closed=c(FALSE, TRUE))

within(1:5, lower=5:1, upper=10:14)


S Ellison
LGC

***
This email and any attachments are confidential. Any use...{{dropped:8}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 stop set.seed() besides exiting out of R?

2013-03-18 Thread William Dunlap
I am not sure of this, but I think you can unset the seed by removing the 
dataset .Random.seed
from the global environment.  E.g.,

 set.seed(1)
 runif(5)
[1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819
 rm(list=.Random.seed, envir=globalenv())
 runif(5)
[1] 0.5952379 0.3355091 0.8820192 0.7633754 0.8064312


 set.seed(1)
 runif(5) # same as before
[1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819
 rm(list=.Random.seed, envir=globalenv())
 runif(5) # different
[1] 0.52023610 0.73407695 0.08824484 0.26977430 0.80089250

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 Nordlund, Dan (DSHS/RDA)
 Sent: Monday, March 18, 2013 12:44 PM
 To: C W
 Cc: r-help
 Subject: Re: [R] How to stop set.seed() besides exiting out of R?
 
 No, you cannot unset the seed.  You can set it to a different value, but a 
 the random
 number generators always need a starting seed.  If you don’t set one, R will 
 set one for
 you , you just won’t know what it is.  And as a practical matter, given a 
 sequence of
 random numbers you can’t tell what the starting seed was.  That is the point 
 of good
 random number generators.  Each sequence of random numbers for most intents 
 and
 purposes can be considered independent from previous sets of numbers.
 
 Hope this is helpful,
 
 Dan
 
 Daniel J. Nordlund
 Washington State Department of Social and Health Services
 Planning, Performance, and Accountability
 Research and Data Analysis Division
 Olympia, WA 98504-5204
 
 From: C W [mailto:tmrs...@gmail.com]
 Sent: Monday, March 18, 2013 12:19 PM
 To: Nordlund, Dan (DSHS/RDA)
 Cc: r-help
 Subject: Re: [R] How to stop set.seed() besides exiting out of R?
 
 Yes, I agree with you.  I guess what I was really looking for is a function 
 like
 UNset.seed()?
 
 By having set.seed(), I can have reproducible code.  But what if I want to 
 check my work
 against what's produced from set.seed(100)?
 
 I really want to escape from the shadow of set.seed(), can I unset it?
 On Mon, Mar 18, 2013 at 3:07 PM, Nordlund, Dan (DSHS/RDA)
 nord...@dshs.wa.govmailto:nord...@dshs.wa.gov wrote:
 As I understand it, how R “‘normally” does it is to use the system clock to 
 set the seed
 once per session, unless you use set.seed() to set a new seed. You chose to 
 set the seed
 to a different value.  But from that point on, the pseudo random number 
 generation
 continues  in the same way it “normally” does.  In your code below, each of 
 your 100
 histograms will be different.  If you then execute the for loop again (but 
 not the
 set.seed(100) statement), you will get a different set of histograms.  The 
 only way you
 would be “confined to set.seed(100)” is if you keep resetting the seed to 100.
 
 Dan
 
 Daniel J. Nordlund
 Washington State Department of Social and Health Services
 Planning, Performance, and Accountability
 Research and Data Analysis Division
 Olympia, WA 98504-5204
 From: C W [mailto:tmrs...@gmail.commailto:tmrs...@gmail.com]
 Sent: Monday, March 18, 2013 11:50 AM
 To: Nordlund, Dan (DSHS/RDA)
 Cc: r-help
 Subject: Re: [R] How to stop set.seed() besides exiting out of R?
 
 set.seed(100)
 for (i in 1:100){
 a - rnorm(1000, mean=0, sd=1)
 hist(a)
 }
 
 #Now say, I want to simulate without being confined to set.seed(100), I just 
 want to get
 a simulation like how R normally does it.
 
 Mike
 On Mon, Mar 18, 2013 at 2:40 PM, Nordlund, Dan (DSHS/RDA)
 nord...@dshs.wa.govmailto:nord...@dshs.wa.govmailto:nord...@dshs.wa.gov
 mailto:nord...@dshs.wa.gov wrote:
  -Original Message-
  From: 
  r-help-boun...@r-project.orgmailto:r-help-boun...@r-project.orgmailto:r-
 help-boun...@r-project.orgmailto:r-help-boun...@r-project.org 
 [mailto:r-help-
 bounces@r-mailto:r-help-bounces@r-mailto:r-help-bounces@r-mailto:r-help-
 bounces@r-
  project.orghttp://project.orghttp://project.org] On Behalf Of C W
  Sent: Monday, March 18, 2013 11:27 AM
  To: r-help
  Subject: [R] How to stop set.seed() besides exiting out of R?
 
  Hi list,
 
  I am curious how to stop the set.seed(), I don't want the same repeated
  random number.  I know I can set it to a different seed, but I don't
  want
  to go through the trouble of setting different seed every time.
 
  Thanks,
  Mike
 
 Can you show us how you are using set.seed() that results in getting the same 
 sequence
 repeatedly?  If you are doing simulations in a loop, then set the seed once, 
 outside the
 loop.  Otherwise, I am not sure what you are doing that causes problems.  A 
 reproducible
 example would really help.
 
 Dan
 
 Daniel J. Nordlund
 Washington State Department of Social and Health Services
 Planning, Performance, and Accountability
 Research and Data Analysis Division
 Olympia, WA 98504-5204
 
 __
 R-help@r-project.orgmailto:R-help@r-project.orgmailto:R-help@r-
 project.orgmailto:R-help@r-project.org 

Re: [R] How to stop set.seed() besides exiting out of R?

2013-03-18 Thread Jeff Newmiller
Your request is meaningless.  The seed itself is effectively overwritten each 
time a random number is requested. It is only the repeatability of the sequence 
of random numbers following the set.seed call that is reproducible. You can set 
the seed to something else, but there is always a seed and it changes as 
numbers are requested.
---
Jeff NewmillerThe .   .  Go Live...
DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live Go...
  Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
--- 
Sent from my phone. Please excuse my brevity.

C W tmrs...@gmail.com wrote:

Yes, I agree with you.  I guess what I was really looking for is a
function
like UNset.seed()?

By having set.seed(), I can have reproducible code.  But what if I want
to
check my work against what's produced from set.seed(100)?

I really want to escape from the shadow of set.seed(), can I unset it?

On Mon, Mar 18, 2013 at 3:07 PM, Nordlund, Dan (DSHS/RDA) 
nord...@dshs.wa.gov wrote:

 As I understand it, how R ��normally� does it is to use the system
clock
 to set the seed once per session, unless you use set.seed() to set a
new
 seed. You chose to set the seed to a different value.  But from that
point
 on, the pseudo random number generation continues  in the same way it
 �normally� does.  In your code below, each of your 100 histograms
will be
 different.  If you then execute the for loop again (but not the
 set.seed(100) statement), you will get a different set of histograms.
 The
 only way you would be �confined to set.seed(100)� is if you keep
resetting
 the seed to 100.

 Dan

 Daniel J. Nordlund
 Washington State Department of Social and Health Services
 Planning, Performance, and Accountability
 Research and Data Analysis Division
 Olympia, WA 98504-5204

 From: C W [mailto:tmrs...@gmail.com]
 Sent: Monday, March 18, 2013 11:50 AM
 To: Nordlund, Dan (DSHS/RDA)
 Cc: r-help
 Subject: Re: [R] How to stop set.seed() besides exiting out of R?

 set.seed(100)
 for (i in 1:100){
 a - rnorm(1000, mean=0, sd=1)
 hist(a)
 }

 #Now say, I want to simulate without being confined to set.seed(100),
I
 just want to get a simulation like how R normally does it.

 Mike

 On Mon, Mar 18, 2013 at 2:40 PM, Nordlund, Dan (DSHS/RDA) 
 nord...@dshs.wa.govmailto:nord...@dshs.wa.gov wrote:
  -Original Message-
  From:
r-help-boun...@r-project.orgmailto:r-help-boun...@r-project.org
 [mailto:r-help-bounces@r-mailto:r-help-bounces@r-
  project.orghttp://project.org] On Behalf Of C W
  Sent: Monday, March 18, 2013 11:27 AM
  To: r-help
  Subject: [R] How to stop set.seed() besides exiting out of R?
 
  Hi list,
 
  I am curious how to stop the set.seed(), I don't want the same
repeated
  random number.  I know I can set it to a different seed, but I
don't
  want
  to go through the trouble of setting different seed every time.
 
  Thanks,
  Mike
 
 Can you show us how you are using set.seed() that results in getting
the
 same sequence repeatedly?  If you are doing simulations in a loop,
then set
 the seed once, outside the loop.  Otherwise, I am not sure what you
are
 doing that causes problems.  A reproducible example would really
help.

 Dan

 Daniel J. Nordlund
 Washington State Department of Social and Health Services
 Planning, Performance, and Accountability
 Research and Data Analysis Division
 Olympia, WA 98504-5204

 __
 R-help@r-project.orgmailto:R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


 [[alternative HTML version deleted]]


 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.



   [[alternative HTML version deleted]]





__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible 

Re: [R] How to stop set.seed() besides exiting out of R?

2013-03-18 Thread C W
Thanks, William.  You read my mind.
Cheers,
Mike

On Mon, Mar 18, 2013 at 4:00 PM, Jeff Newmiller jdnew...@dcn.davis.ca.uswrote:

 Your request is meaningless.  The seed itself is effectively overwritten
 each time a random number is requested. It is only the repeatability of the
 sequence of random numbers following the set.seed call that is
 reproducible. You can set the seed to something else, but there is always a
 seed and it changes as numbers are requested.
 ---
 Jeff NewmillerThe .   .  Go Live...
 DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live
 Go...
   Live:   OO#.. Dead: OO#..  Playing
 Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
 /Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
 ---
 Sent from my phone. Please excuse my brevity.

 C W tmrs...@gmail.com wrote:

 Yes, I agree with you.  I guess what I was really looking for is a
 function
 like UNset.seed()?
 
 By having set.seed(), I can have reproducible code.  But what if I want
 to
 check my work against what's produced from set.seed(100)?
 
 I really want to escape from the shadow of set.seed(), can I unset it?
 
 On Mon, Mar 18, 2013 at 3:07 PM, Nordlund, Dan (DSHS/RDA) 
 nord...@dshs.wa.gov wrote:
 
  As I understand it, how R ��normally� does it is to use the system
 clock
  to set the seed once per session, unless you use set.seed() to set a
 new
  seed. You chose to set the seed to a different value.  But from that
 point
  on, the pseudo random number generation continues  in the same way it
  �normally� does.  In your code below, each of your 100 histograms
 will be
  different.  If you then execute the for loop again (but not the
  set.seed(100) statement), you will get a different set of histograms.
  The
  only way you would be �confined to set.seed(100)� is if you keep
 resetting
  the seed to 100.
 
  Dan
 
  Daniel J. Nordlund
  Washington State Department of Social and Health Services
  Planning, Performance, and Accountability
  Research and Data Analysis Division
  Olympia, WA 98504-5204
 
  From: C W [mailto:tmrs...@gmail.com]
  Sent: Monday, March 18, 2013 11:50 AM
  To: Nordlund, Dan (DSHS/RDA)
  Cc: r-help
  Subject: Re: [R] How to stop set.seed() besides exiting out of R?
 
  set.seed(100)
  for (i in 1:100){
  a - rnorm(1000, mean=0, sd=1)
  hist(a)
  }
 
  #Now say, I want to simulate without being confined to set.seed(100),
 I
  just want to get a simulation like how R normally does it.
 
  Mike
 
  On Mon, Mar 18, 2013 at 2:40 PM, Nordlund, Dan (DSHS/RDA) 
  nord...@dshs.wa.govmailto:nord...@dshs.wa.gov wrote:
   -Original Message-
   From:
 r-help-boun...@r-project.orgmailto:r-help-boun...@r-project.org
  [mailto:r-help-bounces@r-mailto:r-help-bounces@r-
   project.orghttp://project.org] On Behalf Of C W
   Sent: Monday, March 18, 2013 11:27 AM
   To: r-help
   Subject: [R] How to stop set.seed() besides exiting out of R?
  
   Hi list,
  
   I am curious how to stop the set.seed(), I don't want the same
 repeated
   random number.  I know I can set it to a different seed, but I
 don't
   want
   to go through the trouble of setting different seed every time.
  
   Thanks,
   Mike
  
  Can you show us how you are using set.seed() that results in getting
 the
  same sequence repeatedly?  If you are doing simulations in a loop,
 then set
  the seed once, outside the loop.  Otherwise, I am not sure what you
 are
  doing that causes problems.  A reproducible example would really
 help.
 
  Dan
 
  Daniel J. Nordlund
  Washington State Department of Social and Health Services
  Planning, Performance, and Accountability
  Research and Data Analysis Division
  Olympia, WA 98504-5204
 
  __
  R-help@r-project.orgmailto:R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide
  http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
 
 
  [[alternative HTML version deleted]]
 
 
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide
  http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
 
 
 
[[alternative HTML version deleted]]
 
 
 
 
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, 

Re: [R] data.frame with NA

2013-03-18 Thread David L Carlson
Try this

Open the spreadsheet in Excel. Select all of the data click Copy. Don't
close Excel.

Open R and type the following command:

 Foglio1 - read.table(clipboard-128, header=TRUE, sep=\t)

Now take a look at the structure of the data.frame

 str(Foglio1)
'data.frame':   1489 obs. of  15 variables:
 $ Date: Factor w/ 1489 levels 1/10/2002,1/10/2003,..: 1275 1291 1295
1299 1304 1309 1321 1325 1329 1337 ...
 $ a   : num  202 201 202 201 202 ...
 $ b   : num  231 230 230 230 232 ...
 $ c   : num  177 179 181 180 182 ...
 $ d   : num  277 277 276 276 275 ...
 $ e   : num  NA NA NA NA NA NA NA NA NA NA ...
 $ f   : num  275 277 279 279 279 ...
 $ g   : num  91.7 90.7 90.8 91.1 91 ...
 $ h   : num  11446 11258 11280 11396 11127 ...
 $ i   : num  388 389 393 392 393 ...
 $ l   : num  93.2 94 92.4 93.4 93.1 ...
 $ m   : num  128 127 126 129 130 ...
 $ n   : num  NA NA NA NA NA NA NA NA NA NA ...
 $ o   : num  133 133 133 133 133 ...
 $ p   : num  107 107 107 107 107 ...

--
David L Carlson
Associate Professor of Anthropology
Texas AM University
College Station, TX 77843-4352

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
 project.org] On Behalf Of Pietro
 Sent: Monday, March 18, 2013 1:57 PM
 To: Berend Hasselman
 Cc: r-h...@stat.math.ethz.ch
 Subject: Re: [R] data.frame with NA
 
 Yes, it's true Berend!
 
 What i do is simply use read.xlsx  function
 
 db - read.xlsx2(c:/mydb.xlsx,1,as.data.frame=T)
 
 This is excel file i use:
 http://dl.dropbox.com/u/102669/mydb.xlsx
 
 I can't find  a way to import as numeric.
 My objective is to be able to work (in R) with my NA's
 
 
 At 18.46 18/03/2013, Berend Hasselman wrote:
 
 On 18-03-2013, at 16:49, Pete freeri...@gmail.com wrote:
 
  
   I have this little data.frame
  
   http://dl.dropbox.com/u/102669/nanotna.rdata
  
   Two column contains NA, so the best thing to do is use na.locf
  function (with
   fromLast = T)
  
   But locf function doesn't work because NA in my data.frame are
  not recognized as
   real NA.
  
   Is there a way to substitute fake NA with real NA? In this case
  na.locf function
   should work
  
 
 Your data are all characters. Do
 
 str(db)
 
 to see that. What is probably supposed to be numeric is also
 character,
 Somehow you have managed to read in data that R thinks is all chr.
 Your NA are NA in reality: a character string NA.
 
 You will have to review the method you used to get the data into R.
 And make sure that what you want to be numeric is indeed numeric.
 Then you can start to think about doing something about the NA's.
 
 Berend
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] How to stop set.seed() besides exiting out of R?

2013-03-18 Thread Prof Brian Ripley

On 18/03/2013 19:59, William Dunlap wrote:

I am not sure of this, but I think you can unset the seed by removing the 
dataset .Random.seed
from the global environment.  E.g.,


set.seed(1)
runif(5)

[1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819

rm(list=.Random.seed, envir=globalenv())
runif(5)

[1] 0.5952379 0.3355091 0.8820192 0.7633754 0.8064312



set.seed(1)
runif(5) # same as before

[1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819

rm(list=.Random.seed, envir=globalenv())
runif(5) # different

[1] 0.52023610 0.73407695 0.08824484 0.26977430 0.80089250


Yes, almost all the time.  R does keep a copy in memory when it is using 
it and writes it back out when done with it: so assuming nothing 
asynchronous is going on (e.g. a GUI callback) removing .Random.seed 
will cause the RNG to be re-initialized at next use.




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 Nordlund, Dan (DSHS/RDA)
Sent: Monday, March 18, 2013 12:44 PM
To: C W
Cc: r-help
Subject: Re: [R] How to stop set.seed() besides exiting out of R?

No, you cannot unset the seed.  You can set it to a different value, but a the 
random
number generators always need a starting seed.  If you don’t set one, R will 
set one for
you , you just won’t know what it is.  And as a practical matter, given a 
sequence of
random numbers you can’t tell what the starting seed was.  That is the point of 
good
random number generators.  Each sequence of random numbers for most intents and
purposes can be considered independent from previous sets of numbers.

Hope this is helpful,

Dan

Daniel J. Nordlund
Washington State Department of Social and Health Services
Planning, Performance, and Accountability
Research and Data Analysis Division
Olympia, WA 98504-5204

From: C W [mailto:tmrs...@gmail.com]
Sent: Monday, March 18, 2013 12:19 PM
To: Nordlund, Dan (DSHS/RDA)
Cc: r-help
Subject: Re: [R] How to stop set.seed() besides exiting out of R?

Yes, I agree with you.  I guess what I was really looking for is a function like
UNset.seed()?

By having set.seed(), I can have reproducible code.  But what if I want to 
check my work
against what's produced from set.seed(100)?

I really want to escape from the shadow of set.seed(), can I unset it?
On Mon, Mar 18, 2013 at 3:07 PM, Nordlund, Dan (DSHS/RDA)
nord...@dshs.wa.govmailto:nord...@dshs.wa.gov wrote:
As I understand it, how R “‘normally” does it is to use the system clock to set 
the seed
once per session, unless you use set.seed() to set a new seed. You chose to set 
the seed
to a different value.  But from that point on, the pseudo random number 
generation
continues  in the same way it “normally” does.  In your code below, each of 
your 100
histograms will be different.  If you then execute the for loop again (but not 
the
set.seed(100) statement), you will get a different set of histograms.  The only 
way you
would be “confined to set.seed(100)” is if you keep resetting the seed to 100.

Dan

Daniel J. Nordlund
Washington State Department of Social and Health Services
Planning, Performance, and Accountability
Research and Data Analysis Division
Olympia, WA 98504-5204
From: C W [mailto:tmrs...@gmail.commailto:tmrs...@gmail.com]
Sent: Monday, March 18, 2013 11:50 AM
To: Nordlund, Dan (DSHS/RDA)
Cc: r-help
Subject: Re: [R] How to stop set.seed() besides exiting out of R?

set.seed(100)
for (i in 1:100){
 a - rnorm(1000, mean=0, sd=1)
 hist(a)
}

#Now say, I want to simulate without being confined to set.seed(100), I just 
want to get
a simulation like how R normally does it.

Mike
On Mon, Mar 18, 2013 at 2:40 PM, Nordlund, Dan (DSHS/RDA)
nord...@dshs.wa.govmailto:nord...@dshs.wa.govmailto:nord...@dshs.wa.gov
mailto:nord...@dshs.wa.gov wrote:

-Original Message-
From: 
r-help-boun...@r-project.orgmailto:r-help-boun...@r-project.orgmailto:r-

help-boun...@r-project.orgmailto:r-help-boun...@r-project.org [mailto:r-help-
bounces@r-mailto:r-help-bounces@r-mailto:r-help-bounces@r-mailto:r-help-
bounces@r-

project.orghttp://project.orghttp://project.org] On Behalf Of C W
Sent: Monday, March 18, 2013 11:27 AM
To: r-help
Subject: [R] How to stop set.seed() besides exiting out of R?

Hi list,

I am curious how to stop the set.seed(), I don't want the same repeated
random number.  I know I can set it to a different seed, but I don't
want
to go through the trouble of setting different seed every time.

Thanks,
Mike


Can you show us how you are using set.seed() that results in getting the same 
sequence
repeatedly?  If you are doing simulations in a loop, then set the seed once, 
outside the
loop.  Otherwise, I am not sure what you are doing that causes problems.  A 
reproducible
example would really help.

Dan

Daniel J. Nordlund
Washington State Department of Social and Health Services
Planning, Performance, and Accountability
Research and Data Analysis Division

Re: [R] how to plot u-v wind by R?

2013-03-18 Thread Greg Snow
The my.symbols and ms.arrows functions in the TeachingDemos package may
help.


On Mon, Mar 18, 2013 at 3:50 AM, Jie Tang totang...@gmail.com wrote:

 hi R users:
I have a dataset including u wind in x-axis and v wind in y-axis.
  How can I plot the u,v wind data in vector or  barb figure?
  which command ?
 thank you .

 --
 TANG Jie

 [[alternative HTML version deleted]]

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.




-- 
Gregory (Greg) L. Snow Ph.D.
538...@gmail.com

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Superscript followed by number then superscript in text

2013-03-18 Thread Benjamin Gillespie
Hi all,

I'm having problems finding the correct format for a command.

I would like to write some text on a plot.

I'm using the following command:

text(x,y,text here, srt=90)

I would like the text to read:

capacity 10^3 m^3

(with ^ denoting superscript (i.e. each '3' as superscript).

I've tried fiddling around with expression(paste(etc to no avail. I receive 
errors such as: Error: unexpected numeric constant...

Anyone had experience with this before? Any suggestions would be great.

Many thanks in advance,

Ben Gillespie
Research Postgraduate
 
School of Geography
University of Leeds
Leeds
LS2 9JT
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 with NA

2013-03-18 Thread David L Carlson
It appears that you MUST use the colClasses= argument with read.xlsx2:

Foglio1 - read.xlsx2(mydb.xlsx, 1, colClasses=c(Date, rep(numeric,
14)))

However, e and n are converted to NaN not NA so you would need to convert
those columns (at least, I didn't check for missing values in the other
columns):

 Foglio1$e - ifelse(is.nan(Foglio1$e), NA, Foglio1$e)
 Foglio1$n - ifelse(is.nan(Foglio1$n), NA, Foglio1$n)
 str(Foglio1)
'data.frame':   1489 obs. of  15 variables:
 $ Date: Date, format: 2001-08-17 2001-08-20 ...
 $ a   : num  202 201 202 201 202 ...
 $ b   : num  231 230 230 230 232 ...
 $ c   : num  177 179 181 180 182 ...
 $ d   : num  277 277 276 276 275 ...
 $ e   : num  NA NA NA NA NA NA NA NA NA NA ...
 $ f   : num  275 277 279 279 279 ...
 $ g   : num  91.7 90.7 90.8 91.1 91 ...
 $ h   : num  11446 11258 11280 11396 11127 ...
 $ i   : num  388 389 393 392 393 ...
 $ l   : num  93.2 94 92.4 93.4 93.1 ...
 $ m   : num  128 127 126 129 130 ...
 $ n   : num  NA NA NA NA NA NA NA NA NA NA ...
 $ o   : num  133 133 133 133 133 ...
 $ p   : num  107 107 107 107 107 ...

---
David


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
 project.org] On Behalf Of David L Carlson
 Sent: Monday, March 18, 2013 3:22 PM
 To: 'Pietro'; 'Berend Hasselman'
 Cc: r-h...@stat.math.ethz.ch
 Subject: Re: [R] data.frame with NA
 
 Try this
 
 Open the spreadsheet in Excel. Select all of the data click Copy. Don't
 close Excel.
 
 Open R and type the following command:
 
  Foglio1 - read.table(clipboard-128, header=TRUE, sep=\t)
 
 Now take a look at the structure of the data.frame
 
  str(Foglio1)
 'data.frame':   1489 obs. of  15 variables:
  $ Date: Factor w/ 1489 levels 1/10/2002,1/10/2003,..: 1275 1291
 1295
 1299 1304 1309 1321 1325 1329 1337 ...
  $ a   : num  202 201 202 201 202 ...
  $ b   : num  231 230 230 230 232 ...
  $ c   : num  177 179 181 180 182 ...
  $ d   : num  277 277 276 276 275 ...
  $ e   : num  NA NA NA NA NA NA NA NA NA NA ...
  $ f   : num  275 277 279 279 279 ...
  $ g   : num  91.7 90.7 90.8 91.1 91 ...
  $ h   : num  11446 11258 11280 11396 11127 ...
  $ i   : num  388 389 393 392 393 ...
  $ l   : num  93.2 94 92.4 93.4 93.1 ...
  $ m   : num  128 127 126 129 130 ...
  $ n   : num  NA NA NA NA NA NA NA NA NA NA ...
  $ o   : num  133 133 133 133 133 ...
  $ p   : num  107 107 107 107 107 ...
 
 --
 David L Carlson
 Associate Professor of Anthropology
 Texas AM University
 College Station, TX 77843-4352
 
  -Original Message-
  From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
  project.org] On Behalf Of Pietro
  Sent: Monday, March 18, 2013 1:57 PM
  To: Berend Hasselman
  Cc: r-h...@stat.math.ethz.ch
  Subject: Re: [R] data.frame with NA
 
  Yes, it's true Berend!
 
  What i do is simply use read.xlsx  function
 
  db - read.xlsx2(c:/mydb.xlsx,1,as.data.frame=T)
 
  This is excel file i use:
  http://dl.dropbox.com/u/102669/mydb.xlsx
 
  I can't find  a way to import as numeric.
  My objective is to be able to work (in R) with my NA's
 
 
  At 18.46 18/03/2013, Berend Hasselman wrote:
 
  On 18-03-2013, at 16:49, Pete freeri...@gmail.com wrote:
  
   
I have this little data.frame
   
http://dl.dropbox.com/u/102669/nanotna.rdata
   
Two column contains NA, so the best thing to do is use na.locf
   function (with
fromLast = T)
   
But locf function doesn't work because NA in my data.frame are
   not recognized as
real NA.
   
Is there a way to substitute fake NA with real NA? In this case
   na.locf function
should work
   
  
  Your data are all characters. Do
  
  str(db)
  
  to see that. What is probably supposed to be numeric is also
  character,
  Somehow you have managed to read in data that R thinks is all chr.
  Your NA are NA in reality: a character string NA.
  
  You will have to review the method you used to get the data into R.
  And make sure that what you want to be numeric is indeed numeric.
  Then you can start to think about doing something about the NA's.
  
  Berend
 
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide http://www.R-project.org/posting-
  guide.html
  and provide commented, minimal, self-contained, reproducible code.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-
 guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 stop set.seed() besides exiting out of R?

2013-03-18 Thread William Dunlap
 -Original Message-
 From: Prof Brian Ripley [mailto:rip...@stats.ox.ac.uk]
 Sent: Monday, March 18, 2013 1:22 PM
 To: William Dunlap
 Cc: r-help
 Subject: Re: [R] How to stop set.seed() besides exiting out of R?
 
 On 18/03/2013 19:59, William Dunlap wrote:
  I am not sure of this, but I think you can unset the seed by removing the 
  dataset
 .Random.seed
  from the global environment.  E.g.,
 
  set.seed(1)
  runif(5)
  [1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819
  rm(list=.Random.seed, envir=globalenv())
  runif(5)
  [1] 0.5952379 0.3355091 0.8820192 0.7633754 0.8064312
 
 
  set.seed(1)
  runif(5) # same as before
  [1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819
  rm(list=.Random.seed, envir=globalenv())
  runif(5) # different
  [1] 0.52023610 0.73407695 0.08824484 0.26977430 0.80089250
 
 Yes, almost all the time.  R does keep a copy in memory when it is using
 it and writes it back out when done with it: so assuming nothing
 asynchronous is going on (e.g. a GUI callback) removing .Random.seed
 will cause the RNG to be re-initialized at next use.

Thanks Brian,

S+'s set.seed() generates a new seed if it is called with no arguments so
as not to tempt people to directly fiddle around with magic datasets
like .Random.seed.

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com

 
 
  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 Nordlund, Dan (DSHS/RDA)
  Sent: Monday, March 18, 2013 12:44 PM
  To: C W
  Cc: r-help
  Subject: Re: [R] How to stop set.seed() besides exiting out of R?
 
  No, you cannot unset the seed.  You can set it to a different value, but a 
  the random
  number generators always need a starting seed.  If you don’t set one, R 
  will set one
 for
  you , you just won’t know what it is.  And as a practical matter, given a 
  sequence of
  random numbers you can’t tell what the starting seed was.  That is the 
  point of good
  random number generators.  Each sequence of random numbers for most 
  intents and
  purposes can be considered independent from previous sets of numbers.
 
  Hope this is helpful,
 
  Dan
 
  Daniel J. Nordlund
  Washington State Department of Social and Health Services
  Planning, Performance, and Accountability
  Research and Data Analysis Division
  Olympia, WA 98504-5204
 
  From: C W [mailto:tmrs...@gmail.com]
  Sent: Monday, March 18, 2013 12:19 PM
  To: Nordlund, Dan (DSHS/RDA)
  Cc: r-help
  Subject: Re: [R] How to stop set.seed() besides exiting out of R?
 
  Yes, I agree with you.  I guess what I was really looking for is a 
  function like
  UNset.seed()?
 
  By having set.seed(), I can have reproducible code.  But what if I want to 
  check my
 work
  against what's produced from set.seed(100)?
 
  I really want to escape from the shadow of set.seed(), can I unset it?
  On Mon, Mar 18, 2013 at 3:07 PM, Nordlund, Dan (DSHS/RDA)
  nord...@dshs.wa.govmailto:nord...@dshs.wa.gov wrote:
  As I understand it, how R “‘normally” does it is to use the system clock 
  to set the seed
  once per session, unless you use set.seed() to set a new seed. You chose 
  to set the
 seed
  to a different value.  But from that point on, the pseudo random number 
  generation
  continues  in the same way it “normally” does.  In your code below, each 
  of your 100
  histograms will be different.  If you then execute the for loop again (but 
  not the
  set.seed(100) statement), you will get a different set of histograms.  The 
  only way you
  would be “confined to set.seed(100)” is if you keep resetting the seed to 
  100.
 
  Dan
 
  Daniel J. Nordlund
  Washington State Department of Social and Health Services
  Planning, Performance, and Accountability
  Research and Data Analysis Division
  Olympia, WA 98504-5204
  From: C W [mailto:tmrs...@gmail.commailto:tmrs...@gmail.com]
  Sent: Monday, March 18, 2013 11:50 AM
  To: Nordlund, Dan (DSHS/RDA)
  Cc: r-help
  Subject: Re: [R] How to stop set.seed() besides exiting out of R?
 
  set.seed(100)
  for (i in 1:100){
   a - rnorm(1000, mean=0, sd=1)
   hist(a)
  }
 
  #Now say, I want to simulate without being confined to set.seed(100), I 
  just want to
 get
  a simulation like how R normally does it.
 
  Mike
  On Mon, Mar 18, 2013 at 2:40 PM, Nordlund, Dan (DSHS/RDA)
 
 nord...@dshs.wa.govmailto:nord...@dshs.wa.govmailto:nord...@dshs.wa.gov
  mailto:nord...@dshs.wa.gov wrote:
  -Original Message-
  From: r-help-boun...@r-project.orgmailto:r-help-bounces@r-
 project.orgmailto:r-
  help-boun...@r-project.orgmailto:r-help-boun...@r-project.org 
  [mailto:r-help-
  bounces@r-mailto:r-help-bounces@r-mailto:r-help-bounces@r-mailto:r-help-
  bounces@r-
  project.orghttp://project.orghttp://project.org] On Behalf Of C W
  Sent: Monday, March 18, 2013 11:27 AM
  To: r-help
  Subject: [R] How to stop set.seed() besides exiting out of R?
 
  Hi list,
 
  I 

Re: [R] Superscript followed by number then superscript in text

2013-03-18 Thread Thomas Lumley
On Tue, Mar 19, 2013 at 9:29 AM, Benjamin Gillespie gy...@leeds.ac.ukwrote:

 Hi all,

 I'm having problems finding the correct format for a command.

 I would like to write some text on a plot.

 I'm using the following command:

 text(x,y,text here, srt=90)

 I would like the text to read:

 capacity 10^3 m^3

 (with ^ denoting superscript (i.e. each '3' as superscript).


What did you try?

Anyhow, this works

text(1,1,expression(capacity~10^3~m^3))

  -thomas

-- 
Thomas Lumley
Professor of Biostatistics
University of Auckland

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Superscript followed by number then superscript in text

2013-03-18 Thread Rui Barradas

Hello,

Something like this?

plot(1, type = n)
text(1,1, expression(capacity ~ 10^3 ~ m^3))


Hope this helps,

Rui Barradas

Em 18-03-2013 20:29, Benjamin Gillespie escreveu:

Hi all,

I'm having problems finding the correct format for a command.

I would like to write some text on a plot.

I'm using the following command:

text(x,y,text here, srt=90)

I would like the text to read:

capacity 10^3 m^3

(with ^ denoting superscript (i.e. each '3' as superscript).

I've tried fiddling around with expression(paste(etc to no avail. I receive errors 
such as: Error: unexpected numeric constant...

Anyone had experience with this before? Any suggestions would be great.

Many thanks in advance,

Ben Gillespie
Research Postgraduate

School of Geography
University of Leeds
Leeds
LS2 9JT
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] melt with complications

2013-03-18 Thread arun
HI,
Not sure whether this is what you wanted.
library(reshape2)
result.dcast-dcast(meltTest,A+D~H,value.var=M)
attr(result.reshape,reshapeWide)-NULL
 row.names(result.reshape)-1:nrow(result.reshape)
 identical(result.reshape,result.dcast)
#[1] TRUE
result.dcast
#  A D  I  J  K  L
#1 B E  1  2  3  4
#2 B F  5  6  7  8
#3 B G  9 10 11 12
#4 C E 13 14 15 16
#5 C F 17 18 19 20
#6 C G 21 22 23 24

A.K.





- Original Message -
From: Richard M. Heiberger r...@temple.edu
To: r-help r-help@r-project.org; Hadley Wickham h.wick...@gmail.com
Cc: 
Sent: Monday, March 18, 2013 1:15 PM
Subject: [R] melt with complications

## Can someone suggest a simpler expression than either of these, with the
goal
## of taking a long matrix into a wide one with exactly one of the factors
converted to
## columns and all the rest retained as factors.  I want something that
generalizes beyond
## the three factors illustrated here.

## Rich

meltTest - data.frame(A=rep(c(B,C), each=12),
                       D=rep(c(E,F,G), each=4, times=2),
                       H=rep(c(I,J,K,L), times=6),
                       M=1:24)

meltTest

result.melt - do.call(rbind, {
  tmp - cast(D ~ H | A, value=M, data=meltTest)
  lapply(names(tmp), function(x) cbind(A=x, tmp[[x]])) ## explicit use of
name A
})
result.melt

result.reshape - reshape(meltTest, direction=wide, timevar=H,
idvar=c(A,D))
names(result.reshape)[3:6] - unique(as.character(meltTest$H))  ## explicit
use of name H
result.reshape

    [[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] melt with complications

2013-03-18 Thread Charles Berry
Richard M. Heiberger rmh at temple.edu writes:

 
 ## Can someone suggest a simpler expression than either of these, with the
 goal
 ## of taking a long matrix into a wide one with exactly one of the factors
 converted to
 ## columns and all the rest retained as factors.  I want something that
 generalizes beyond
 ## the three factors illustrated here.
 
 ## Rich
 
 meltTest - data.frame(A=rep(c(B,C), each=12),
D=rep(c(E,F,G), each=4, times=2),
H=rep(c(I,J,K,L), times=6),
M=1:24)


amat - ftable( xtabs( M ~ A + D + H, meltTest ),row.vars=1:2 )

amat is such a matrix with a few attributes added.

HTH,

Chuck

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Fit a mixture of lognormal and normal distributions

2013-03-18 Thread To . .

Hello

I am trying to find an automated way of fitting a mixture of normal and 
log-normal distributions to data which is clearly bimodal.
Here's a simulated example:
x.1-rnorm(6000, 2.4, 0.6)
x.2-rlnorm(1, 1.3,0.1)
X-c(x.1, x.2)
hist(X,100,freq=FALSE, ylim=c(0,1.5))lines(density(x.1), lty=2, lwd=2)
lines(density(x.2), lty=2, lwd=2)
lines(density(X), lty=4)

Currently i am using mixtools and just run:
library(mixtools)
mixmdl = normalmixEM(X, k=2, epsilon = 1e-08, maxit = 1000, maxrestarts=20, 
verb = TRUE, fast=FALSE, ECM = FALSE, arbmean = TRUE, arbvar = TRUE)
plot(mixmdl,which=2)lines(density(X), lty=2, lwd=2)

This is obviously not the best way of doing this. 
The estimates it gives are:
mu 3.6595737(x.2 log()=1.29) 2.3113135(x.1)

These are not too far off but I was wondering if someone knows of a better R 
package/way of doing this or has any hints that would help me coding it from 
scratch (EM+writing down the formulae for ML? but... would this really be 
better than using a more-advanced-already-available-package-made-by-pros?).
The objective is to estimate threshold values at specific FDRs. (some help with 
that would also be most helpful.)

Thanks to all in advance!
To'   
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Fit a mixture of lognormal and normal distributions

2013-03-18 Thread Bert Gunter
You posted this earlier. Do not repost, please. It was seen.

-- Bert

On Mon, Mar 18, 2013 at 3:28 PM, To . . kid...@hotmail.com wrote:


 Hello

 I am trying to find an automated way of fitting a mixture of normal and
 log-normal distributions to data which is clearly bimodal.
 Here's a simulated example:
 x.1-rnorm(6000, 2.4, 0.6)
 x.2-rlnorm(1, 1.3,0.1)
 X-c(x.1, x.2)
 hist(X,100,freq=FALSE, ylim=c(0,1.5))lines(density(x.1), lty=2, lwd=2)
 lines(density(x.2), lty=2, lwd=2)
 lines(density(X), lty=4)

 Currently i am using mixtools and just run:
 library(mixtools)
 mixmdl = normalmixEM(X, k=2, epsilon = 1e-08, maxit = 1000,
 maxrestarts=20, verb = TRUE, fast=FALSE, ECM = FALSE, arbmean = TRUE,
 arbvar = TRUE)
 plot(mixmdl,which=2)lines(density(X), lty=2, lwd=2)

 This is obviously not the best way of doing this.
 The estimates it gives are:
 mu 3.6595737(x.2 log()=1.29) 2.3113135(x.1)

 These are not too far off but I was wondering if someone knows of a better
 R package/way of doing this or has any hints that would help me coding it
 from scratch (EM+writing down the formulae for ML? but... would this really
 be better than using a
 more-advanced-already-available-package-made-by-pros?).
 The objective is to estimate threshold values at specific FDRs. (some help
 with that would also be most helpful.)

 Thanks to all in advance!
 To'
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.




-- 

Bert Gunter
Genentech Nonclinical Biostatistics

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm

[[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] R web application development

2013-03-18 Thread John linux-user
Dear all,

I am wondering if what would be the simple way to develop a simple web 
application that runs R. That is, the web application allows any user upload a 
dataframe as a variable to my web server, a linux-based apache, and then run a 
R package (my package) on the variable that should ideally be handled  as a 
variable in memory instead of saving to the disk in my server , for security 
concern. After running, the result would be returned to the web.  Any 
suggestion will be appreciated. 

Best,

John  
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] R web application development

2013-03-18 Thread Stephen Sefick
I can't offer any advice, but I feel like you could probably get a good 
start on this by looking through the archives.


On Mon 18 Mar 2013 06:55:33 PM CDT, John linux-user wrote:


Dear all,

I am wondering if what would be the simple way to develop a simple web 
application that runs R. That is, the web application allows any user 
upload a dataframe as a variable to my web server, a linux-based 
apache, and then run a R package (my package) on the variable that 
should ideally be handled as a variable in memory instead of saving to 
the disk in my server , for security concern. After running, the 
result would be returned to the web. Any suggestion will be appreciated.


Best,

John
[[alternative HTML version deleted]]



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html

and provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] How to automate this model selection algorithm?

2013-03-18 Thread Andrew Crane-Droesch
I've got a complicated semi-parametric model that I'm fitting with 
mgcv.  I start with a model based on theory.  Its got lots of 
interaction terms.  I want to winnow it down: removing each interaction 
term or un-interacted main effect one by one, checking the AIC, and 
retaining the model that gives me the lowest AIC.  I then want to repeat 
the procedure on the retained model.


Here is a simple example:

set.seed(42)
library(mgcv)
N=220
fac = ceiling(runif(N)*2)
a = rnorm(N); b = rnorm(N); c = rnorm(N); d = runif(N); e = rnorm(N);
y = a^fac + b*e + d/(e+1)

m1 = gam(y~ as.factor(fac)
+s(a)
+s(b)
+s(c)
+s(d)
+te(a,b,c)
+te(a,d,by=as.factor(fac))
)

m2 = gam(y~ as.factor(fac)
+s(a)
+s(b)
+s(c)
+s(d)
+te(b,c)
+te(a,d,by=as.factor(fac))
)

m3 = gam(y~ as.factor(fac)
+s(a)
+s(b)
+s(c)
+s(d)
+te(a,c)
+te(a,d,by=as.factor(fac))
)

m4 = gam(y~ as.factor(fac)
+s(a)
+s(b)
+s(c)
+s(d)
+te(a,b)
+te(a,d,by=as.factor(fac))
)

m5 = gam(y~ as.factor(fac)
+s(a)
+s(b)
+s(c)
+s(d)
+te(a,b,c)
+te(d,by=as.factor(fac))
)

m6 = gam(y~ as.factor(fac)
+s(a)
+s(b)
+s(c)
+s(d)
+te(a,b,c)
+te(a,by=as.factor(fac))
)

m7 = gam(y~ as.factor(fac)
+s(a)
+s(b)
+s(c)
+s(d)
+te(a,b,c)
+te(a,d)
)

selection = AIC(m1,m2,m3,m4,m5,m6,m7)
selection

df  AIC

m1 14.53435 1626.611

m2 12.52501 1622.635

m3 12.54566 1622.615

m4 12.52652 1622.633

m5 13.14108 1620.759

m6 10.99684 1621.156

m7 10.98136 1622.229

m5 is the best

m5 = gam(y~ as.factor(fac)
+s(a)
+s(b)
+s(c)
+s(d)
+te(a,b,c)
+te(d,by=as.factor(fac))
)

m5.1 = gam(y~ as.factor(fac)
+s(a)
+s(b)
+s(c)
+s(d)
+te(b,c)
+te(d,by=as.factor(fac))
)

m5.2 = gam(y~ as.factor(fac)
+s(a)
+s(b)
+s(c)
+s(d)
+te(a,c)
+te(d,by=as.factor(fac))
)
m5.3 = gam(y~ as.factor(fac)
+s(a)
+s(b)
+s(c)
+s(d)
+te(a,b)
+te(d,by=as.factor(fac))
)

m5.4 = gam(y~ as.factor(fac)
+s(a)
+s(b)
+s(c)
+s(d)
+te(a,b,c)
#+te(d,by=as.factor(fac))
)

selection.2 = AIC(m5,m5.1,m5.2,m5.3,m5.4)
selection.2

df  AIC

m5   13.363029 1621.183

m5.1  9.671656 1617.641

m5.2  9.730047 1617.549

m5.3  9.706424 1617.569

m5.4  9.857504 1620.028

5.2 is the best

...and so on.  I'd next try out each interaction term or un-interacted 
main effect in m5.2.  Question is how I can automate this?  To start 
with a model (m1, in my case), and have R give me the best model after 
running this algorithm to the point where there are no longer any 
better models according to this algorithm?


Right now I can do it manually, but as I add data over time, model 
selection may change.


Note the mgcv's select=TRUE functionality doesn't quite work for me 
(at least not directly), because I want to see if single parts of 
three-way interactions can/should be removed.


Thanks in advance for any tips, and apologies for cross-posting (on 
stack overflow).


Best,
Andrew

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 web application development

2013-03-18 Thread John linux-user
Thanks for reply, but which archives? 

Thanks again.
 John



 From: Stephen Sefick ssef...@gmail.com

Cc: r-help@r-project.org r-help@r-project.org 
Sent: Monday, March 18, 2013 8:02 PM
Subject: Re: [R] R web application development

I can't offer any advice, but I feel like you could probably get a good start 
on this by looking through the archives.

On Mon 18 Mar 2013 06:55:33 PM CDT, John linux-user wrote:
 
 Dear all,
 
 I am wondering if what would be the simple way to develop a simple web 
 application that runs R. That is, the web application allows any user upload 
 a dataframe as a variable to my web server, a linux-based apache, and then 
 run a R package (my package) on the variable that should ideally be handled 
 as a variable in memory instead of saving to the disk in my server , for 
 security concern. After running, the result would be returned to the web. Any 
 suggestion will be appreciated.
 
 Best,
 
 John
 [[alternative HTML version deleted]]
 
 
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] melt with complications

2013-03-18 Thread Richard M. Heiberger
thank you all.

the dcast formulation is what I am looking for.  That version of the
formula A+D ~ H
also works with the original cast function.

result.cast - cast(A+D ~ H, value=M, data=meltTest)

It generalizes to a four-factor example, as requested.

This ftable solution didn't do what I wanted because it didn't retain A and
D as factors.

On Mon, Mar 18, 2013 at 6:00 PM, Charles Berry ccbe...@ucsd.edu wrote:

 Richard M. Heiberger rmh at temple.edu writes:

 
  ## Can someone suggest a simpler expression than either of these, with
 the
  goal
  ## of taking a long matrix into a wide one with exactly one of the
 factors
  converted to
  ## columns and all the rest retained as factors.  I want something that
  generalizes beyond
  ## the three factors illustrated here.
 
  ## Rich
 
  meltTest - data.frame(A=rep(c(B,C), each=12),
 D=rep(c(E,F,G), each=4, times=2),
 H=rep(c(I,J,K,L), times=6),
 M=1:24)
 

 amat - ftable( xtabs( M ~ A + D + H, meltTest ),row.vars=1:2 )

 amat is such a matrix with a few attributes added.

 HTH,

 Chuck

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] How to automate this model selection algorithm?

2013-03-18 Thread Andrew Crane-Droesch
Just to add, I've been playing around with select=TRUE in mgcv, and it 
does seem that it could work if I were to specify all of the nested 
two-way interactions in my three-way interactions (see the toy example 
below).  But the problem is that I don't have enough degrees of freedom 
to feed such a model into GAM using my main dataset.


N=200
a = rnorm(N)
b = rnorm(N)
c = rnorm(N)

y = rnorm(N)+a+b+c+a*b

m = gam(y~s(a)+s(b)+s(c)+te(a,b)+te(a,b,c))
msel = gam(y~s(a)+s(b)+s(c)+te(a,b)+te(a,b,c),select=TRUE)
mdrop = gam(y~s(a)+s(b)+s(c)+te(a,b))
summary(m)
summary(msel)
summary(mdrop)
plot(density(m$fitted.values))
lines(density(msel$fitted.values),col=red)
lines(density(mdrop$fitted.values),col=blue)

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 web application development

2013-03-18 Thread David Winsemius
I generally use the Markmail search engine but others use the Newcastle 
resource. Some of other ones can be found with perusal of the Posting Guide and 
Mailing list info page.

http://markmail.org/search/?q=list%3Aorg.r-project.r-help


On Mar 18, 2013, at 5:11 PM, John linux-user wrote:

 Thanks for reply, but which archives? 
 
 Thanks again.
  John
 
 
 
 From: Stephen Sefick ssef...@gmail.com
 
 Cc: r-help@r-project.org r-help@r-project.org 
 Sent: Monday, March 18, 2013 8:02 PM
 Subject: Re: [R] R web application development
 
 I can't offer any advice, but I feel like you could probably get a good start 
 on this by looking through the archives.
 
 On Mon 18 Mar 2013 06:55:33 PM CDT, John linux-user wrote:
 
 Dear all,
 
 I am wondering if what would be the simple way to develop a simple web 
 application that runs R. That is, the web application allows any user upload 
 a dataframe as a variable to my web server, a linux-based apache, and then 
 run a R package (my package) on the variable that should ideally be handled 
 as a variable in memory instead of saving to the disk in my server , for 
 security concern. After running, the result would be returned to the web. 
 Any suggestion will be appreciated.
 
 Best,
 
 John
 [[alternative HTML version deleted]]
 
 
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

David Winsemius
Alameda, CA, USA

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] RGL 3D plots are flat. Please Help

2013-03-18 Thread Noah Silverman
Hello,

I have a matrix of simulated data that I want to plot as a 3D surface using 
RGL.  Have followed the documentation carefully, but my plot only contains a 
weird single slice of the data.  

The actual RGL library and dependance seem fine as all of the demo code plots 
beautifully.

The actual command I'm using is:  rgl.surface(x,z,d, 
color=blue,alpha=0.5,shininess=128)

The length of the x and z vectors match the dimensions of the y matrix.  

Can't figure out why I'm only getting a slice and not a full 3D surface.

Without digging into too much detail, here is what my data looks like.

x
 [1]0  100  200  300  400  500  600  700  800  900 1000 1100 1200 1300 1400 
1500
[17] 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500

z
 [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 
1.8 1.9
[21] 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0

y (Only first bit included here.)
  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
 [1,]000000000 0 0 0 0 0
 [2,]023568   10   11   131416181921
 [3,]036   10   13   16   19   22   262932353842
 [4,]05   10   14   19   24   29   34   384348535862
 [5,]06   13   19   26   32   38   45   515864707783
 [6,]08   16   24   32   40   48   56   6472808896   104
 [7,]0   10   19   29   38   48   58   67   778696   106   115   125
 [8,]0   11   22   34   45   56   67   78   90   101   112   123   134   146
 [9,]0   13   26   38   51   64   77   90  102   115   128   141   154   166
[10,]0   14   29   43   58   72   86  101  115   130   144   158   173   187
[11,]0   16   32   48   64   80   96  112  128   144   160   176   192   208
[12,]0   18   35   53   70   88  106  123  141   158   176   194   211   229
[13,]0   19   38   58   77   96  115  134  154   173   192   211   230   250
[14,]0   21   42   62   83  104  125  146  166   187   208   229   250   270
[15,]0   22   45   67   90  112  134  157  179   202   224   246   269   291
[16,]0   24   48   72   96  120  144  168  192   216   240   264   288   312
[17,]0   26   51   77  102  128  154  179  205   230   256   282   307   333
[18,]0   27   54   82  109  136  163  190  218   245   272   299   326   354
[19,]0   29   58   86  115  144  173  202  230   259   288   317   346   374
[20,]0   30   61   91  122  152  182  213  243   274   304   334   365   395
[21,]0   32   64   96  128  160  192  224  256   288   320   352   384   416
[22,]0   34   67  101  134  168  202  235  269   302   336   370   403   437
[23,]0   35   70  106  141  176  211  246  282   317   352   387   422   458
[24,]0   37   74  110  147  184  221  258  294   331   368   405   442   478
[25,]0   38   77  115  154  192  230  269  307   346   384   422   461   499
[26,]0   40   80  120  160  200  240  280  320   360   400   440   480   520


--
Noah Silverman, M.S.
UCLA Department of Statistics
8117 Math Sciences Building
Los Angeles, CA 90095

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] RGL 3D plots are flat. Please Help

2013-03-18 Thread Pascal Oettli

Hi,

Your example is not reproducible.

What is d in your command?

Regards,
Pascal


On 19/03/13 10:21, Noah Silverman wrote:

Hello,

I have a matrix of simulated data that I want to plot as a 3D surface using RGL.  Have 
followed the documentation carefully, but my plot only contains a weird single 
slice of the data.

The actual RGL library and dependance seem fine as all of the demo code plots 
beautifully.

The actual command I'm using is:  rgl.surface(x,z,d, 
color=blue,alpha=0.5,shininess=128)

The length of the x and z vectors match the dimensions of the y matrix.

Can't figure out why I'm only getting a slice and not a full 3D surface.

Without digging into too much detail, here is what my data looks like.

x
  [1]0  100  200  300  400  500  600  700  800  900 1000 1100 1200 1300 
1400 1500
[17] 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500

z
  [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 
1.8 1.9
[21] 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0

y (Only first bit included here.)
   [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] 
[,14]
  [1,]000000000 0 0 0 0 0
  [2,]023568   10   11   1314161819
21
  [3,]036   10   13   16   19   22   2629323538
42
  [4,]05   10   14   19   24   29   34   3843485358
62
  [5,]06   13   19   26   32   38   45   5158647077
83
  [6,]08   16   24   32   40   48   56   6472808896   
104
  [7,]0   10   19   29   38   48   58   67   778696   106   115   
125
  [8,]0   11   22   34   45   56   67   78   90   101   112   123   134   
146
  [9,]0   13   26   38   51   64   77   90  102   115   128   141   154   
166
[10,]0   14   29   43   58   72   86  101  115   130   144   158   173   187
[11,]0   16   32   48   64   80   96  112  128   144   160   176   192   208
[12,]0   18   35   53   70   88  106  123  141   158   176   194   211   229
[13,]0   19   38   58   77   96  115  134  154   173   192   211   230   250
[14,]0   21   42   62   83  104  125  146  166   187   208   229   250   270
[15,]0   22   45   67   90  112  134  157  179   202   224   246   269   291
[16,]0   24   48   72   96  120  144  168  192   216   240   264   288   312
[17,]0   26   51   77  102  128  154  179  205   230   256   282   307   333
[18,]0   27   54   82  109  136  163  190  218   245   272   299   326   354
[19,]0   29   58   86  115  144  173  202  230   259   288   317   346   374
[20,]0   30   61   91  122  152  182  213  243   274   304   334   365   395
[21,]0   32   64   96  128  160  192  224  256   288   320   352   384   416
[22,]0   34   67  101  134  168  202  235  269   302   336   370   403   437
[23,]0   35   70  106  141  176  211  246  282   317   352   387   422   458
[24,]0   37   74  110  147  184  221  258  294   331   368   405   442   478
[25,]0   38   77  115  154  192  230  269  307   346   384   422   461   499
[26,]0   40   80  120  160  200  240  280  320   360   400   440   480   520


--
Noah Silverman, M.S.
UCLA Department of Statistics
8117 Math Sciences Building
Los Angeles, CA 90095

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] RGL 3D plots are flat. Please Help

2013-03-18 Thread Noah Silverman
Oops, that was a type.

d is just y.

If you want to reproduce, I put the entire matrix, as a csv here:   
http://pastebin.com/gniyD4Rc


Thanks,


--
Noah Silverman, M.S.
UCLA Department of Statistics
8117 Math Sciences Building
Los Angeles, CA 90095

On Mar 18, 2013, at 6:58 PM, Pascal Oettli kri...@ymail.com wrote:

 Hi,
 
 Your example is not reproducible.
 
 What is d in your command?
 
 Regards,
 Pascal
 
 
 On 19/03/13 10:21, Noah Silverman wrote:
 Hello,
 
 I have a matrix of simulated data that I want to plot as a 3D surface using 
 RGL.  Have followed the documentation carefully, but my plot only contains a 
 weird single slice of the data.
 
 The actual RGL library and dependance seem fine as all of the demo code 
 plots beautifully.
 
 The actual command I'm using is:  rgl.surface(x,z,d, 
 color=blue,alpha=0.5,shininess=128)
 
 The length of the x and z vectors match the dimensions of the y matrix.
 
 Can't figure out why I'm only getting a slice and not a full 3D surface.
 
 Without digging into too much detail, here is what my data looks like.
 
 x
  [1]0  100  200  300  400  500  600  700  800  900 1000 1100 1200 1300 
 1400 1500
 [17] 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500
 
 z
  [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 
 1.8 1.9
 [21] 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0
 
 y (Only first bit included here.)
   [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] 
 [,14]
  [1,]000000000 0 0 0 0   
   0
  [2,]023568   10   11   1314161819   
  21
  [3,]036   10   13   16   19   22   2629323538   
  42
  [4,]05   10   14   19   24   29   34   3843485358   
  62
  [5,]06   13   19   26   32   38   45   5158647077   
  83
  [6,]08   16   24   32   40   48   56   6472808896   
 104
  [7,]0   10   19   29   38   48   58   67   778696   106   115   
 125
  [8,]0   11   22   34   45   56   67   78   90   101   112   123   134   
 146
  [9,]0   13   26   38   51   64   77   90  102   115   128   141   154   
 166
 [10,]0   14   29   43   58   72   86  101  115   130   144   158   173   
 187
 [11,]0   16   32   48   64   80   96  112  128   144   160   176   192   
 208
 [12,]0   18   35   53   70   88  106  123  141   158   176   194   211   
 229
 [13,]0   19   38   58   77   96  115  134  154   173   192   211   230   
 250
 [14,]0   21   42   62   83  104  125  146  166   187   208   229   250   
 270
 [15,]0   22   45   67   90  112  134  157  179   202   224   246   269   
 291
 [16,]0   24   48   72   96  120  144  168  192   216   240   264   288   
 312
 [17,]0   26   51   77  102  128  154  179  205   230   256   282   307   
 333
 [18,]0   27   54   82  109  136  163  190  218   245   272   299   326   
 354
 [19,]0   29   58   86  115  144  173  202  230   259   288   317   346   
 374
 [20,]0   30   61   91  122  152  182  213  243   274   304   334   365   
 395
 [21,]0   32   64   96  128  160  192  224  256   288   320   352   384   
 416
 [22,]0   34   67  101  134  168  202  235  269   302   336   370   403   
 437
 [23,]0   35   70  106  141  176  211  246  282   317   352   387   422   
 458
 [24,]0   37   74  110  147  184  221  258  294   331   368   405   442   
 478
 [25,]0   38   77  115  154  192  230  269  307   346   384   422   461   
 499
 [26,]0   40   80  120  160  200  240  280  320   360   400   440   480   
 520
 
 
 --
 Noah Silverman, M.S.
 UCLA Department of Statistics
 8117 Math Sciences Building
 Los Angeles, CA 90095
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Convert to date and time of the year

2013-03-18 Thread Janesh Devkota
Dear R Users, 

I have data for more than 3 years. For each year I want to find the day
corresponding to Jaunary 1 of that year. For example:

 x - c('5/5/2007','12/31/2007','1/2/2008')

 #Convert to day of year (julian date) -

 strptime(x,%m/%d/%Y)$yday+1

[1] 125 365   2

I want to know how to do the same thing but with time added. But I still get
the day not time. Can anyone suggest what is the better way to find the
julian date with date and time ?

 x1 - c('5/5/2007 02:00','12/31/2007 05:58','1/2/2008 16:25')

 #Convert to day of year (julian date) -

 strptime(x1,%m/%d/%Y %H:%M)$yday+1

[1] 125 365   2

Thank you so much.

Janesh


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Convert to date and time of the year

2013-03-18 Thread Ben Tupper
Hi Janesh,

On Mar 18, 2013, at 10:49 PM, Janesh Devkota wrote:

 Dear R Users, 
 
 I have data for more than 3 years. For each year I want to find the day
 corresponding to Jaunary 1 of that year. For example:
 
 x - c('5/5/2007','12/31/2007','1/2/2008')
 
 #Convert to day of year (julian date) -
 
 strptime(x,%m/%d/%Y)$yday+1
 
 [1] 125 365   2
 
 I want to know how to do the same thing but with time added. But I still get
 the day not time. Can anyone suggest what is the better way to find the
 julian date with date and time ?
 
 x1 - c('5/5/2007 02:00','12/31/2007 05:58','1/2/2008 16:25')
 
 #Convert to day of year (julian date) -
 
 strptime(x1,%m/%d/%Y %H:%M)$yday+1
 
 [1] 125 365   2
 

Would this do it for you?

x1 - c('5/5/2007 02:00','12/31/2007 05:58','1/2/2008 16:25')
x1 - as.POSIXct(x1, format = %m/%d/%Y %H:%M)

day - as.numeric(format(x1, %j)) 
hour - as.numeric(format(x1, %H))/24 
minute - as.numeric(format(x1, %M))/(60*24) 
second - as.numeric(format(x1, %S))/(60*60*24)

day + hour + minute + second

[1] 125.08 365.248611   2.684028

Cheers,
Ben




 Thank you so much.
 
 Janesh
 
 
   [[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.

Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.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.


  1   2   >