Re: [R] How to Plot With Different Marker ( ‘x ’ and ‘o’) Based on Condition in R

2010-07-06 Thread Peter Ehlers

On 2010-07-05 23:05, Gundala Viswanath wrote:

Dear Expert,

I have a data that looks like this:

for_y_axis-c(0.49534,0.80796,0.93970,0.8)
for_x_axis-c(1,2,3,4)
count-c(0,33,0,4)

What I want to do is to plot the graph using for_x_axis and
for_y_axis but will mark
each point with o if the value is equal to 0(zero) and with x if
the count value is greater than zero.

Is there a simple way to achieve that in R?


Here's one way:

 .pch - ifelse(count  0, x, o)
 plot(for_x_axis, for_y_axis, pch = .pch)

or you might prefer

 plot(for_x_axis, for_y_axis, pch = .pch, family = mono)


  -Peter Ehlers

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

2010-07-06 Thread Rubén Roa
 -Mensaje original-
 De: r-help-boun...@r-project.org 
 [mailto:r-help-boun...@r-project.org] En nombre de Fabian Gehring
 Enviado el: lunes, 05 de julio de 2010 21:53
 Para: r-help@r-project.org
 Asunto: [R] selection of optim parameters
 
 Hi all,
 
 I am trying to rebuild the results of a study using a 
 different data set. I'm using about 450 observations. The 
 code I've written seems to work well, but I have some 
 troubles minimizing the negative of the LogLikelyhood 
 function using 5 free parameters.
 
 As starting values I am using the result of the paper I am rebuiling. 
 The system.time of the calculation of the function is about 0.65 sec. 
 Since the free parameters should be within some boundaries I 
 am using the following command:
 
 optim(fn=calculateLogLikelyhood, c(0.4, 2, 0.4, 8, 0.8), 
 lower=c(0,0,0,0,0), upper=c(1, 5, Inf, Inf, 1), control=list(trace=1,
 maxit=1000))
 
 Unfortunately the result doesn't seem to be reasonable. 3 of 
 the optimized parameters are on the boundaries.

1) Your parameters seem to vary over several orders of magnitude. The control 
argument has a parscale parameter that can be used to re-scale all parameters 
to the same order of magnitude. Alternatively, your could estimate the log of 
your parameters, say par=c(log(0.4), log(2), log(0.4), log(8), log(0.8) 
(and equivalent changes in lower and upper), and in your function, instead of 
the parameter value, use exp(parameter value9. That way the _numerical 
optimization_ occurs in the log space whereas the _function evaluation_ occurs 
in the original space. This transformation approach would make your parameter 
estimates and their hessian matrix (in case you are interested in the hessian) 
be output in the transformed space, so estimates and their covariance matrix 
will have to be back-transformed. For estimates just use exp(), whereas for the 
covariance matrix you might have to use something like Taylor series.
2) Did you use L-BFGS-B in the method argument of optim()? This method admits 
box-constrained optimization whereas the default (which you seem to be using, 
Nelder-Mead) in unconstrained, as far as I know.

 Unfortunately I don't have much experience using 
 optimizatzion methods. 
 That's why I am asking you.
 Do you have any hints for me what should be taken into 
 account when doing such an optimization.
 
 Is there a good way to implement the boundaries into the code 
 (instead of doing it while optimizing)? I've read about 
 parscale in the help-section. Unfortunately I don't really 
 know how to use it. And anyways, could this help? What other 
 points/controls should be taken into account?

About using parscale, see Ravi Varadhan's post Re: optim() not finding optimal 
values the 27th of June.

HTH

Rubén



Dr. Rubén Roa-Ureta
AZTI - Tecnalia / Marine Research Unit
Txatxarramendi Ugartea z/g
48395 Sukarrieta (Bizkaia)
SPAIN

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Conditional Splitting a Vectors into Two Vectors

2010-07-06 Thread Gundala Viswanath
Suppose I have two vectors of same dimensions:

x -c(0.49534,0.80796,0.93970,0.8)
count  -c(0,33,0,4)

How can I group the vectors 'x' into two vectors:

 1. Vector `grzero` that contain value in x with `count` value greater
than 0 and
 2. Vector `eqzero` with value in x with `count` value equal to zero.

Yielding

 print(grzero)
 [1] 0.80796 0.8
 print(eqzero)
 [1] 0.49534 0.93970

Regards,
G.V.

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

2010-07-06 Thread Peter Ehlers

On 2010-07-02 13:23, Erin Hodgess wrote:

Dear R People:

I have found the starma.c program in R.

But now I need to find the R_setup_starma and R_free_starma subroutines as well.

Where would I go about finding them, please?


These are objects of class NativeSymbolInfo. You can
see their list elements with

 getNativeSymbolInfo(setup_starma)
 getNativeSymbolInfo(free_starma)

or with

 stats:::R_setup_starma


  -Peter Ehlers



Thanks for any help.

Have a great weekend.

Sincerely,
Erin


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Conditional Splitting a Vectors into Two Vectors

2010-07-06 Thread Peter Ehlers

On 2010-07-06 0:39, Gundala Viswanath wrote:

Suppose I have two vectors of same dimensions:

 x-c(0.49534,0.80796,0.93970,0.8)
 count-c(0,33,0,4)

How can I group the vectors 'x' into two vectors:

  1. Vector `grzero` that contain value in x with `count` value greater
than 0 and
  2. Vector `eqzero` with value in x with `count` value equal to zero.

Yielding

   print(grzero)
   [1] 0.80796 0.8
   print(eqzero)
   [1] 0.49534 0.93970

Regards,
G.V.



It might be time to work your way through
'An Introduction to R'.

Anyway, here you want conditional extraction:

 x[count  0]
 x[!(count  0)]  ##or equivalent


  -Peter Ehlers

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


Re: [R] Function to compute the multinomial beta function?

2010-07-06 Thread Robin Hankin

It's usually better to build vectorization in to functions:

 beta3- function (n1, n2, n3) 
exp(lgamma(n1)+lgamma(n2)+lgamma(n3)-lgamma(n1+n2+n3))

 f - function(x){exp(sum(lgamma(x))-lgamma(sum(x)))}
 beta3(5,3,8)
[1] 1.850002e-07
 f(c(5,3,8))
[1] 1.850002e-07


rksh


On 07/06/2010 01:54 AM, Robert A LaBudde wrote:

At 05:10 PM 7/5/2010, Gregory Gentlemen wrote:

Dear R-users,

Is there an R function to compute the multinomial beta function? That 
is, the normalizing constant that arises in a Dirichlet distribution. 
For example, with three parameters the beta function is 
Beta(n1,n2,n2) = Gamma(n1)*Gamma(n2)*Gamma(n3)/Gamma(n1+n2+n3)


 beta3- function (n1, n2, n3) 
exp(lgamma(n1)+lgamma(n2)+lgamma(n3)-lgamma(n1+n2+n3))

 beta3(5,3,8)
[1] 1.850002e-07




--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877

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

2010-07-06 Thread nuncio m
Hi list,
  I have the following code to compute the acf of a time series
acfresid - acf(residfit), where residfit is the series
when I type acfresid at the prompt the follwoing is displayed

Autocorrelations of series ‘residfit’, by lag

0. 0.0833 0.1667 0.2500 0. 0.4167 0.5000 0.5833 0.6667 0.7500 0.8333

 1.000 -0.015  0.010  0.099  0.048 -0.014 -0.039 -0.019  0.040  0.018  0.042

0.9167 1. 1.0833 1.1667 1.2500 1. 1.4167 1.5000 1.5833 1.6667 1.7500

 0.078 -0.029  0.028 -0.016 -0.021 -0.109  0.000 -0.038 -0.006  0.015 -0.032

1.8333 1.9167 2. 2.0833
-0.002  0.014 -0.226 -0.030
Residfit is a timeseries object at monthly interval (0.0833), Here I
understand R computed the correlation at lags 0 to 2 years.

What is surprising to me is
if I type acfresidfit at the prompt the following is displayed

Autocorrelations of series ‘residfit’, by lag

 0  1  2  3  4  5  6  7  8  9 10

 1.000 -0.004  0.011  0.041 -0.056  0.019 -0.052 -0.027 -0.008 -0.012 -0.034

11 12 13 14 15 16 17 18 19 20 21

 0.024 -0.005  0.006 -0.045  0.031 -0.035 -0.011 -0.021 -0.020 -0.010 -0.007

22 23 24 25
-0.038  0.017  0.051  0.038
From the header I understand both are autocorrelation computed at the same
lags. but the correlations are different

where am I going wrong and which is the correct one.

file residfit is also attached(filename-fileree2_test_out.txt)
Thanks
nuncio
-- 
Nuncio.M
Research Scientist
National Center for Antarctic and Ocean research
Head land Sada
Vasco da Gamma
Goa-403804
4.54540234232334
-14.4778008999506
-3.79668140611868
-7.81347830768482
-6.27293225798647
-6.87201981207487
-6.64965905122317
-6.75123982158051
18.7798275931915
6.81254237499438
17.533220743665
11.8179723199377
-13.8382453401278
3.54961036332585
-3.97203313203956
-1.11029123677042
-2.21102643268432
-1.78925830375029
-1.95531252764598
-1.90169073408017
-16.2251504389829
31.1513180684656
-5.0374856234746
3.64237537817929
35.7825921539956
-24.3131812339976
7.67437033129054
-4.57640838126481
0.13663691723
-1.678854866888
-0.980572231720623
-1.25035922801876
-23.7161201125191
5.99944860265195
6.61868629783817
-1.3221234388819
10.2414420460275
-8.47848871031044
1.63078138547826
-2.26306522145922
-0.764479820758242
-1.34245520108361
-1.12077197260524
-1.20702766721386
5.98843679322213
7.33766331906203
-2.39623468909737
0.0189595505058051
-9.7878766373797
3.90933597855122
-3.15027883357282
-0.432594130850131
-1.48003000602628
-1.07756318111729
-1.23343646676257
-1.17429968377673
13.8000693717299
6.02003887937457
-8.34195666853458
-1.11099261410632
1.56296156653197
-2.75269349211589
-0.595469082283298
-1.42707843208285
-1.10772380913903
-1.23159158432425
-1.18477973889978
-1.20369456481332
21.7239849952369
-9.66196813638387
-0.891521279149956
9.25210349671729
10.3427639242166
-9.80695249149604
2.11038354823627
-2.47972753146567
-0.713019047280778
-1.39424540580337
-1.13280146608686
31.4656314670951
-0.108909595891454
-0.451065307346077
10.56546397663
-0.647946708663838
18.1367401912301
-3.53869009989366
-2.65626056469235
-0.654807910127964
-1.42643142842822
-1.13017662852014
-1.24514893937744
-13.9208705459327
-5.42722370392831
11.5554183170163
-1.11287843661796
-3.79345702201642
-16.8470908324909
4.77680059742864
-2.61594762546361
-0.680989698404531
-1.42700689601552
-1.14061278256847
-1.25178786153672
-19.2156665052095
3.59214434723623
9.68031595123078
-0.905773257642835
11.6697859140499
-10.438203868059
-3.04292317412673
0.76333844274339
-1.99296992395689
-0.932436833094549
-1.34172353169702
-1.18500028791065
-1.2462406728489
-33.1551884101196
19.0913314431478
-6.7817942558301
5.00083912161537
-3.23440161285666
-0.80515939314185
-1.3987298644782
-1.17104115381219
-1.25960948746092
-1.22639103288181
-1.2400712072929
46.5643087168762
17.4522540549196
-10.9812432669912
10.0524952224354
-1.45992013300943
4.90518728893209
-5.35388787161763
0.342273064087802
-1.85213822475519
-1.00798494789867
-1.33394637973212
-1.20931069510199
-19.8506519887964
-28.6308445373069
5.59082974922492
6.79381059715144
-2.5896614051167
-1.29096368912288
-1.29288147655713
-1.2322294690926
-1.25647397513093
-1.2480258537729
-1.25216727045333
-1.25146060728689
-27.5730296981935
5.19481715163884
5.228853569497
-3.52768178324466
2.26390585067962
10.4196518274518
-8.93836873407259
1.70129806057068
-2.39679794286931
-0.819558553443901
-1.42782263987167
-1.19447559456423
-1.28522288490369
-14.7482003454544
14.9839484682792
-3.23392273754434
3.74619714728151
-11.0770500672262
4.03870484048016
-3.30667482507862
-0.478945510500303
-1.56875824361717
-1.14997259983422
-1.31213006657225
44.8494266790537
45.874661297
-25.3865845300823
11.7813342281731
-1.19733031473043
17.5712159736525
-13.2582237918684
3.34350682483415
-3.050508327763
-0.589138057355864
-1.53787032402965
-1.17341307292178
-8.8458687102
14.3506481320503
-4.07628087468823

[R] Error in affypdnn package

2010-07-06 Thread Hasan, Ahmed Ryadh - hasar001
Dear all,
I am a PhD student working with Affymetrix HGU133atag array for analyzing the 
Latin square experiment.

I was trying to generate gene expression index for hgu133atag array for PDNN 
model. While extracting the chiptype specific data structure, I got the 
following error-


 library(affypdnn)
Loading required package: affy
Loading required package: Biobase
Welcome to Bioconductor
  Vignettes contain introductory material. To view, type
  'openVignette()'. To cite Bioconductor, see
  'citation(Biobase)' and for packages 'citation(pkgname)'.
registering new summary method 'pdnn'.
registering new pmcorrect method 'pdnn' and 'pdnnpredict'.
 library(hgu133atagprobe)
Loading required package: AnnotationDbi
 energy.file - system.file(exampleData, 
 pdnn-energy-parameter_hg-u133a.txt, package = affypdnn)
 params.chiptype - pdnn.params.chiptype(energy.file, probes.pack = 
 hgu133aprobe)
Calculating chip type specific parameters, (may take some time)...
||
|Error in object$call : $ operator not defined for this S4 class
 
After this, I also tried to do it for hgu95av2 array as described in the 
package manual-

 library(affypdnn)
Loading required package: affy
Loading required package: Biobase
Welcome to Bioconductor
  Vignettes contain introductory material. To view, type
  'openVignette()'. To cite Bioconductor, see
  'citation(Biobase)' and for packages 'citation(pkgname)'.
registering new summary method 'pdnn'.
registering new pmcorrect method 'pdnn' and 'pdnnpredict'.
 library(hgu95av2probe)
Loading required package: AnnotationDbi
 energy.file - system.file(exampleData, 
 pdnn-energy-parameter_hg-u95av2.txt, package = affypdnn)
 params.chiptype - pdnn.params.chiptype(energy.file, probes.pack = 
 hgu95av2probe)
Calculating chip type specific parameters, (may take some time)...
||
|Error in object$call : $ operator not defined for this S4 class
 
Please help me to get rid of this error(Error in object$call : $ operator not 
defined for this S4 class). 

Best regards,

Ahmed Ryadh Hasan

PhD Student,
Knowledge-Based Intelligent Engineering Systems centre (KES)
Mawson Lakes Campus, University of South Australia,
Mawson Lakes, SA 5095

Phone: +61 8 8302 8332 (Office)
   +61 8 8359 5747 (Home)
   +61 4 3006 4390 (Cell)

[[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] Assign Formulas to Arrays or Matrices?

2010-07-06 Thread McLovin

Hi,

I am very new to R.  I am hoping to create formulas and assign them to
locations within an array (or matrix, if it will work).

Here's a simplified example of what I'm trying to do:

form.arr - array(31,5,3)
for (i in seq(from=1, to=31, by=1)) {
for (j in seq(from=1, to=5, by=1)) {
form.arr[i,j,] - as.formula(y~1+2)
}
}

which results in this error:
Error in form.arr[i, j, ] - as.formula(y ~ 1 + 2) : 
  incorrect number of subscripts

The reason I had made the 3rd dimension of the array size 3 is because
that's the length R tells me that formula is.

When I had tried to do this using a matrix, using this code:

form.mat - matrix(31,5,3)
for (i in seq(from=1, to=31, by=1)) {
for (j in seq(from=1, to=5, by=1)) {
form.mat[i,j] = as.formula(y~1+2)
}
}

I was told:

Error in form.mat[i, j] = as.formula(y ~ 1 + 2) : 
  number of items to replace is not a multiple of replacement length

My question is: is it possible to assign formulas within a matrix or array? 
If so, how?  thanks@real.com
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Assign-Formulas-to-Arrays-or-Matrices-tp2279136p2279136.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Assign Formulas to Arrays or Matrices?

2010-07-06 Thread Peter Ehlers

On 2010-07-06 1:13, McLovin wrote:


Hi,

I am very new to R.  I am hoping to create formulas and assign them to
locations within an array (or matrix, if it will work).

Here's a simplified example of what I'm trying to do:

form.arr- array(31,5,3)
for (i in seq(from=1, to=31, by=1)) {
for (j in seq(from=1, to=5, by=1)) {
form.arr[i,j,]- as.formula(y~1+2)
}
}

which results in this error:
Error in form.arr[i, j, ]- as.formula(y ~ 1 + 2) :
   incorrect number of subscripts

The reason I had made the 3rd dimension of the array size 3 is because
that's the length R tells me that formula is.

When I had tried to do this using a matrix, using this code:

form.mat- matrix(31,5,3)
for (i in seq(from=1, to=31, by=1)) {
for (j in seq(from=1, to=5, by=1)) {
form.mat[i,j] = as.formula(y~1+2)
}
}

I was told:

Error in form.mat[i, j] = as.formula(y ~ 1 + 2) :
   number of items to replace is not a multiple of replacement length

My question is: is it possible to assign formulas within a matrix or array?
If so, how?  thanks@real.com


I don't think it's possible in the way you're trying
to do it. A formula is not the same thing as a 3-element
vector.

Why not use a list?

ps. for (i in seq(from=1, to=31, by=1)) is equivalent
to for(i in seq_len(31))

  -Peter Ehlers

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Selection with changing number of columns

2010-07-06 Thread Kunzler, Andreas
Dear list,

I'm looking for a way to select rows of a data.frame with changing number of 
columns (constraint) involved.

Assume a data (d) structure like


Var.1 Var.2 Var.3
9   2   1
2   9   5
1   2   1

I know the number of involved columns.

Is there a way to generate the following selection automatically (maybe for 
loop), so that it makes no difference if there are two or ten columns involved.
 
Selection:
d[d$Var.1==9 | d$Var.1==9 | d$Var.1==9  ,]


Does anybody know a way?

Thanks

Mit freundlichen Grüßen

Andreas Kunzler

Bundeszahnärztekammer (BZÄK)
Chausseestraße 13
10115 Berlin

Tel.: 030 40005-113
Fax:  030 40005-119

E-Mail: a.kunz...@bzaek.de 

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

2010-07-06 Thread Joris Meys
Your question is a bit confusing. acfresidfit is an object, of which
we don't know the origin. with your test file, I arrive at the first
correlations (but with integer headings) :

 residfit - read.table(fileree2_test_out.txt)
 acf(residfit)
 acfresid - acf(residfit)
 acfresid

Autocorrelations of series ‘residfit’, by lag

 0  1  2  3  4  5  6  7  8  9
   10 11 12 13 14 15 16
 1.000 -0.015  0.010  0.099  0.048 -0.014 -0.039 -0.019  0.040  0.018
0.042  0.078 -0.029  0.028 -0.016 -0.021 -0.109
17 18 19 20 21 22 23 24 25
 0.000 -0.038 -0.006  0.015 -0.032 -0.002  0.014 -0.226 -0.030

Could you please check where the object acfresidfit is coming from and
how you generated it?
Cheers
Joris

On Tue, Jul 6, 2010 at 9:47 AM, nuncio m nunci...@gmail.com wrote:
 Hi list,
          I have the following code to compute the acf of a time series
 acfresid - acf(residfit), where residfit is the series
 when I type acfresid at the prompt the follwoing is displayed

 Autocorrelations of series ‘residfit’, by lag

 0. 0.0833 0.1667 0.2500 0. 0.4167 0.5000 0.5833 0.6667 0.7500 0.8333

  1.000 -0.015  0.010  0.099  0.048 -0.014 -0.039 -0.019  0.040  0.018  0.042

 0.9167 1. 1.0833 1.1667 1.2500 1. 1.4167 1.5000 1.5833 1.6667 1.7500

  0.078 -0.029  0.028 -0.016 -0.021 -0.109  0.000 -0.038 -0.006  0.015 -0.032

 1.8333 1.9167 2. 2.0833
 -0.002  0.014 -0.226 -0.030
 Residfit is a timeseries object at monthly interval (0.0833), Here I
 understand R computed the correlation at lags 0 to 2 years.

 What is surprising to me is
 if I type acfresidfit at the prompt the following is displayed

 Autocorrelations of series ‘residfit’, by lag

     0      1      2      3      4      5      6      7      8      9     10

  1.000 -0.004  0.011  0.041 -0.056  0.019 -0.052 -0.027 -0.008 -0.012 -0.034

    11     12     13     14     15     16     17     18     19     20     21

  0.024 -0.005  0.006 -0.045  0.031 -0.035 -0.011 -0.021 -0.020 -0.010 -0.007

    22     23     24     25
 -0.038  0.017  0.051  0.038
 From the header I understand both are autocorrelation computed at the same
 lags. but the correlations are different

 where am I going wrong and which is the correct one.

 file residfit is also attached(filename-fileree2_test_out.txt)
 Thanks
 nuncio
 --
 Nuncio.M
 Research Scientist
 National Center for Antarctic and Ocean research
 Head land Sada
 Vasco da Gamma
 Goa-403804

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





-- 
Joris Meys
Statistical consultant

Ghent University
Faculty of Bioscience Engineering
Department of Applied mathematics, biometrics and process control

tel : +32 9 264 59 87
joris.m...@ugent.be
---
Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] metafor and meta-analysis at arm-level

2010-07-06 Thread Angelo Franchini
Hello Wolfgang,

Thank you very much for your response.
When you mentionthe appropriate design matrix, do you mean by that the
'n1i, n2i, m1i, m2i, sd1i, sd2i' arguments of the rma function, or am I
missing something?
I read the documentation on metafor (introduction), rma/rma.uni and
escalc, and that was the only way that I could find for the package to use
information at the arm-level rather than the trial one.

As for the complexity of possible correlations between effects, that is
something to be considered for the network analysis case, correct?

Many thanks.

Best regards,
Angelo



On Sun, July 4, 2010 6:06 am, Wolfgang Viechtbauer wrote:
 Hello Angelo,

 You can either supply the arm-level outcomes and corresponding sampling
 variances directly (via the yi and vi arguments) or supply the necessary
 information so that the escalc() or rma() functions can calculate an
 appropriate arm-level outcome (such as the log odds). See the
 documentation of the escalc() function and in particular the part about
 proportions and tranaformations thereof as possible outcome measures.

 This is the easy part. Then you need to set up an appropriate design
 matrix to code what arm each observed outcome corresponds to. And finally
 comes the tricky/problematic part. The rma() function assumes independent
 sampling errors and independent random effects for each observed outcome.
 Independent sampling errors is (usually) ok when using arm-level outcomes,
 but the independent random errors part may not be appropriate. This is why
 I am working on functions that do not make this independence assumption.
 With those functions, you can then carry out multivariate and network-type
 meta-analyses. These functions will become part of the metafor package in
 the future.

 Best,

 --
 Wolfgang Viechtbauer
 http://www.wvbauer.com

 Angelo Franchini angelo.franch...@bristol.ac.uk wrote:

Hi,

I have been looking for an R package which allowed to do meta-analysis
(both pairwise and network/mixed-treatment) at arm-level rather than at
trial-level, the latter being the common way in which meta-analysis is
done.
By arm-level meta-analysis I mean one that accounts for data provided at
the level of the individual arms of each trial and that does not simply
derive the difference between arms and do the meta-analysis on that.

I am not sure metafor can do that, but hopefully someone more experienced
on it can clarify that to me. I can see that it can take data in both
forms, arm and trial level, but it looks as if the arm-level information
would be converted into trial one through escalc and the latter then used
for the meta-analysis. Is that right?

Many thanks.

Angelo


--
NIHR Research Methods Training Fellow,
Department of Community Based Medicine
University of Bristol
25 Belgrave Road
Bristol BS8 2AA

Tel. 0779 265-6552

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




-- 
NIHR Research Methods Training Fellow,
Department of Community Based Medicine
University of Bristol
25 Belgrave Road
Bristol BS8 2AA

Tel. 0779 265-6552

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] multiple time series plot with dual 'y' axes

2010-07-06 Thread Jorge Hernandez
Hello.


I would like to know how to generate dual 'y' axes on a multiple time 
series plot. 


I am using  ts.plot()  to get the multiple time series plot, but I would 
like a second 
vertical axis on the right to include another time series on a different 
scale.


Thanks for any help.

Cheers.

Jorge
**

PRIVATE  CONFIDENTIAL: This email and any files transmitted with it are 
confidential. Any unauthorized use of the information contained in t
his email or its attachments is prohibited.  If this email is received in 
error, please contact the sender and delete the material from your
computer systems. Do not use, copy, or disclose the contents of this email or 
any attachments.

Abu Dhabi Investment Authority (ADIA) does not enter into contracts or provide 
undertakings by email.  ADIA accepts no responsibility for the
 content of this mail to the extent that it is unrelated to its activities or 
the same consists of statements or opinions  which are the send
er's own and not made on behalf of ADIA.  ADIA does not accept any liability 
for any errors or omissions in the content of this email caused
by electronic and technical failures. Although ADIA has taken reasonable 
precautions to ensure that no viruses are present in this email, ADI
A accepts no responsibility for any loss or damage arising from the use of this 
email or its attachments.

**

[[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] timeseries

2010-07-06 Thread Joris Meys
Difficult to guess why, but I reckon you should use ts() instead of
as.ts. Otherwise set the tsp-attribute correctly. Eg :

 x - cumsum(1+round(rnorm(20),2))
 as.ts(x)
Time Series:
Start = 1
End = 20
Frequency = 1
 [1]  0.87  3.51  4.08  4.20  3.25  4.63  6.30  6.89  9.28  9.93 10.19
 9.97 10.20 11.51 11.52 11.53 12.97 14.04 17.02 18.47

 tseries - ts(x,frequency=12,start=c(2004,3))
 tseries
   Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
2004  0.87  3.51  4.08  4.20  3.25  4.63  6.30  6.89  9.28  9.93
2005 10.19  9.97 10.20 11.51 11.52 11.53 12.97 14.04 17.02 18.47

 tsp(x) - c(2004+2/12,2005.75,12)
 as.ts(x)
   Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
2004  0.87  3.51  4.08  4.20  3.25  4.63  6.30  6.89  9.28  9.93
2005 10.19  9.97 10.20 11.51 11.52 11.53 12.97 14.04 17.02 18.47

See ?ts to get the options right. I suggest to use the function ts()
instead of assigning the tsp attribute.
)
Cheers
Joris


On Mon, Jul 5, 2010 at 9:35 AM, nuncio m nunci...@gmail.com wrote:
 Dear useRs,
 I am trying to construct a time series using as.ts function, surprisingly
 when I plot
 the data the x axis do not show the time in years, however if I use
 ts(data), time in years are shown in the
 x axis.  Why such difference in the results of both the commands
 Thanks
 nuncio


 --
 Nuncio.M
 Research Scientist
 National Center for Antarctic and Ocean research
 Head land Sada
 Vasco da Gamma
 Goa-403804

        [[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.




-- 
Joris Meys
Statistical consultant

Ghent University
Faculty of Bioscience Engineering
Department of Applied mathematics, biometrics and process control

tel : +32 9 264 59 87
joris.m...@ugent.be
---
Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] calculation on series with different time-steps

2010-07-06 Thread Joris Meys
Look at ?ifelse en ?abs, eg :

data_frame$new_column_in_dataframe - ifelse(stage$julian_day ==
baro$julian_day  abs(stage$time -
baro$hour) = 30,
stage$stage.cm - baro$pressure, NA )

Cheers
Joris


On Thu, Jul 1, 2010 at 10:09 PM, Jeana Lee jeana@colorado.edu wrote:
 Hello,

 I have two series, one with stream stage measurements every 5 minutes, and
 the other with barometric pressure measurements every hour.  I want to
 subtract each barometric pressure measurement from the 12 stage measurements
 closest in time to it (6 stage measurements on either side of the hour).

 I want to do something like the following, but I don't know the syntax.

 If the Julian day of the stage measurement is equal to the Julian day of
 the pressure measurement, AND the absolute value of the difference between
 the time of the stage measurement and the hour of the pressure measurement
 is less than or equal to 30 minutes, then subtract the pressure measurement
 from the stage measurement (and put it in a new column in the stage data
 frame).

                     if ( stage$julian_day = baro$julian_day  |stage$time -
 baro$hour| = 30 )
                     then (stage$stage.cm - baro$pressure)

 Can you help me?

 Thanks,
 JL

        [[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.




-- 
Joris Meys
Statistical consultant

Ghent University
Faculty of Bioscience Engineering
Department of Applied mathematics, biometrics and process control

tel : +32 9 264 59 87
joris.m...@ugent.be
---
Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php

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

2010-07-06 Thread n.via...@libero.it

Dear list,
I'm trying to generate a latex Document in which there are a lot of tables. I'm 
using the Sweave function in the package utils, but I'm having a lot of 
problems with the format. This is my code:
\documentclass[a4paper]{amsbook}
\title{Schema di bilancio}
\begin{document}
\maketitle
echo=F,results=hide=
report=Bilanci
mynames-names(report)
mynames[mynames==AA01]-Immobilizzazioni tecniche nette
mynames[mynames==AA01I]-Immobilizzazioni imm. nette
mynames[mynames==AA01M]-Immobilizzazioni mat. nette
mynames[mynames==AA02]-Partecipazioni e crediti fin.
mynames[mynames==AA02B]-Attivita' fin. a breve
mynames[mynames==AA02L]-Immobilizzazioni finan.
mynames[mynames==AA03]-Magazzino
mynames[mynames==AA04]-Crediti commerciali
mynames[mynames==AA05]-Liquidita'
mynames[mynames==AA06]-Altre attivita'
mynames[mynames==AA07]-Tot attivita'
mynames[mynames==AL01]-Capitale netto
mynames[mynames==AL02]-Fondo tfr
mynames[mynames==AL03]-Altri fondi
mynames[mynames==AL04]-Debiti commerciali
mynames[mynames==AL04A]-Anticipi di clienti
mynames[mynames==AL04B]-Debiti vs fornitori
mynames[mynames==AL05]-Debiti fin. tot.
mynames[mynames==AL05B]-Debiti fin. a breve
mynames[mynames==AL05L]-Debiti fin. a medio/lungo
mynames[mynames==AL99]-Altre passivita'
mynames[mynames==AL06]-Tot passivita'
mynames[mynames==EC01]-Ricavi netti
mynames[mynames==EC02]-Produzione int. capitalizzate
mynames[mynames==EC03]=Variazione scorte prod finiti
mynames[mynames==EC04]-Acquisti
mynames[mynames==EC05]-Variazioni scorte mat. prime
mynames[mynames==EC06]-Costi per servizi e god. di beni terzi
mynames[mynames==EC07]-Costo del lavoro tot
mynames[mynames==EC08]-Ammortamenti e accantonamenti
mynames[mynames==EC08A]-Ammortamenti
mynames[mynames==EC08B]-Accantonamenti e utilizzi di riserve
mynames[mynames==EC09]-Oneri fin.
mynames[mynames==EC10]-Proventi fin.
mynames[mynames==EC11]-Ricavi diversi netti
mynames[mynames==EC11A]-Altri ricavi netti ord
mynames[mynames==EC11C]-Contributi in conto esercizio
mynames[mynames==EC12]-Proventi straord. netti
mynames[mynames==EC13]-Imposte
mynames[mynames==EC14]-Utile netto rettificato
mynames[mynames==EC15]-Rettifiche
mynames[mynames==EC16]-Utile dell'esercizio
names(report)-mynames
report=split(report,report$CFISCALE)
report1=lapply(report,function(x){
t(x)})
@
echo=F,results=tex=
report2=lapply(report1, function(x) {
print(xtable(x,na.string=-))})
@
\end{document}


Even if I put the code referring to the title, in my pdf document I don't get 
it and I don't know why. Secondly I get the following error message:(\end 
occurred when \ifnum on line 2150 was incomplete)
the results of this error is that I loose a lot of tables. Instead of having 
500 tables I have just 250 tables. Another problem is that the format is not 
what I would like to get, my tables appear at the center of the page and I 
would like to put them at the left (the result is that my table are cut), and I 
don't know how to do that, I've tried to put in the xtable function the option 
table.placement=H but it seems that it doesn't work.An example of what I get 
by using the split function and then the xtable function is:
amp; 49 amp; 48 amp; 47 \\ 
  \hline
CFISCALE amp; 5060157 amp; 5060157 amp; 5060157 \\ 
  RAGSOCB amp; GIUSEPPE TARENZI S.R.L. amp; GIUSEPPE TARENZI S.R.L. amp; 
GIUSEPPE TARENZI S.R.L. \\ 
  ANNO amp; 2005 amp; 2006 amp; 2007 \\ 
  Ricavi netti amp;  77 amp;  98 amp; 124 \\ 
  Produzione int. capitalizzate amp; 0 amp; 0 amp; 0 \\ 
  Variazione scorte prod finiti amp; 2059 amp; 2105 amp; 2120 \\ 
  Acquisti amp; 1542 amp; 1564 amp; 1576 \\ 
  Costi per servizi e god. di beni terzi amp; 122 amp; 135 amp; 121 \\ 
  Costo del lavoro tot amp; 273 amp; 281 amp; 301 \\ 
  Ammortamenti e accantonamenti amp; 11 amp;  5 amp;  7 \\ 
  Ammortamenti amp; 9.9 amp; 4.5 amp; 6.3 \\ 
  Accantonamenti e utilizzi di riserve amp; 1.1 amp; 0.5 amp; 0.7 \\ 
  Oneri fin. amp; 38 amp; 42 amp; 35 \\ 
  Proventi fin. amp; 1 amp; 0 amp; 1 \\ 
  Ricavi diversi netti amp;   0 amp; -13 amp; -33 \\ 
  Altri ricavi netti ord amp;   0 amp; -13 amp; -33 \\ 
  Contributi in conto esercizio amp; 0 amp; 0 amp; 0 \\ 
  Proventi straord. netti amp; 1 amp; 0 amp; 0 \\ 
  Imposte amp; 73 amp; 78 amp; 80 \\ 
  Utile netto rettificato amp; 79 amp; 85 amp; 92 \\ 
  Utile dell'esercizio amp; 79 amp; 85 amp; 92 \\ 
  Immobilizzazioni tecniche nette amp; 269 amp; 267 amp; 274 \\ 
  Partecipazioni e crediti fin. amp; 0 amp; 3 amp; 0 \\ 
  Magazzino amp; 592 amp; 623 amp; 656 \\ 
  Crediti commerciali amp; 17.56393 amp; 28.15887 amp; 26.14891 \\ 
  Liquidita' amp;  14 amp;  amp; 177 \\ 
  Altre attivita' amp;  -892.564 amp;  amp; -1133.149 \\ 
  Tot attivita' amp; 0.0 amp; 0.005924951 amp; 0.0 \\ 
  Immobilizzazioni imm. nette amp; 53.8 amp; 53.4 amp; 54.8 \\ 
  Immobilizzazioni mat. nette amp; 215.2 amp; 213.6 amp; 219.2 \\ 
  Attivita' fin. a breve amp; 0.0 amp; 1.5 amp; 0.0 \\ 
  Immobilizzazioni finan. amp; 0.0 amp; 1.5 amp; 0.0 \\ 
  Capitale 

Re: [R] Help in the legend()

2010-07-06 Thread Dennis Murphy
Hi:

You didn't specify what you meant in your original post by
following that I had used didn't gave desired result (sic).
I noted two potential problems: (i) the n in the plot is not what
you expected it to be, and (ii) the legend didn't render in the
range of values established by the graph you produced.

Problem (i) arises because
 seq(4:13)
 [1]  1  2  3  4  5  6  7  8  9 10
The sequence you actually wanted can be obtained either
by 4:13 or seq(4, 13):
 seq(4, 13)
 [1]  4  5  6  7  8  9 10 11 12 13
 4:13
 [1]  4  5  6  7  8  9 10 11 12 13

Re the second problem, the x-axis of your graph
(when plotted correctly) runs from 4 to 13, but the x-coordinate
of your legend command starts at x = 30, so it won't appear
in the graphics window, but rather in the 'great beyond' :)

One way to get the graph and legend is to use matplot(),
because the x-axis will be the same in each curve rendered;
I observed that when using your code,
lines(*s, lty = ?)
used the x-values 1:10 rather than the 4:13 that you wanted.
An alternative is shown below.

I decided to put everything in a data frame, although this is not
strictly necessary:

dd - data.frame(n = 4:13, pg, gs, ps)
# Plot all three curves together
with(dd, matplot(n, as.matrix(dd2[, -1]), type = 'l', lty = c(1, 2, 5), col
= 1,
   xlab = 'n', ylab = 'Differences of the
variances') )
# Set up legend and insert it:
lx - expression(var(t^3) - var(t^2), var(t^2) - var(t^1),
   var(t^3) -  var(t^1))
legend(10, 0.0021, legend = lx, lty = c(1, 2, 5))

Using your code, this appears to work (pg, ps and gs same as yours):
n - 4:13
plot(n, pg, type=l,xlab=n,
  ylab=Differences of the variances, ylim=c(-0.0012,0.0023) );
lines(n, gs,lty = 2)
lines(n, ps,lty=5)
legend(10, 0.0021, legend = lx, lty = c(1, 2, 5))

HTH,
Dennis


On Mon, Jul 5, 2010 at 9:41 PM, Shant Ch sha1...@yahoo.com wrote:

 Thanks Dr. Winsemius. Here's the toy data set.

 Basically pg = var(t^(3))-var(t^(2), gs = var(t^(2))-var(t^(1))and
 ps=var(t^(3))-var(t^(1)). The revised code and the data set is as follows:

 n-seq(4:13)
 pg-c(-1.241394e-03, -9.738079e-04, -7.158755e-04, -5.343962e-04,
 -4.088778e-04, -3.202068e-04, -2.558709e-04, -2.079914e-04, -1.715435e-04,
  -1.432430e-04)
 gs-c(0.0022520038, 0.0020060234, 0.0017601434, 0.0015519810,
 0.0013810851,0.0012407732, 0.0011245410, 0.0010271681, 0.0009446642,
 0.0008740083)
 ps-c( 0.0010106098, 0.0010322155, 0.0010442678, 0.0010175848,
 0.0009722074,0.0009205665, 0.0008686700, 0.0008191768, 0.0007731207,
 0.0007307653)

 plot(n, pg, type=l,xlab=n,ylab=Differences of the
 variances,ylim=c(-0.0012,0.0023) );
 lines(gs,lty = 2)
 lines(ps,lty=5)

 legend(30, 0.0021, expression( c ( var(t^(3))-var(t^(2)),
 var(t^(2))-var(t^(1))), var(t^(3))-var(t^(1)) ) ), lty=c(1,2,5)).



 
 From: David Winsemius dwinsem...@comcast.net

 Cc: r-help@r-project.org
 Sent: Mon, July 5, 2010 9:43:19 PM
 Subject: Re: [R] Help in the legend()


 On Jul 5, 2010, at 8:06 PM, Shant Ch wrote:

  Hi R-users,
 
  I was plotting the differences of the variances of the three estimators-
 T^(1), T^(2), T^(3), ofcourse taking two at a time. I was using the
 expression() in the legend function in order to show which line correspond
 to which of the difference, but the following that I had used didn't gave
 desired result. I would be grateful, if you help me out.
 
  plot(n, pg, type=l,xlab=n,ylab=Differences of the
 variances,ylim=c(-0.0012,0.0023), xlim=c(0,60));
  lines(gs,lty = 2)
  lines(ps,lty=5)
 
  legend(30, 0.0021, expression( c ( var(t^(3))-var(t^(2)),
 var(t^(2))-var(t^(1))), var(t^(3))-var(t^(1)) ) ), lty=c(1,2,5))
 

 Have you consider offering a toy set of objects which defines t, n, and
 pg.

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

 David Winsemius, MD
 West Hartford, CT



[[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] Fast String operations in R ? Cost of String operations

2010-07-06 Thread Martin Maechler
 RB == Ralf B ralf.bie...@gmail.com
 on Mon, 5 Jul 2010 02:33:13 -0400 writes:

RB Hi experts,
RB currently developing some code that checks a large amount of Strings
RB for the existence of sub-strings and pattern (detecting sub-strings
RB within URLs). I wonder if there is information about how well
RB particular String operations work in R together with comparisons. Are
RB there  recommendations (based on such information) regarding what
RB operations should be used and what should be avoided? Are there
RB libraries and functions that provide optimized String operations for
RB such needs or is R simply not the right choice for that?

Why not?
R itself is nowadays used to build R packages, the need for Perl
basically gone, so why do you assume that R's string operations
are not fast enough for your task at hand?

Please read and adhere to the posting guide,
see the footer of *every* R-help message, cited below...
((and then work, and may be then ask a more specific question))

Martin Maechler, ETH Zurich

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Function for gruping similar variables?

2010-07-06 Thread Timo W
Hi,

I have a matrix of results of multiple 2x2 chi^2 tests, non-
significant tests are marked as TRUE. Is there a function for grouping
those variables in a similar way LSD.test from agricolae library does?
I reviewed LSD.test's source but it's not helpful for me.

This is my matrix:

   [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]
1  TRUE FALSE  TRUE  TRUE  TRUE FALSE  TRUE
2 FALSE  TRUE FALSE FALSE FALSE FALSE FALSE
3  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE
4  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE
5  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE
6 FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE
7  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE

And I'd like to get similar output to this (LSD.test):

Groups, Treatments and means
aoo  36.9
aff  36.3
 b   cc  24.4
  c  fc  12.86667

For example, there's a group of variables 4, 5, 6,  7 that don't
significantly differ from each other. Is there any _simple_ way to
mark all of them as a group?

If you have any ideas, I'd be grateful.

Thanks in advance,
Timo

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


Re: [R] Error in affypdnn package

2010-07-06 Thread Henrik Bengtsson
Hi,

repost your question to the Bioconductor mailing list instead, and
your chances for getting help will be much greater:

  http://bioconductor.org/docs/mailList.html

/Henrik

On Tue, Jul 6, 2010 at 9:06 AM, Hasan, Ahmed Ryadh - hasar001
ahmed.ha...@postgrads.unisa.edu.au wrote:
 Dear all,
 I am a PhD student working with Affymetrix HGU133atag array for analyzing the 
 Latin square experiment.

 I was trying to generate gene expression index for hgu133atag array for PDNN 
 model. While extracting the chiptype specific data structure, I got the 
 following error-


 library(affypdnn)
 Loading required package: affy
 Loading required package: Biobase
 Welcome to Bioconductor
  Vignettes contain introductory material. To view, type
  'openVignette()'. To cite Bioconductor, see
  'citation(Biobase)' and for packages 'citation(pkgname)'.
 registering new summary method 'pdnn'.
 registering new pmcorrect method 'pdnn' and 'pdnnpredict'.
 library(hgu133atagprobe)
 Loading required package: AnnotationDbi
 energy.file - system.file(exampleData, 
 pdnn-energy-parameter_hg-u133a.txt, package = affypdnn)
 params.chiptype - pdnn.params.chiptype(energy.file, probes.pack = 
 hgu133aprobe)
 Calculating chip type specific parameters, (may take some time)...
 |                    |
 |Error in object$call : $ operator not defined for this S4 class

 After this, I also tried to do it for hgu95av2 array as described in the 
 package manual-

 library(affypdnn)
 Loading required package: affy
 Loading required package: Biobase
 Welcome to Bioconductor
  Vignettes contain introductory material. To view, type
  'openVignette()'. To cite Bioconductor, see
  'citation(Biobase)' and for packages 'citation(pkgname)'.
 registering new summary method 'pdnn'.
 registering new pmcorrect method 'pdnn' and 'pdnnpredict'.
 library(hgu95av2probe)
 Loading required package: AnnotationDbi
 energy.file - system.file(exampleData, 
 pdnn-energy-parameter_hg-u95av2.txt, package = affypdnn)
 params.chiptype - pdnn.params.chiptype(energy.file, probes.pack = 
 hgu95av2probe)
 Calculating chip type specific parameters, (may take some time)...
 |                    |
 |Error in object$call : $ operator not defined for this S4 class

 Please help me to get rid of this error(Error in object$call : $ operator not 
 defined for this S4 class).

 Best regards,

 Ahmed Ryadh Hasan

 PhD Student,
 Knowledge-Based Intelligent Engineering Systems centre (KES)
 Mawson Lakes Campus, University of South Australia,
 Mawson Lakes, SA 5095

 Phone: +61 8 8302 8332 (Office)
       +61 8 8359 5747 (Home)
       +61 4 3006 4390 (Cell)

        [[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] xyplot: filtering out empty plots

2010-07-06 Thread Coen van Hasselt
Hello,

I would like to know how I can filter out empty plots in xyplot, when
stratifying on some variables.

Example:

I have a dataset in which I plot CONC ~ TIME, stratified for patient
ID(1,2,..,100), FORM(1,2) and BOOST (1,2).
Some patients (ID's) do not have values for all stratification
conditions. I.e. one patient may have values for FORM=1 and BOOST=1,2,
while others may have data on alle combinations.

xyplot(CONC ~ TIME | ID + FORM + BOOST, dat=data)

If I do this I also get empty plots for all the combinations for which
I do not have data for. Is it possible for xyplot NOT to make plots
for combinations of ID, FORM and BOOST that do not contain data ?

Thanks in advance!

Coen

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Assign Formulas to Arrays or Matrices?

2010-07-06 Thread Dennis Murphy
Hi:

See inline...

On Tue, Jul 6, 2010 at 12:13 AM, McLovin dave_dec...@hotmail.com wrote:


 Hi,

 I am very new to R.  I am hoping to create formulas and assign them to
 locations within an array (or matrix, if it will work).

 Here's a simplified example of what I'm trying to do:

 form.arr - array(31,5,3)
 for (i in seq(from=1, to=31, by=1)) {
for (j in seq(from=1, to=5, by=1)) {
form.arr[i,j,] - as.formula(y~1+2)
}
 }

 which results in this error:
 Error in form.arr[i, j, ] - as.formula(y ~ 1 + 2) :
  incorrect number of subscripts


 fm - y ~ 1 + 2
 class(fm)
[1] formula

# Let's investigate the fm object a little further...
 str(fm)
Class 'formula' length 3 y ~ 1 + 2
  ..- attr(*, .Environment)=environment: R_GlobalEnv
 fm[1]
`~`()
 fm[2]
y()
 fm[3]
1 + 2()

The reason I had made the 3rd dimension of the array size 3 is because
 that's the length R tells me that formula is.


True, but...

 sapply(fm, class)
[1] name name call

so the classes of the components of the formula are not the same.
An array requires that each element of the array have the same class.

When I had tried to do this using a matrix, using this code:

 form.mat - matrix(31,5,3)
 for (i in seq(from=1, to=31, by=1)) {
for (j in seq(from=1, to=5, by=1)) {
form.mat[i,j] = as.formula(y~1+2)
}
 }

 I was told:

 Error in form.mat[i, j] = as.formula(y ~ 1 + 2) :
  number of items to replace is not a multiple of replacement length


You're trying to assign a three-element object to a single subscript
of a matrix - R is (correctly) telling you that doesn't compute.


 My question is: is it possible to assign formulas within a matrix or
 array?

If so, how?  thanks@real.com


I think you might have better luck trying to assign formulas to list
components,
since the classes of the components of a formula are not the same. For
this reason, you can't cram them into arrays or matrices for the reason
given above.

HTH,
Dennis


 --
 View this message in context:
 http://r.789695.n4.nabble.com/Assign-Formulas-to-Arrays-or-Matrices-tp2279136p2279136.html
 Sent from the R help mailing list archive at Nabble.com.

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


[[alternative HTML version deleted]]

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


Re: [R] Sweave function

2010-07-06 Thread Eik Vettorazzi
Hi,
I don't know why, but your resulting table is not latex but html. (e.g.
amp; is not recognized in latex).

NA.string is an argument for print.xtable, not for xtable itself and is
case sensitive.

So

echo=F,results=tex=
report2=lapply(report1, function(x) {
print(xtable(x),NA.string=-)})
@

should work.
hth.


Am 06.07.2010 11:13, schrieb n.via...@libero.it:
 Dear list,
 I'm trying to generate a latex Document in which there are a lot of tables. 
 I'm using the Sweave function in the package utils, but I'm having a lot of 
 problems with the format. This is my code:
 \documentclass[a4paper]{amsbook}
 \title{Schema di bilancio}
 \begin{document}
 \maketitle
 echo=F,results=hide=
 report=Bilanci
 mynames-names(report)
 mynames[mynames==AA01]-Immobilizzazioni tecniche nette
 mynames[mynames==AA01I]-Immobilizzazioni imm. nette
 mynames[mynames==AA01M]-Immobilizzazioni mat. nette
 mynames[mynames==AA02]-Partecipazioni e crediti fin.
 mynames[mynames==AA02B]-Attivita' fin. a breve
 mynames[mynames==AA02L]-Immobilizzazioni finan.
 mynames[mynames==AA03]-Magazzino
 mynames[mynames==AA04]-Crediti commerciali
 mynames[mynames==AA05]-Liquidita'
 mynames[mynames==AA06]-Altre attivita'
 mynames[mynames==AA07]-Tot attivita'
 mynames[mynames==AL01]-Capitale netto
 mynames[mynames==AL02]-Fondo tfr
 mynames[mynames==AL03]-Altri fondi
 mynames[mynames==AL04]-Debiti commerciali
 mynames[mynames==AL04A]-Anticipi di clienti
 mynames[mynames==AL04B]-Debiti vs fornitori
 mynames[mynames==AL05]-Debiti fin. tot.
 mynames[mynames==AL05B]-Debiti fin. a breve
 mynames[mynames==AL05L]-Debiti fin. a medio/lungo
 mynames[mynames==AL99]-Altre passivita'
 mynames[mynames==AL06]-Tot passivita'
 mynames[mynames==EC01]-Ricavi netti
 mynames[mynames==EC02]-Produzione int. capitalizzate
 mynames[mynames==EC03]=Variazione scorte prod finiti
 mynames[mynames==EC04]-Acquisti
 mynames[mynames==EC05]-Variazioni scorte mat. prime
 mynames[mynames==EC06]-Costi per servizi e god. di beni terzi
 mynames[mynames==EC07]-Costo del lavoro tot
 mynames[mynames==EC08]-Ammortamenti e accantonamenti
 mynames[mynames==EC08A]-Ammortamenti
 mynames[mynames==EC08B]-Accantonamenti e utilizzi di riserve
 mynames[mynames==EC09]-Oneri fin.
 mynames[mynames==EC10]-Proventi fin.
 mynames[mynames==EC11]-Ricavi diversi netti
 mynames[mynames==EC11A]-Altri ricavi netti ord
 mynames[mynames==EC11C]-Contributi in conto esercizio
 mynames[mynames==EC12]-Proventi straord. netti
 mynames[mynames==EC13]-Imposte
 mynames[mynames==EC14]-Utile netto rettificato
 mynames[mynames==EC15]-Rettifiche
 mynames[mynames==EC16]-Utile dell'esercizio
 names(report)-mynames
 report=split(report,report$CFISCALE)
 report1=lapply(report,function(x){
 t(x)})
 @
 echo=F,results=tex=
 report2=lapply(report1, function(x) {
 print(xtable(x,na.string=-))})
 @
 \end{document}


 Even if I put the code referring to the title, in my pdf document I don't get 
 it and I don't know why. Secondly I get the following error message:(\end 
 occurred when \ifnum on line 2150 was incomplete)
 the results of this error is that I loose a lot of tables. Instead of having 
 500 tables I have just 250 tables. Another problem is that the format is not 
 what I would like to get, my tables appear at the center of the page and I 
 would like to put them at the left (the result is that my table are cut), and 
 I don't know how to do that, I've tried to put in the xtable function the 
 option table.placement=H but it seems that it doesn't work.An example of 
 what I get by using the split function and then the xtable function is:
 amp; 49 amp; 48 amp; 47 \\ 
   \hline
 CFISCALE amp; 5060157 amp; 5060157 amp; 5060157 \\ 
   RAGSOCB amp; GIUSEPPE TARENZI S.R.L. amp; GIUSEPPE TARENZI S.R.L. amp; 
 GIUSEPPE TARENZI S.R.L. \\ 
   ANNO amp; 2005 amp; 2006 amp; 2007 \\ 
   Ricavi netti amp;  77 amp;  98 amp; 124 \\ 
   Produzione int. capitalizzate amp; 0 amp; 0 amp; 0 \\ 
   Variazione scorte prod finiti amp; 2059 amp; 2105 amp; 2120 \\ 
   Acquisti amp; 1542 amp; 1564 amp; 1576 \\ 
   Costi per servizi e god. di beni terzi amp; 122 amp; 135 amp; 121 \\ 
   Costo del lavoro tot amp; 273 amp; 281 amp; 301 \\ 
   Ammortamenti e accantonamenti amp; 11 amp;  5 amp;  7 \\ 
   Ammortamenti amp; 9.9 amp; 4.5 amp; 6.3 \\ 
   Accantonamenti e utilizzi di riserve amp; 1.1 amp; 0.5 amp; 0.7 \\ 
   Oneri fin. amp; 38 amp; 42 amp; 35 \\ 
   Proventi fin. amp; 1 amp; 0 amp; 1 \\ 
   Ricavi diversi netti amp;   0 amp; -13 amp; -33 \\ 
   Altri ricavi netti ord amp;   0 amp; -13 amp; -33 \\ 
   Contributi in conto esercizio amp; 0 amp; 0 amp; 0 \\ 
   Proventi straord. netti amp; 1 amp; 0 amp; 0 \\ 
   Imposte amp; 73 amp; 78 amp; 80 \\ 
   Utile netto rettificato amp; 79 amp; 85 amp; 92 \\ 
   Utile dell'esercizio amp; 79 amp; 85 amp; 92 \\ 
   Immobilizzazioni tecniche nette amp; 269 amp; 267 amp; 274 \\ 
   Partecipazioni e crediti fin. amp; 0 amp; 3 amp; 0 \\ 
   Magazzino amp; 592 amp; 623 amp; 656 \\ 
   Crediti commerciali amp; 

[R] Could not find createData function

2010-07-06 Thread shabnam k
Hi,

  I am using *Maanova* package to do anova. I have created *datafile* with
probeID as the first column, which is a tab limited text file and also
created *designfile*. I have created *readma object* which is named as abf1.
From that readma object, i have to create data object by using
*createData*function and also i hav to create model object by using
*makemodel* function, then i have to fit the model of anova.But, the problem
is it could not find createData function. Am pasting the commands which i
used below.Please give me the solution to my problem, as am unabl;e to
proceed further.


R version 2.11.1 (2010-05-31)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0


R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

 library(affy)
Loading required package: Biobase

Welcome to Bioconductor


  'openVignette()'. To cite Bioconductor, see
  'citation(Biobase)' and for packages 'citation(pkgname)'.

 library(maanova)

Attaching package: 'maanova'

The following object(s) are masked from 'package:base':

norm

 abf1.raw - read.madata(gcrmadata.expr.
dat, designfile=design.dat,
+ probeID=1, pmt=2, spotflag=F)
Reading one color array.
 Otherwise change arrayType='twoColor' then read the data again
Warning messages:
1: In read.madata(gcrmadata.expr.dat, designfile = design.dat,  :
  Assume that the first column is probeid. If you have probeid specify it,
otherwise set 'probeid=0' then read the data again

2: In read.madata(gcrmadata.expr.dat, designfile = design.dat,  :
  Assume that intensity value is saved from the second column. Otherwise
provide 'intensity' (first column storing intensity) information, and read
the data again

 abf1.raw

Summary for this experiment

Number of dyes:  1
Number of arrays:4
Number of genes: 8799
Number of replicates:1
Transformation method:   None
Replicate collapsed: FALSE
 data(abf1)
 abf1 - *createData*(abf1, 1, log.trans=F)
*Error: *could not find function* createData




*

[[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] Assign Formulas to Arrays or Matrices?

2010-07-06 Thread Jens Oehlschlägel
for assigning formulas to arrays use an array of list

nr   form.arr[[31,5]]y ~ 1 + 2
Jens Oehlschlägel


-Ursprüngliche Nachricht-
Von: McLovin  
Gesendet: Jul 6, 2010 9:13:49 AM
An: r-help@r-project.org
Betreff: [R] Assign Formulas to Arrays or Matrices?


Hi,

I am very new to R.  I am hoping to create formulas and assign them to
locations within an array (or matrix, if it will work).

Here's a simplified example of what I'm trying to do:

form.arr  for (i in seq(from=1, to=31, by=1)) {
   for (j in seq(from=1, to=5, by=1)) {
   form.arr[i,j,]  }
}

which results in this error:
Error in form.arr[i, j, ]incorrect number of subscripts

The reason I had made the 3rd dimension of the array size 3 is because
that's the length R tells me that formula is.

When I had tried to do this using a matrix, using this code:

form.mat  for (i in seq(from=1, to=31, by=1)) {
   for (j in seq(from=1, to=5, by=1)) {
   form.mat[i,j] = as.formula(y~1+2)
   }
}

I was told:

Error in form.mat[i, j] = as.formula(y ~ 1 + 2) : 
  number of items to replace is not a multiple of replacement length

My question is: is it possible to assign formulas within a matrix or array? 
If so, how?  thanks@real.com
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Assign-Formulas-to-Arrays-or-Matrices-tp2279136p2279136.html
Sent from the R help mailing list archive at Nabble.com.

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

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


Re: [R] Assign Formulas to Arrays or Matrices?

2010-07-06 Thread Patrick Burns

On 06/07/2010 08:13, McLovin wrote:


Hi,

I am very new to R.  I am hoping to create formulas and assign them to
locations within an array (or matrix, if it will work).

Here's a simplified example of what I'm trying to do:

form.arr- array(31,5,3)
for (i in seq(from=1, to=31, by=1)) {
for (j in seq(from=1, to=5, by=1)) {
form.arr[i,j,]- as.formula(y~1+2)
}
}

which results in this error:
Error in form.arr[i, j, ]- as.formula(y ~ 1 + 2) :
   incorrect number of subscripts

The reason I had made the 3rd dimension of the array size 3 is because
that's the length R tells me that formula is.

When I had tried to do this using a matrix, using this code:

form.mat- matrix(31,5,3)
for (i in seq(from=1, to=31, by=1)) {
for (j in seq(from=1, to=5, by=1)) {
form.mat[i,j] = as.formula(y~1+2)
}
}

I was told:

Error in form.mat[i, j] = as.formula(y ~ 1 + 2) :
   number of items to replace is not a multiple of replacement length

My question is: is it possible to assign formulas within a matrix or array?
If so, how?  thanks@real.com


You can create a matrix of mode list to do that.
An example is on page 39 of 'The R Inferno'.

--
Patrick Burns
pbu...@pburns.seanet.com
http://www.burns-stat.com
(home of 'Some hints for the R beginner'
and 'The R Inferno')

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


[R] Error in createData function

2010-07-06 Thread shabnam k
Hi,

  I am using *Maanova* package to do anova. I have created *datafile* with
probeID as the first column, which is a tab limited text file and also
created *designfile*. I have created *readma object* which is named as abf1.
From that readma object, i have to create data object by using
*createData*function and also i hav to create model object by using
*makemodel* function, then i have to fit the model of anova.But, the problem
is it could not find createData function. Am pasting the commands which i
used below.Please give me the solution to my problem, as am unabl;e to
proceed further.


R version 2.11.1 (2010-05-31)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0


R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

 library(affy)
Loading required package: Biobase

Welcome to Bioconductor


  'openVignette()'. To cite Bioconductor, see
  'citation(Biobase)' and for packages 'citation(pkgname)'.

 library(maanova)

Attaching package: 'maanova'

The following object(s) are masked from 'package:base':

norm

 abf1.raw - read.madata(gcrmadata.expr.
dat, designfile=design.dat,
+ probeID=1, pmt=2, spotflag=F)
Reading one color array.
 Otherwise change arrayType='twoColor' then read the data again
Warning messages:
1: In read.madata(gcrmadata.expr.dat, designfile = design.dat,  :
  Assume that the first column is probeid. If you have probeid specify it,
otherwise set 'probeid=0' then read the data again

2: In read.madata(gcrmadata.expr.dat, designfile = design.dat,  :
  Assume that intensity value is saved from the second column. Otherwise
provide 'intensity' (first column storing intensity) information, and read
the data again

 abf1.raw

Summary for this experiment

Number of dyes:  1
Number of arrays:4
Number of genes: 8799
Number of replicates:1
Transformation method:   None
Replicate collapsed: FALSE
 data(abf1)
 abf1 - *createData*(abf1, 1, log.trans=F)
*Error: *could not find function* createData




*





Shabnam K.

[[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 simplify a table

2010-07-06 Thread Simone Gabbriellini
Hello List,

as a result of a table(), I have this:

 0  1   3  4  5  6  7  9 11 12 13 14 16 17 18 19 20 21 22 24 27 28
27  2  2  2  2  2  1  2  1   3   2   1   2  1   1   3  1   1   1   2  1   1

I would like to simplify it by grouping like:

0   0-5 6-10 11-15 16-20 21-25 26-30
25  84   78   4   2

is it possible? could someone point me to the right function to do this?

best regards,
Simone

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

2010-07-06 Thread Timo W
see:

?cut

On 6 Lip, 12:49, Simone Gabbriellini simone.gabbriell...@gmail.com
wrote:
 Hello List,

 as a result of a table(), I have this:

  0  1   3  4  5  6  7  9 11 12 13 14 16 17 18 19 20 21 22 24 27 28
 27  2  2  2  2  2  1  2  1   3   2   1   2  1   1   3  1   1   1   2  1   1

 I would like to simplify it by grouping like:

 0   0-5 6-10 11-15 16-20 21-25 26-30
 25  8    4       7        8       4       2

 is it possible? could someone point me to the right function to do this?

 best regards,
 Simone

 __
 r-h...@r-project.org mailing listhttps://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guidehttp://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] locfit

2010-07-06 Thread Arnab Maity
Hi,

Can you provide me an example in R to estimate the density using locfit package 
with the help of multi dimensional explanatory variables and one dimensional 
dependent variable?
Thank you.

Arnab Kumar Maity
Research Assistant
Indian School of Business
Hyderabad, India

DISCLAIMER:\ This e-mail (including any attachments) is ...{{dropped:12}}

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

2010-07-06 Thread Martin Morgan
Hi --

On 07/06/2010 03:21 AM, shabnam k wrote:
 Hi,
 
   I am using *Maanova* package to do anova. I have created *datafile* with
 probeID as the first column, which is a tab limited text file and also
 created *designfile*. I have created *readma object* which is named as abf1.
From that readma object, i have to create data object by using
 *createData*function and also i hav to create model object by using
 *makemodel* function, then i have to fit the model of anova.But, the problem
 is it could not find createData function. Am pasting the commands which i
 used below.Please give me the solution to my problem, as am unabl;e to
 proceed further.

maanova is a Bioconductor package so please ask on the Bioconductor
mailing list

  http://bioconductor.org/docs/mailList.html

The vignette for this package

  browseVignettes('maanova')

does not mention createData; what instructions are you following?

Martin

 
 
 R version 2.11.1 (2010-05-31)
 Copyright (C) 2010 The R Foundation for Statistical Computing
 ISBN 3-900051-07-0
 
 
 R is free software and comes with ABSOLUTELY NO WARRANTY.
 You are welcome to redistribute it under certain conditions.
 Type 'license()' or 'licence()' for distribution details.
 
   Natural language support but running in an English locale
 
 R is a collaborative project with many contributors.
 Type 'contributors()' for more information and
 'citation()' on how to cite R or R packages in publications.
 
 Type 'demo()' for some demos, 'help()' for on-line help, or
 'help.start()' for an HTML browser interface to help.
 Type 'q()' to quit R.
 
 library(affy)
 Loading required package: Biobase
 
 Welcome to Bioconductor
 
 
   'openVignette()'. To cite Bioconductor, see
   'citation(Biobase)' and for packages 'citation(pkgname)'.
 
 library(maanova)
 
 Attaching package: 'maanova'
 
 The following object(s) are masked from 'package:base':
 
 norm
 
 abf1.raw - read.madata(gcrmadata.expr.
 dat, designfile=design.dat,
 + probeID=1, pmt=2, spotflag=F)
 Reading one color array.
  Otherwise change arrayType='twoColor' then read the data again
 Warning messages:
 1: In read.madata(gcrmadata.expr.dat, designfile = design.dat,  :
   Assume that the first column is probeid. If you have probeid specify it,
 otherwise set 'probeid=0' then read the data again
 
 2: In read.madata(gcrmadata.expr.dat, designfile = design.dat,  :
   Assume that intensity value is saved from the second column. Otherwise
 provide 'intensity' (first column storing intensity) information, and read
 the data again
 
 abf1.raw
 
 Summary for this experiment
 
 Number of dyes:  1
 Number of arrays:4
 Number of genes: 8799
 Number of replicates:1
 Transformation method:   None
 Replicate collapsed: FALSE
 data(abf1)
 abf1 - *createData*(abf1, 1, log.trans=F)
 *Error: *could not find function* createData
 
 
 
 
 *
 
   [[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.


-- 
Martin Morgan
Computational Biology / Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N.
PO Box 19024 Seattle, WA 98109

Location: Arnold Building M1 B861
Phone: (206) 667-2793

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


Re: [R] multiple time series plot with dual 'y' axes

2010-07-06 Thread Jim Lemon

On 07/06/2010 06:45 PM, Jorge Hernandez wrote:

Hello.


I would like to know how to generate dual 'y' axes on a multiple time
series plot.


I am using  ts.plot()  to get the multiple time series plot, but I would
like a second
vertical axis on the right to include another time series on a different
scale.


Hi Jorge,
The code in the twoord.plot function in the plotrix package might help you.

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.


[R] Pseudo F statistics with index.G1

2010-07-06 Thread Kennedy

Hello,

I have done some clustering with Agnes and want to calculate the pseudo F
statistics using index.G1. It works for a low number of clusters but when i
increase the number of clusters i eventually get the following message:

Error in apply(x[cl == i, ], 2, mean) : 
  dim(X) must have a positive length

The following code produces an example comparable to the data i am
clustering:

library(cluster)
library(ade4)
library(R2HTML)
library(e1071)
library(class)
library(rgl)
library(MASS)
library(clusterSim)

# Create a symmetric matrix with ones on the diagonal
mat - diag(1,27,27)
f - runif(sum(26:1),0,1)
mat[lower.tri(mat)] - f
mat - t(mat)
mat[lower.tri(mat)] - f

# Cluster with Agnes
A - agnes(mat,diss=T,method=average)
C - cutree(A,k=7)   # Value of k = the number of clusters
F - index.G1(mat,C)

The code above works for k=2:6 but then the error message appears. 


Sincerely

Henrik
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Pseudo-F-statistics-with-index-G1-tp2279216p2279216.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] information reduction-database management question

2010-07-06 Thread Paul Chatfield

If you redefine your NAs as below to be detected as some arbitrary large
number, then the code should work through.  Any 5's left in your dataset can
be replaced just as easily by NAs again.  Not elegant, but effective.

site - c(s1, s1, s1, s2,s2, s2)
pref - c(1, 2, 3, 1, 2, 3)
R1 - c(NA, NA, 1, NA,NA,NA)
R2 - c(NA, 0, 1, 1, NA, 1)
R3 - c(NA, 1, 1, NA, 1, 1)
R4 - c(0, NA, 0, 1, NA, 0)
R5 - c(NA, 0, 1, NA, 1, 1)
datum - data.frame(site, pref, R1, R2, R3, R4, R5)

## For 1 column;
datum$R1[is.na(datum$R1)==T]-5
tapply(datum$R1, datum$site, min, na.rm=T)

## Can loop this over all columns;

new-matrix(0,5,2)
for (i in 3:7)
{datum[,i][is.na(datum[,i])==T]-5
new[i-2,]-tapply(datum[,i], datum$site, min, na.rm=T)}

-- 
View this message in context: 
http://r.789695.n4.nabble.com/information-reduction-database-management-question-tp2278863p2279385.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Linux-Windows problem

2010-07-06 Thread Bos, Roger
Uwe,

I suspect I might be having a similar problem as the R code I generate
in Windows editor (Tinn-R) doesn't always open properly in my Linux
editor (Rkward), but I don't know how to change the default encoding in
either one.  

In R on both machines, getOption(encoding) returns native.enc, so
the problem is not with R, but rather at the OS or editor level.  If you
or anyone could help me out that would be appreciated.

Thanks,

Roger 

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
On Behalf Of Uwe Ligges
Sent: Monday, July 05, 2010 10:04 AM
To: Ildiko Varga
Cc: r-help@r-project.org
Subject: Re: [R] Linux-Windows problem



On 05.07.2010 14:31, Ildiko Varga wrote:
 Dear All,

 I faced the following problem. With the same data.frame the results 
 are different under Linux and Windows.
 Could you help on this topic?

I guess you read in the data differently since you have different
default encodings on both platforms (e.g. latin1 vs. UTF-8) and you data
is probably not plain ASCII.

Best,
Uwe Ligges

 Thanks in advance,
 Ildiko

 Linux:
 d = read.csv(CRP.csv)
 d$drugCode = as.numeric(d$drug)
 cor(d, use=pairwise.complete.obs)
 PATIENT BL.CRP   X24HR.CRP   X48HR.CRP drug   drugCode
 PATIENTNA NA  NA  NA   NA NA
 BL.CRP NA  1.000  0.84324880 -0.05699590   NA -0.3367147
 X24HR.CRP  NA  0.8432488  1. -0.06162383   NA -0.3557316
 X48HR.CRP  NA -0.0569959 -0.06162383  1.   NA  0.1553356
 drug   NA NA  NA  NA   NA NA
 drugCode   NA -0.3367147 -0.35573159  0.15533562   NA  1.000
 Warning message:
 In cor(d, use = pairwise.complete.obs) : NAs introduced by coercion
 str(d)
 'data.frame':   41 obs. of  6 variables:
$ PATIENT  : Factor w/ 41 levels RV13,RV14,..: 2 3 4 6 7 12 13 
 14
 15 17 ...
$ BL.CRP   : num  7.3 31.2 4.2 6.7 1.6 7.7 5.3 38.9 1 7.3 ...
$ X24HR.CRP: num  6.1 24.9 11.1 4.9 1 5 3.7 18 1 7.3 ...
$ X48HR.CRP: num  121.5 40 28.4 34.5 33.3 ...
$ drug : Factor w/ 2 levels active,placebo: 1 1 1 1 1 1 1 1
1
 1 ...
$ drugCode : num  1 1 1 1 1 1 1 1 1 1 ...

   Windows:
 d = read.csv(CRP.csv)
 d$drugCode = as.numeric(d$drug)
 cor(d, use=pairwise.complete.obs) Error in cor(d, use = 
 pairwise.complete.obs) : 'x' must be numeric
 str(d)
 'data.frame':   41 obs. of  6 variables:
$ PATIENT  : Factor w/ 41 levels RV13,RV14,..: 2 3 4 6 7 12 13 
 14
 15 17 ...
$ BL.CRP   : num  7.3 31.2 4.2 6.7 1.6 7.7 5.3 38.9 1 7.3 ...
$ X24HR.CRP: num  6.1 24.9 11.1 4.9 1 5 3.7 18 1 7.3 ...
$ X48HR.CRP: num  121.5 40 28.4 34.5 33.3 ...
$ drug : Factor w/ 2 levels active,placebo: 1 1 1 1 1 1 1 1
1
 1 ...
$ drugCode : num  1 1 1 1 1 1 1 1 1 1 ...

   [[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.
***

This message is for the named person's use only. It may\...{{dropped:20}}

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

2010-07-06 Thread Peter Ehlers

On 2010-07-06 4:21, shabnam k wrote:

Hi,

   I am using *Maanova* package to do anova. I have created *datafile* with
probeID as the first column, which is a tab limited text file and also
created *designfile*. I have created *readma object* which is named as abf1.

From that readma object, i have to create data object by using

*createData*function and also i hav to create model object by using
*makemodel* function, then i have to fit the model of anova.But, the problem
is it could not find createData function. Am pasting the commands which i
used below.Please give me the solution to my problem, as am unabl;e to
proceed further.


R version 2.11.1 (2010-05-31)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0


R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

   Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.


library(affy)

Loading required package: Biobase

Welcome to Bioconductor


   'openVignette()'. To cite Bioconductor, see
   'citation(Biobase)' and for packages 'citation(pkgname)'.


library(maanova)


Attaching package: 'maanova'

The following object(s) are masked from 'package:base':

 norm


abf1.raw- read.madata(gcrmadata.expr.

dat, designfile=design.dat,
+ probeID=1, pmt=2, spotflag=F)
Reading one color array.
  Otherwise change arrayType='twoColor' then read the data again
Warning messages:
1: In read.madata(gcrmadata.expr.dat, designfile = design.dat,  :
   Assume that the first column is probeid. If you have probeid specify it,
otherwise set 'probeid=0' then read the data again

2: In read.madata(gcrmadata.expr.dat, designfile = design.dat,  :
   Assume that intensity value is saved from the second column. Otherwise
provide 'intensity' (first column storing intensity) information, and read
the data again


abf1.raw


 Summary for this experiment

Number of dyes:  1
Number of arrays:4
Number of genes: 8799
Number of replicates:1
Transformation method:   None
Replicate collapsed: FALSE

data(abf1)
abf1- *createData*(abf1, 1, log.trans=F)

*Error: *could not find function* createData



I think that you're using code from an ancient demo file with
a more recent version of maanova. There used to be a
createData() function and read.madata used to have
arguments 'probeID' and 'pmt', but those are now 'probeid'
and, I think, 'intensity'. Try the demo in your (*unstated*)
version of maanova.

  -Peter Ehlers

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

2010-07-06 Thread David Winsemius


On Jul 6, 2010, at 7:08 AM, Arnab Maity wrote:


Hi,

Can you provide me an example in R to estimate the density using  
locfit package with the help of multi dimensional explanatory  
variables and one dimensional dependent variable?


When I came up with a solution I posted it:

https://stat.ethz.ch/pipermail/r-help/2008-January/151755.html


Thank you.

Arnab Kumar Maity
Research Assistant
Indian School of Business
Hyderabad, India

DISCLAIMER:\ This e-mail (including any attachments) is ...{{dropped: 
12}}

--

David Winsemius, MD
West Hartford, CT

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


Re: [R] adding a row of names to data.frame

2010-07-06 Thread jamaas

Thanks Pete,

It works fine.  I should have explained, already had the data.frame with
data in it, just wanted to add a string vector into the first column of the
data.frame, could not figure out how to do that.

Jim
-- 
View this message in context: 
http://r.789695.n4.nabble.com/adding-a-row-of-names-to-data-frame-tp2278278p2279421.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Conditional Splitting a Vectors into Two Vectors

2010-07-06 Thread Christos Argyropoulos

One possible way is the following:

x -c(0.49534,0.80796,0.93970,0.8)
count -c(0,33,0,4)

 x[count==0]
[1] 0.49534 0.93970

 x[count0]
[1] 0.80796 0.8

Christos

 Date: Tue, 6 Jul 2010 15:39:08 +0900
 From: gunda...@gmail.com
 To: r-h...@stat.math.ethz.ch
 Subject: [R] Conditional Splitting a Vectors into Two Vectors
 
 Suppose I have two vectors of same dimensions:
 
 x -c(0.49534,0.80796,0.93970,0.8)
 count  -c(0,33,0,4)
 
 How can I group the vectors 'x' into two vectors:
 
  1. Vector `grzero` that contain value in x with `count` value greater
 than 0 and
  2. Vector `eqzero` with value in x with `count` value equal to zero.
 
 Yielding
 
  print(grzero)
  [1] 0.80796 0.8
  print(eqzero)
  [1] 0.49534 0.93970
 
 Regards,
 G.V.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
  
_
Hotmail: Trusted email with powerful SPAM protection.

[[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] adding a row of names to data.frame

2010-07-06 Thread David Winsemius


On Jul 6, 2010, at 8:00 AM, jamaas wrote:



Thanks Pete,

It works fine.  I should have explained, already had the data.frame  
with
data in it, just wanted to add a string vector into the first column  
of the

data.frame, could not figure out how to do that.


?cbind   # generic, has a data.frame method




Jim
--
View this message in context: 
http://r.789695.n4.nabble.com/adding-a-row-of-names-to-data-frame-tp2278278p2279421.html
Sent from the R help mailing list archive at Nabble.com.

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


David Winsemius, MD
West Hartford, CT

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


Re: [R] A question about conducting crossed random effects in R

2010-07-06 Thread Paul Chatfield

Sounds distinctly like an assignment you've been set for which we wouldn't
help.  All I'll say is crossed random effects can be dealt with effectively
in lmer.  See that for more.
-- 
View this message in context: 
http://r.789695.n4.nabble.com/A-question-about-conducting-crossed-random-effects-in-R-tp2278443p2279508.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] How do I plat timestamed (in seconds) data in R?

2010-07-06 Thread LosemindL

Hi all,

I have an Nx2 array, where the first column contains the timestamps and the
second column contains the corresponding data.

second  | ts
| --
14:25:00| 18
14:25:02| 14
14:25:04| 11
14:25:06| 4
14:25:08| 24
14:25:10| 13
14:25:12| 12
14:25:14| 6
14:25:16| 21
14:25:18| 37
14:25:20| 21
14:25:22| 9

Since I don't know how to do this in R yet, please allow me to show data in
Matlab:

In Matlab, the first column becomes:


   0.60069496185
   0.600752314785495
   0.600810185191222
   0.600868055596948
   0.600925925886258
   0.600983796291985
   0.60104197711
   0.601099536987022
   0.601157407392748

How do I plot such data with X-axis showing the timestamps (seconds) and the
Y-axis showing the data?

Thanks a lot!
-- 
View this message in context: 
http://r.789695.n4.nabble.com/How-do-I-plat-timestamed-in-seconds-data-in-R-tp2279519p2279519.html
Sent from the R help mailing list archive at Nabble.com.

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


[R] Help With ANOVA

2010-07-06 Thread Amit Patel
Hi I needed some help with ANOVA

I have a problem with My ANOVA
analysis. I have a dataset with a known ANOVA p-value, however I can
not seem to re-create it in R.

I have created a list (zzzanova) which contains
1)Intensity Values
2)Group Number (6 Different Groups)
3)Sample Number (54 different samples)
this is created by the script in Appendix 1

I then conduct ANOVA with the command
 zzz.aov - aov(Intensity ~ Group, data = zzzanova)

I get a p-value of
Pr(F)1 
0.9483218 

The
expected p-value is 0.00490 so I feel I maybe using ANOVA incorrectly
or have put in a wrong formula. I am trying to do an ANOVA analysis
across all 6 Groups. Is there something wrong with my formula. But I think I
have made a mistake in the formula rather than anything else.




APPENDIX 1

datalist - c(-4.60517, -4.60517, -4.60517, -4.60517, -4.60517, -4.60517, 
-4.60517, 3.003749, -4.60517, 
2.045314, 2.482557, -4.60517, -4.60517, -4.60517, -4.60517, 1.592743, 
-4.60517,
-4.60517, 0.91328, -4.60517, -4.60517, 1.827744, 2.457795, 0.355075, 
-4.60517, 2.39127,
2.016987, 2.319903, 1.146683, -4.60517, -4.60517, -4.60517, 1.846162, 
-4.60517, 2.121427, 1.973118,
-4.60517, 2.251568, -4.60517, 2.270724, 0.70338, 0.963816, -4.60517,  
0.023703, -4.60517,
2.043382, 1.070586, 2.768289, 1.085169, 0.959334, -0.02428, -4.60517, 
1.371895, 1.533227)

zzzanova -
structure(list(Intensity = c(t(Samp1), t(Samp2), t(Samp3), t(Samp4)), 
Group = structure(c(1,1,1,1,1,1,1,1,1,
 2,2,2,2,2,2,2,2,
 3,3,3,3,3,3,3,3,3,
 4,4,4,4,4,4,4,4,4,4,
 5,5,5,5,5,5,5,5,5,
 6,6,6,6,6,6,6,6,6), .Label = c(Group1, Group2, Group3, Group4, 
Group5, Group6), class = factor), 
Sample = structure(c( 1, 2, 3, 4, 5, 6, 7, 8, 9, 
10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 
30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 
40,41,42,43,44,45,46,47,48,49,50,51,52,53,54)
))
, .Names = c(Intensity, 
Group, Sample), row.names = 
c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
51, 52, 53, 54),class = data.frame)




__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Conditional Splitting a Vectors into Two Vectors

2010-07-06 Thread Henrique Dallazuanna
You can't try this also:

split(x, count  0)

On Tue, Jul 6, 2010 at 3:39 AM, Gundala Viswanath gunda...@gmail.comwrote:

 Suppose I have two vectors of same dimensions:

x -c(0.49534,0.80796,0.93970,0.8)
count  -c(0,33,0,4)

 How can I group the vectors 'x' into two vectors:

  1. Vector `grzero` that contain value in x with `count` value greater
 than 0 and
  2. Vector `eqzero` with value in x with `count` value equal to zero.

 Yielding

 print(grzero)
 [1] 0.80796 0.8
 print(eqzero)
 [1] 0.49534 0.93970

 Regards,
 G.V.

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




-- 
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40 S 49° 16' 22 O

[[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] Selection with changing number of columns

2010-07-06 Thread Allan Engelhardt
I'm not sure your question completely makes sense, but perhaps this will 
help:


set.seed(1)
d- data.frame(matrix(floor(runif(100, max=10)), 10))  # Example data
d[apply(d == 9, 1, any), ]  # Select rows with 9 in any column

## Or more generally:
d[ , c(1, 2, 3)] == c(2, 2, 9)   #  Or maybe d == 0:9 to select on 
all columns
apply(d[ , c(1, 2, 3)] == c(2, 2, 9), 1, any)  # any() being the general `|` 
function here
d[apply(d[ ,c(1, 2, 3)] == c(2, 2, 9), 1, any), ]  # Finally: the rows we were 
looking for


A bit over-engineered, perhaps, but gets you to use some generally 
useful functions.


Hope this helps.

Allan



On 06/07/10 09:33, Kunzler, Andreas wrote:

Dear list,

I'm looking for a way to select rows of a data.frame with changing number of 
columns (constraint) involved.

Assume a data (d) structure like


Var.1 Var.2 Var.3
9   2   1
2   9   5
1   2   1

I know the number of involved columns.

Is there a way to generate the following selection automatically (maybe for 
loop), so that it makes no difference if there are two or ten columns involved.

Selection:
d[d$Var.1==9 | d$Var.1==9 | d$Var.1==9  ,]


Does anybody know a way?

Thanks

Mit freundlichen Grüßen

Andreas Kunzler

Bundeszahnärztekammer (BZÄK)
Chausseestraße 13
10115 Berlin

Tel.: 030 40005-113
Fax:  030 40005-119

E-Mail: a.kunz...@bzaek.de

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 With ANOVA (corrected please ignore last email)

2010-07-06 Thread Amit Patel
Sorry i had a misprint in the appendix code in the last email


Hi I needed some help with ANOVA

I have a problem with My ANOVA
analysis. I have a dataset with a known ANOVA p-value, however I can
not seem to re-create it in R.

I have created a list (zzzanova) which contains
1)Intensity Values
2)Group Number (6 Different Groups)
3)Sample Number (54 different samples)
this is created by the script in Appendix 1

I then conduct ANOVA with the command
 zzz.aov - aov(Intensity ~ Group, data = zzzanova)

I get a p-value of
Pr(F)1 
0.9483218 

The
expected p-value is 0.00490 so I feel I maybe using ANOVA incorrectly
or have put in a wrong formula. I am trying to do an ANOVA analysis
across all 6 Groups. Is there something wrong with my formula. But I think I
have made a mistake in the formula rather than anything else.




APPENDIX 1

datalist - c(-4.60517, -4.60517, -4.60517, -4.60517, -4.60517, -4.60517, 
-4.60517, 3.003749, -4.60517, 
2.045314, 2.482557, -4.60517, -4.60517, -4.60517, -4.60517, 1.592743, 
-4.60517,
-4.60517, 0.91328, -4.60517, -4.60517, 1.827744, 2.457795, 0.355075, 
-4.60517, 2.39127,
2.016987, 2.319903, 1.146683, -4.60517, -4.60517, -4.60517, 1.846162, 
-4.60517, 2.121427, 1.973118,
-4.60517, 2.251568, -4.60517, 2.270724, 0.70338, 0.963816, -4.60517,  
0.023703, -4.60517,
2.043382, 1.070586, 2.768289, 1.085169, 0.959334, -0.02428, -4.60517, 
1.371895, 1.533227)

zzzanova -
structure(list(Intensity = datalist, 
Group = structure(c(1,1,1,1,1,1,1,1,1,
 2,2,2,2,2,2,2,2,
 3,3,3,3,3,3,3,3,3,
 4,4,4,4,4,4,4,4,4,4,
 5,5,5,5,5,5,5,5,5,
 6,6,6,6,6,6,6,6,6), .Label = c(Group1, Group2, Group3, Group4, 
Group5, Group6), class = factor), 
Sample = structure(c( 1, 2, 3, 4, 5, 6, 7, 8, 9, 
10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 
30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 
40,41,42,43,44,45,46,47,48,49,50,51,52,53,54)
))
, .Names = c(Intensity, 
Group, Sample), row.names = 
c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
51, 52, 53, 54),class = data.frame)




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

2010-07-06 Thread Henrik Bengtsson
Sorry for asking the obvious, but have you confirmed that you are
running the same version of R on both systems?

/H

On Tue, Jul 6, 2010 at 2:11 PM, Bos, Roger roger@rothschild.com wrote:
 Uwe,

 I suspect I might be having a similar problem as the R code I generate
 in Windows editor (Tinn-R) doesn't always open properly in my Linux
 editor (Rkward), but I don't know how to change the default encoding in
 either one.

 In R on both machines, getOption(encoding) returns native.enc, so
 the problem is not with R, but rather at the OS or editor level.  If you
 or anyone could help me out that would be appreciated.

 Thanks,

 Roger

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
 On Behalf Of Uwe Ligges
 Sent: Monday, July 05, 2010 10:04 AM
 To: Ildiko Varga
 Cc: r-help@r-project.org
 Subject: Re: [R] Linux-Windows problem



 On 05.07.2010 14:31, Ildiko Varga wrote:
 Dear All,

 I faced the following problem. With the same data.frame the results
 are different under Linux and Windows.
 Could you help on this topic?

 I guess you read in the data differently since you have different
 default encodings on both platforms (e.g. latin1 vs. UTF-8) and you data
 is probably not plain ASCII.

 Best,
 Uwe Ligges

 Thanks in advance,
 Ildiko

 Linux:
     d = read.csv(CRP.csv)
     d$drugCode = as.numeric(d$drug)
     cor(d, use=pairwise.complete.obs)
             PATIENT     BL.CRP   X24HR.CRP   X48HR.CRP drug   drugCode
 PATIENT        NA         NA          NA          NA   NA         NA
 BL.CRP         NA  1.000  0.84324880 -0.05699590   NA -0.3367147
 X24HR.CRP      NA  0.8432488  1. -0.06162383   NA -0.3557316
 X48HR.CRP      NA -0.0569959 -0.06162383  1.   NA  0.1553356
 drug           NA         NA          NA          NA   NA         NA
 drugCode       NA -0.3367147 -0.35573159  0.15533562   NA  1.000
 Warning message:
 In cor(d, use = pairwise.complete.obs) : NAs introduced by coercion
     str(d)
 'data.frame':   41 obs. of  6 variables:
    $ PATIENT  : Factor w/ 41 levels RV13,RV14,..: 2 3 4 6 7 12 13
 14
 15 17 ...
    $ BL.CRP   : num  7.3 31.2 4.2 6.7 1.6 7.7 5.3 38.9 1 7.3 ...
    $ X24HR.CRP: num  6.1 24.9 11.1 4.9 1 5 3.7 18 1 7.3 ...
    $ X48HR.CRP: num  121.5 40 28.4 34.5 33.3 ...
    $ drug     : Factor w/ 2 levels active,placebo: 1 1 1 1 1 1 1 1
 1
 1 ...
    $ drugCode : num  1 1 1 1 1 1 1 1 1 1 ...

           Windows:
     d = read.csv(CRP.csv)
     d$drugCode = as.numeric(d$drug)
     cor(d, use=pairwise.complete.obs) Error in cor(d, use =
 pairwise.complete.obs) : 'x' must be numeric
     str(d)
 'data.frame':   41 obs. of  6 variables:
    $ PATIENT  : Factor w/ 41 levels RV13,RV14,..: 2 3 4 6 7 12 13
 14
 15 17 ...
    $ BL.CRP   : num  7.3 31.2 4.2 6.7 1.6 7.7 5.3 38.9 1 7.3 ...
    $ X24HR.CRP: num  6.1 24.9 11.1 4.9 1 5 3.7 18 1 7.3 ...
    $ X48HR.CRP: num  121.5 40 28.4 34.5 33.3 ...
    $ drug     : Factor w/ 2 levels active,placebo: 1 1 1 1 1 1 1 1
 1
 1 ...
    $ drugCode : num  1 1 1 1 1 1 1 1 1 1 ...

       [[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.
 ***

 This message is for the named person's use only. It may\...{{dropped:20}}

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] Selection with changing number of columns

2010-07-06 Thread Henrique Dallazuanna
Try this:

d[colSums(d == 9)]

On Tue, Jul 6, 2010 at 5:33 AM, Kunzler, Andreas a.kunz...@bzaek.de wrote:

 Dear list,

 I'm looking for a way to select rows of a data.frame with changing number
 of columns (constraint) involved.

 Assume a data (d) structure like


 Var.1 Var.2 Var.3
 9   2   1
 2   9   5
 1   2   1

 I know the number of involved columns.

 Is there a way to generate the following selection automatically (maybe for
 loop), so that it makes no difference if there are two or ten columns
 involved.

 Selection:
 d[d$Var.1==9 | d$Var.1==9 | d$Var.1==9  ,]


 Does anybody know a way?

 Thanks

 Mit freundlichen Grüßen

 Andreas Kunzler
 
 Bundeszahnärztekammer (BZÄK)
 Chausseestraße 13
 10115 Berlin

 Tel.: 030 40005-113
 Fax:  030 40005-119

 E-Mail: a.kunz...@bzaek.de

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




-- 
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40 S 49° 16' 22 O

[[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 do I plat timestamed (in seconds) data in R?

2010-07-06 Thread David Winsemius


On Jul 6, 2010, at 8:58 AM, LosemindL wrote:



Hi all,

I have an Nx2 array, where the first column contains the timestamps  
and the

second column contains the corresponding data.

second  | ts
| --
14:25:00| 18
14:25:02| 14
14:25:04| 11
14:25:06| 4
14:25:08| 24
14:25:10| 13
14:25:12| 12
14:25:14| 6
14:25:16| 21
14:25:18| 37
14:25:20| 21
14:25:22| 9

Since I don't know how to do this in R yet, please allow me to show  
data in

Matlab:


 tz - read.table(textConnection(second  | ts
+ 14:25:00| 18
+ 14:25:02| 14
+ 14:25:04| 11
+ 14:25:06| 4
+ 14:25:08| 24
+ 14:25:10| 13
+ 14:25:12| 12
+ 14:25:14| 6
+ 14:25:16| 21
+ 14:25:18| 37
+ 14:25:20| 21
+ 14:25:22| 9), header=TRUE, sep=|)


In Matlab, the first column becomes:


  0.60069496185
  0.600752314785495
  0.600810185191222
  0.600868055596948
  0.600925925886258
  0.600983796291985
  0.60104197711
  0.601099536987022
  0.601157407392748

How do I plot such data with X-axis showing the timestamps (seconds)  
and the

Y-axis showing the data?



 plot (zoo(tz[,2], order.by=tz[,1])  )

--
David Winsemius, MD
West Hartford, CT

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


Re: [R] Pseudo F statistics with index.G1

2010-07-06 Thread Allan Engelhardt

Always use set.seed in your examples.

Running your code after set.seed(1) works fine but gives the error after 
set.seed(2).  The problem in the latter case being that there is only 
one value 7 in C and you need two or more for the index.G1 code to make 
sense.


Hope this helps a little

Allan

Example:


 set.seed(1)
 mat- diag(1,27,27)
 f- runif(sum(26:1),0,1)
 mat[lower.tri(mat)]- f
 mat- t(mat)
 mat[lower.tri(mat)]- f

 # Cluster with Agnes
 A- agnes(mat,diss=T,method=average)
 C- cutree(A,k=7)   # Value of k = the number of clusters
 F- index.G1(mat,C)
 table(C)

C
1 2 3 4 5 6 7
6 3 5 2 3 4 4

 set.seed(2)
 mat- diag(1,27,27)
 f- runif(sum(26:1),0,1)
 mat[lower.tri(mat)]- f
 mat- t(mat)
 mat[lower.tri(mat)]- f

 # Cluster with Agnes
 A- agnes(mat,diss=T,method=average)
 C- cutree(A,k=7)   # Value of k = the number of clusters
 F- index.G1(mat,C)

Error in apply(x[cl == i, ], 2, mean) :
  dim(X) must have a positive length

 table(C)

C
1 2 3 4 5 6 7
8 4 5 3 4 2 1



On 06/07/10 09:32, Kennedy wrote:

Hello,

I have done some clustering with Agnes and want to calculate the pseudo F
statistics using index.G1. It works for a low number of clusters but when i
increase the number of clusters i eventually get the following message:

Error in apply(x[cl == i, ], 2, mean) :
   dim(X) must have a positive length

The following code produces an example comparable to the data i am
clustering:

library(cluster)
library(ade4)
library(R2HTML)
library(e1071)
library(class)
library(rgl)
library(MASS)
library(clusterSim)

# Create a symmetric matrix with ones on the diagonal
mat- diag(1,27,27)
f- runif(sum(26:1),0,1)
mat[lower.tri(mat)]- f
mat- t(mat)
mat[lower.tri(mat)]- f

# Cluster with Agnes
A- agnes(mat,diss=T,method=average)
C- cutree(A,k=7)   # Value of k = the number of clusters
F- index.G1(mat,C)

The code above works for k=2:6 but then the error message appears.


Sincerely

Henrik



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 do I plat timestamed (in seconds) data in R?

2010-07-06 Thread Gabor Grothendieck
On Tue, Jul 6, 2010 at 8:58 AM, LosemindL comtech@gmail.com wrote:

 Hi all,

 I have an Nx2 array, where the first column contains the timestamps and the
 second column contains the corresponding data.

 second  | ts
 | --
 14:25:00| 18
 14:25:02| 14
 14:25:04| 11
 14:25:06| 4
 14:25:08| 24
 14:25:10| 13
 14:25:12| 12
 14:25:14| 6
 14:25:16| 21
 14:25:18| 37
 14:25:20| 21
 14:25:22| 9

 How do I plot such data with X-axis showing the timestamps (seconds) and the
 Y-axis showing the data?


Try this (and see the three vignettes, i.e. pdf documents, that come
with zoo, the zoo help files and the article on dates and times in R
News 4/1):

Lines - second  | ts
14:25:00| 18
14:25:02| 14
14:25:04| 11
14:25:06| 4
14:25:08| 24
14:25:10| 13
14:25:12| 12
14:25:14| 6
14:25:16| 21
14:25:18| 37
14:25:20| 21
14:25:22| 9

library(zoo)
library(chron)
z - read.zoo(textConnection(Lines), header = TRUE, sep = |, FUN = times)
plot(z, type = o) # overplot points and lines

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 do I plat timestamed (in seconds) data in R?

2010-07-06 Thread David Winsemius


On Jul 6, 2010, at 9:28 AM, David Winsemius wrote:



On Jul 6, 2010, at 8:58 AM, LosemindL wrote:



Hi all,

I have an Nx2 array, where the first column contains the timestamps  
and the

second column contains the corresponding data.

second  | ts
| --
14:25:00| 18
14:25:02| 14
14:25:04| 11
14:25:06| 4
14:25:08| 24
14:25:10| 13
14:25:12| 12
14:25:14| 6
14:25:16| 21
14:25:18| 37
14:25:20| 21
14:25:22| 9

Since I don't know how to do this in R yet, please allow me to show  
data in

Matlab:


 tz - read.table(textConnection(second  | ts
+ 14:25:00| 18
+ 14:25:02| 14
+ 14:25:04| 11
+ 14:25:06| 4
+ 14:25:08| 24
+ 14:25:10| 13
+ 14:25:12| 12
+ 14:25:14| 6
+ 14:25:16| 21
+ 14:25:18| 37
+ 14:25:20| 21
+ 14:25:22| 9), header=TRUE, sep=|)


In Matlab, the first column becomes:


 0.60069496185
 0.600752314785495
 0.600810185191222
 0.600868055596948
 0.600925925886258
 0.600983796291985
 0.60104197711
 0.601099536987022
 0.601157407392748

How do I plot such data with X-axis showing the timestamps  
(seconds) and the

Y-axis showing the data?




# Should have included the package call and perhaps

install.packages(zoo)
require(zoo)

 plot (zoo(tz[,2], order.by=tz[,1])  )

__


David Winsemius, MD
West Hartford, CT

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


Re: [R] selection of optim parameters

2010-07-06 Thread Prof. John C Nash

Without the data and objective function, it is fairly difficult to tell what is 
going on.

However, we can note:

- no method is specified, but it would have to be L-BFGS-B as this is 
the only
one that can handle box constraints

- the fourth parameter is at a different scale, and you are right to think about 
this, especially as you don't seem to be providing analytic gradients.


Options:
- I would definitely rescale your function i.e., in the code, so your rawpar[4] 
is
1*par[4]. That gets all the scales approx. the same, assuming your start reflects the 
eventual answer

- try bobyqa from minqa package once scaled -- it really doesn't like scale 
differences.

If you have data and a script nicely packaged that can be emailed, I'll take a look at it. 
If analytic gradients can be provided, there are more possibilities too, e.g., Rvmmin. spg 
from BB package also does constraints.


It may, of course, turn out that the solution really is on the boundary.

JN




Message: 47
Date: Mon, 05 Jul 2010 21:53:18 +0200
From: Fabian Gehring fabian.gehr...@bluewin.ch
To: r-help@r-project.org
Subject: [R] selection of optim parameters
Message-ID: 4c32382e.3050...@bluewin.ch
Content-Type: text/plain; charset=ISO-8859-15; format=flowed

Hi all,

I am trying to rebuild the results of a study using a different data
set. I'm using about 450 observations. The code I've written seems to
work well, but I have some troubles minimizing the negative of the
LogLikelyhood function using 5 free parameters.

As starting values I am using the result of the paper I am rebuiling.
The system.time of the calculation of the function is about 0.65 sec.
Since the free parameters should be within some boundaries I am using
the following command:

optim(fn=calculateLogLikelyhood, c(0.4, 2, 0.4, 8, 0.8),
lower=c(0,0,0,0,0), upper=c(1, 5, Inf, Inf, 1), control=list(trace=1,
maxit=1000))

Unfortunately the result doesn't seem to be reasonable. 3 of the
optimized parameters are on the boundaries.

Unfortunately I don't have much experience using optimizatzion methods.
That's why I am asking you.
Do you have any hints for me what should be taken into account when
doing such an optimization.

Is there a good way to implement the boundaries into the code (instead
of doing it while optimizing)? I've read about parscale in the
help-section. Unfortunately I don't really know how to use it. And
anyways, could this help? What other points/controls should be taken
into account?

I know that this might be a bit little information about my current
code. But I don't know what you need to give me some advise. Just let me
know what you need to know.

Thankds




__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Pseudo F statistics with index.G1

2010-07-06 Thread Henrik Aldberg
Thank you Allan,

As i understand it the index.G1 function does not work if one of the
clusters in the partition only contains one object. Is there a way to get
around this in R? In SAS the PSF function seems to ignore the presence of
singleton clusters.


Sincerely

Henrik

[[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] Help With ANOVA (corrected please ignore last email)

2010-07-06 Thread Joshua Wiley
Hi Amit,

When I copy in your data and run

aov(Intensity ~ Group, data = zzzanova)

I get neither the p-value you showed nor the one you expected.  My
suggestions at things to look at would be

1) Where/How did you get the expected p-value?  Another statistics
program (e.g., SPSS or SAS)?  It helps to know where the expected
value came from.  If it was from another program, I would also
recommend searching the R-help archives (for instance for SPSS and
ANOVA)

2) Perhaps something is happening with the unbalanced design (e.g.,
aov() in R handles it differently than you expected)?

3) You only mention that the p-values do not match.  What about other
aspects?  Are the df the same?  I would try to be certain that the
data in R is the same as what was used to calculate the expected
p-value.  summary(zzzanova) will return some nice summary statistics
for each column of your dataframe.

Cheers,


Josh



On Tue, Jul 6, 2010 at 6:12 AM, Amit Patel amitrh...@yahoo.co.uk wrote:
 Sorry i had a misprint in the appendix code in the last email


 Hi I needed some help with ANOVA

 I have a problem with My ANOVA
 analysis. I have a dataset with a known ANOVA p-value, however I can
 not seem to re-create it in R.

 I have created a list (zzzanova) which contains
 1)Intensity Values
 2)Group Number (6 Different Groups)
 3)Sample Number (54 different samples)
 this is created by the script in Appendix 1

 I then conduct ANOVA with the command
 zzz.aov - aov(Intensity ~ Group, data = zzzanova)

 I get a p-value of
 Pr(F)1
 0.9483218

 The
 expected p-value is 0.00490 so I feel I maybe using ANOVA incorrectly
 or have put in a wrong formula. I am trying to do an ANOVA analysis
 across all 6 Groups. Is there something wrong with my formula. But I think I
 have made a mistake in the formula rather than anything else.




 APPENDIX 1

 datalist - c(-4.60517, -4.60517, -4.60517, -4.60517, -4.60517, -4.60517, 
 -4.60517, 3.003749, -4.60517,
    2.045314, 2.482557, -4.60517, -4.60517, -4.60517, -4.60517, 1.592743, 
 -4.60517,
    -4.60517, 0.91328, -4.60517, -4.60517, 1.827744, 2.457795, 0.355075, 
 -4.60517, 2.39127,
    2.016987, 2.319903, 1.146683, -4.60517, -4.60517, -4.60517, 1.846162, 
 -4.60517, 2.121427, 1.973118,
    -4.60517, 2.251568, -4.60517, 2.270724, 0.70338, 0.963816, -4.60517,  
 0.023703, -4.60517,
    2.043382, 1.070586, 2.768289, 1.085169, 0.959334, -0.02428, -4.60517, 
 1.371895, 1.533227)

 zzzanova -
 structure(list(Intensity = datalist,
 Group = structure(c(1,1,1,1,1,1,1,1,1,
         2,2,2,2,2,2,2,2,
         3,3,3,3,3,3,3,3,3,
         4,4,4,4,4,4,4,4,4,4,
         5,5,5,5,5,5,5,5,5,
         6,6,6,6,6,6,6,6,6), .Label = c(Group1, Group2, Group3, 
 Group4, Group5, Group6), class = factor),
    Sample = structure(c( 1, 2, 3, 4, 5, 6, 7, 8, 9,
    10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
    20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
    30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 
 40,41,42,43,44,45,46,47,48,49,50,51,52,53,54)
 ))
 , .Names = c(Intensity,
 Group, Sample), row.names =
 c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
 51, 52, 53, 54),class = data.frame)




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




-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/

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


Re: [R] Linux-Windows problem

2010-07-06 Thread Uwe Ligges



On 06.07.2010 15:24, Henrik Bengtsson wrote:

Sorry for asking the obvious, but have you confirmed that you are
running the same version of R on both systems?



Since the question was about the Editor, I doubt the R version is related.



/H

On Tue, Jul 6, 2010 at 2:11 PM, Bos, Rogerroger@rothschild.com  wrote:

Uwe,

I suspect I might be having a similar problem as the R code I generate
in Windows editor (Tinn-R) doesn't always open properly in my Linux
editor (Rkward), but I don't know how to change the default encoding in
either one.


I am neither using Tinn-R nor Rkward, hence I cannot give any advice here.

Uwe





In R on both machines, getOption(encoding) returns native.enc, so
the problem is not with R, but rather at the OS or editor level.  If you
or anyone could help me out that would be appreciated.

Thanks,

Roger

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
On Behalf Of Uwe Ligges
Sent: Monday, July 05, 2010 10:04 AM
To: Ildiko Varga
Cc: r-help@r-project.org
Subject: Re: [R] Linux-Windows problem



On 05.07.2010 14:31, Ildiko Varga wrote:

Dear All,

I faced the following problem. With the same data.frame the results
are different under Linux and Windows.
Could you help on this topic?


I guess you read in the data differently since you have different
default encodings on both platforms (e.g. latin1 vs. UTF-8) and you data
is probably not plain ASCII.

Best,
Uwe Ligges


Thanks in advance,
Ildiko

Linux:
   d = read.csv(CRP.csv)
   d$drugCode = as.numeric(d$drug)
   cor(d, use=pairwise.complete.obs)
 PATIENT BL.CRP   X24HR.CRP   X48HR.CRP drug   drugCode
PATIENTNA NA  NA  NA   NA NA
BL.CRP NA  1.000  0.84324880 -0.05699590   NA -0.3367147
X24HR.CRP  NA  0.8432488  1. -0.06162383   NA -0.3557316
X48HR.CRP  NA -0.0569959 -0.06162383  1.   NA  0.1553356
drug   NA NA  NA  NA   NA NA
drugCode   NA -0.3367147 -0.35573159  0.15533562   NA  1.000
Warning message:
In cor(d, use = pairwise.complete.obs) : NAs introduced by coercion
   str(d)
'data.frame':   41 obs. of  6 variables:
$ PATIENT  : Factor w/ 41 levels RV13,RV14,..: 2 3 4 6 7 12 13
14
15 17 ...
$ BL.CRP   : num  7.3 31.2 4.2 6.7 1.6 7.7 5.3 38.9 1 7.3 ...
$ X24HR.CRP: num  6.1 24.9 11.1 4.9 1 5 3.7 18 1 7.3 ...
$ X48HR.CRP: num  121.5 40 28.4 34.5 33.3 ...
$ drug : Factor w/ 2 levels active,placebo: 1 1 1 1 1 1 1 1

1

1 ...
$ drugCode : num  1 1 1 1 1 1 1 1 1 1 ...

   Windows:
   d = read.csv(CRP.csv)
   d$drugCode = as.numeric(d$drug)
   cor(d, use=pairwise.complete.obs) Error in cor(d, use =
pairwise.complete.obs) : 'x' must be numeric
   str(d)
'data.frame':   41 obs. of  6 variables:
$ PATIENT  : Factor w/ 41 levels RV13,RV14,..: 2 3 4 6 7 12 13
14
15 17 ...
$ BL.CRP   : num  7.3 31.2 4.2 6.7 1.6 7.7 5.3 38.9 1 7.3 ...
$ X24HR.CRP: num  6.1 24.9 11.1 4.9 1 5 3.7 18 1 7.3 ...
$ X48HR.CRP: num  121.5 40 28.4 34.5 33.3 ...
$ drug : Factor w/ 2 levels active,placebo: 1 1 1 1 1 1 1 1

1

1 ...
$ drugCode : num  1 1 1 1 1 1 1 1 1 1 ...

   [[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.
***

This message is for the named person's use only. It may\...{{dropped:20}}

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


[R] Interpreting NB GLM output - effect sizes?

2010-07-06 Thread Anna Berthinussen

Hi,

I am trying to find out how to interpret the summary output from a neg
bin GLM?

I have 3 significant variables and I can see whether they have a
positive or negative effect, but I can't work out how to calculate the
magnitude of the effect on the mean of the dependent variable. I used
a log link function so I think I might have to use the antilogs of the
coefficients but I have no idea how this relates to the dependent
variable??

Any help would be much appreciated.

My model and output is shown below.

Thanks

Anna

Call:
glm.nb(formula = Pass ~ Dist + Time + Wind, data = bats, link = log,
  init.theta = 0.8510838809)

Deviance Residuals:
  Min   1Q   Median   3Q  Max
-2.2784  -0.9967  -0.3594   0.2603   2.2142

Coefficients:
Estimate Std. Error z value Pr(|z|)
(Intercept)  3.3329718  0.3603909   9.248   2e-16 ***
Dist 0.0008892  0.0002377   3.741 0.000183 ***
Time-0.0159068  0.0034665  -4.589 4.46e-06 ***
Wind-0.1177475  0.0346301  -3.400 0.000673 ***
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

(Dispersion parameter for Negative Binomial(0.8511) family taken to be 1)

  Null deviance: 134.586  on 79  degrees of freedom
Residual deviance:  92.725  on 76  degrees of freedom
AIC: 501.21

Number of Fisher Scoring iterations: 1


Theta:  0.851
Std. Err.:  0.164

   2 x log-likelihood:  -491.211

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] hide ends of line in a density plot

2010-07-06 Thread Albert Vilella
Hi,

(I googled for this answer but didn't find anything)

I am using density plot and I want to trim the ends of the line. eg:

i = rnorm(100,100,2)
j = subset(i,i102  i98)
summary(j)
plot(density(j))

I only want the line to go from 98 to 102. How can I limit the line
(and the axis)
to the values I specify?

Cheers,

Albert.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Profiler for R ? (HFWUtils package)

2010-07-06 Thread Uwe Ligges

or just see

?Rprof

and

?Rprofmem


Uwe Ligges


On 06.07.2010 01:21, Jim Callahan wrote:

Message: 21
Date: Mon, 5 Jul 2010 02:26:29 -0400
From: Ralf Bralf.bie...@gmail.com
To: r-help@r-project.orgr-help@r-project.org
Subject: [R] Profiler for R ?

Hi,

is there such a thing as a profiler for R that informs about a) how
much processing time is used by particular functions and commands and
b) how much memory is used for creating how many objects (or types of
data structures)?


Haven't tried it; but stumbled across Profiling() function in the
HFWUtils package.
Starting at bottom of page 29-30 of HFWUtils package user manual:

profiling
plots tree of execution times

Description
determines how much time a function its and sub-functions (and
sub-functions thereof etc) take to run (‘profiling’). Also draws
picture of this using the interrelations of functions.


HTH,
Jim Callahan
Orlando, FL


In a way I am looking for something similar to the

java profiler (which is started by command line and provides profiling
information collected from the run of a particular program). Is there
such a tool through the R command line or RGUI ? Are there profilers
available for the Eclipse StatET or though another package or
extension?

Thanks,
Ralf


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] hide ends of line in a density plot

2010-07-06 Thread Joshua Wiley
Hello,

If you're just looking to 'zoom in' as it were, this should do it:

plot(density(j), xlim = c(98, 102))

HTH,

Josh

On Tue, Jul 6, 2010 at 8:50 AM, Albert Vilella avile...@gmail.com wrote:
 Hi,

 (I googled for this answer but didn't find anything)

 I am using density plot and I want to trim the ends of the line. eg:

 i = rnorm(100,100,2)
 j = subset(i,i102  i98)
 summary(j)
 plot(density(j))

 I only want the line to go from 98 to 102. How can I limit the line
 (and the axis)
 to the values I specify?

 Cheers,

 Albert.

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




-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/

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


Re: [R] Linux-Windows problem

2010-07-06 Thread Bos, Roger
Its not an R problem, but an editor problem (therefore slightly off-topic I 
admit).  I do most of the coding in Windows and somestimes the Rkward editor 
(in linux) will not open the file and it tells me it cannot determine the 
encoding.  If I try to source that file into R (in the linux version) I 
sometimes get an unexpected character on line such and such, but I don't see 
any problem and the code runs fine in windows.  The fix, I discovered is to 
open the file in notepad (in windows) and resave it.  Then Rkward will open it 
just fine.  

I am sure there is an easy setting I can change to make the two work more 
smoothly together, but I just haven't found it yet.

Thanks,

Roger


-Original Message-
From: henrik.bengts...@gmail.com [mailto:henrik.bengts...@gmail.com] On Behalf 
Of Henrik Bengtsson
Sent: Tuesday, July 06, 2010 9:24 AM
To: Bos, Roger
Cc: Uwe Ligges; Ildiko Varga; r-help
Subject: Re: [R] Linux-Windows problem

Sorry for asking the obvious, but have you confirmed that you are running the 
same version of R on both systems?

/H

On Tue, Jul 6, 2010 at 2:11 PM, Bos, Roger roger@rothschild.com wrote:
 Uwe,

 I suspect I might be having a similar problem as the R code I generate 
 in Windows editor (Tinn-R) doesn't always open properly in my Linux 
 editor (Rkward), but I don't know how to change the default encoding 
 in either one.

 In R on both machines, getOption(encoding) returns native.enc, so 
 the problem is not with R, but rather at the OS or editor level.  If 
 you or anyone could help me out that would be appreciated.

 Thanks,

 Roger

 -Original Message-
 From: r-help-boun...@r-project.org 
 [mailto:r-help-boun...@r-project.org]
 On Behalf Of Uwe Ligges
 Sent: Monday, July 05, 2010 10:04 AM
 To: Ildiko Varga
 Cc: r-help@r-project.org
 Subject: Re: [R] Linux-Windows problem



 On 05.07.2010 14:31, Ildiko Varga wrote:
 Dear All,

 I faced the following problem. With the same data.frame the results 
 are different under Linux and Windows.
 Could you help on this topic?

 I guess you read in the data differently since you have different 
 default encodings on both platforms (e.g. latin1 vs. UTF-8) and you 
 data is probably not plain ASCII.

 Best,
 Uwe Ligges

 Thanks in advance,
 Ildiko

 Linux:
 d = read.csv(CRP.csv)
 d$drugCode = as.numeric(d$drug)
 cor(d, use=pairwise.complete.obs)
 PATIENT BL.CRP   X24HR.CRP   X48HR.CRP drug   
 drugCode PATIENTNA NA  NA  NA   NA 
 NA BL.CRP NA  1.000  0.84324880 -0.05699590   NA 
 -0.3367147 X24HR.CRP  NA  0.8432488  1. -0.06162383   NA 
 -0.3557316 X48HR.CRP  NA -0.0569959 -0.06162383  1.   NA  
 0.1553356 drug   NA NA  NA  NA   NA 
 NA drugCode   NA -0.3367147 -0.35573159  0.15533562   NA  
 1.000 Warning message:
 In cor(d, use = pairwise.complete.obs) : NAs introduced by coercion
 str(d)
 'data.frame':   41 obs. of  6 variables:
$ PATIENT  : Factor w/ 41 levels RV13,RV14,..: 2 3 4 6 7 12 13
 14
 15 17 ...
$ BL.CRP   : num  7.3 31.2 4.2 6.7 1.6 7.7 5.3 38.9 1 7.3 ...
$ X24HR.CRP: num  6.1 24.9 11.1 4.9 1 5 3.7 18 1 7.3 ...
$ X48HR.CRP: num  121.5 40 28.4 34.5 33.3 ...
$ drug : Factor w/ 2 levels active,placebo: 1 1 1 1 1 1 1 
 1
 1
 1 ...
$ drugCode : num  1 1 1 1 1 1 1 1 1 1 ...

   Windows:
 d = read.csv(CRP.csv)
 d$drugCode = as.numeric(d$drug)
 cor(d, use=pairwise.complete.obs) Error in cor(d, use =
 pairwise.complete.obs) : 'x' must be numeric
 str(d)
 'data.frame':   41 obs. of  6 variables:
$ PATIENT  : Factor w/ 41 levels RV13,RV14,..: 2 3 4 6 7 12 13
 14
 15 17 ...
$ BL.CRP   : num  7.3 31.2 4.2 6.7 1.6 7.7 5.3 38.9 1 7.3 ...
$ X24HR.CRP: num  6.1 24.9 11.1 4.9 1 5 3.7 18 1 7.3 ...
$ X48HR.CRP: num  121.5 40 28.4 34.5 33.3 ...
$ drug : Factor w/ 2 levels active,placebo: 1 1 1 1 1 1 1 
 1
 1
 1 ...
$ drugCode : num  1 1 1 1 1 1 1 1 1 1 ...

   [[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.
 ***

 This message is for the named person's use only. It 
 may\...{{dropped:20}}

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

Re: [R] Interpreting NB GLM output - effect sizes?

2010-07-06 Thread Achim Zeileis

On Tue, 6 Jul 2010, Anna Berthinussen wrote:


Hi,

I am trying to find out how to interpret the summary output from a neg
bin GLM?

I have 3 significant variables and I can see whether they have a
positive or negative effect, but I can't work out how to calculate the
magnitude of the effect on the mean of the dependent variable. I used
a log link function so I think I might have to use the antilogs of the
coefficients but I have no idea how this relates to the dependent
variable??


The mean equation is

  log(mu) = x'b

so this is similar in interpretation to a semi-logarithmic linear model. 
Absolute changes in x lead to relative changes in the response. In your 
example below, a sloppy formulation would be: If Time increases by one 
unit, the expected mean Pass decreases by 1.6% (if everything else stays 
the same).


A useful discussion of this is for example in Analysis of Microdata by 
Winkelmann  Boes (2009, Springer). But of course in many other textbooks 
as well.


Another useful approach is to employ effects to visualize the effects, 
e.g.:


  library(effects)
  plot(allEffects(fitted_glm.nb_object), ask = FALSE, rescale = FALSE)

hth,
Z


Any help would be much appreciated.

My model and output is shown below.

Thanks

Anna

Call:
glm.nb(formula = Pass ~ Dist + Time + Wind, data = bats, link = log,
init.theta = 0.8510838809)

Deviance Residuals:
Min   1Q   Median   3Q  Max
-2.2784  -0.9967  -0.3594   0.2603   2.2142

Coefficients:
  Estimate Std. Error z value Pr(|z|)
(Intercept)  3.3329718  0.3603909   9.248   2e-16 ***
Dist 0.0008892  0.0002377   3.741 0.000183 ***
Time-0.0159068  0.0034665  -4.589 4.46e-06 ***
Wind-0.1177475  0.0346301  -3.400 0.000673 ***
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

(Dispersion parameter for Negative Binomial(0.8511) family taken to be 1)

Null deviance: 134.586  on 79  degrees of freedom
Residual deviance:  92.725  on 76  degrees of freedom
AIC: 501.21

Number of Fisher Scoring iterations: 1


  Theta:  0.851
  Std. Err.:  0.164

 2 x log-likelihood:  -491.211

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

2010-07-06 Thread Joris Meys
We're missing the samp1 etc. in order to be able to test the code.
Where did you get the other p-value?
Cheers
Joris

On Tue, Jul 6, 2010 at 3:08 PM, Amit Patel amitrh...@yahoo.co.uk wrote:
 Hi I needed some help with ANOVA

 I have a problem with My ANOVA
 analysis. I have a dataset with a known ANOVA p-value, however I can
 not seem to re-create it in R.

 I have created a list (zzzanova) which contains
 1)Intensity Values
 2)Group Number (6 Different Groups)
 3)Sample Number (54 different samples)
 this is created by the script in Appendix 1

 I then conduct ANOVA with the command
 zzz.aov - aov(Intensity ~ Group, data = zzzanova)

 I get a p-value of
 Pr(F)1
 0.9483218

 The
 expected p-value is 0.00490 so I feel I maybe using ANOVA incorrectly
 or have put in a wrong formula. I am trying to do an ANOVA analysis
 across all 6 Groups. Is there something wrong with my formula. But I think I
 have made a mistake in the formula rather than anything else.




 APPENDIX 1

 datalist - c(-4.60517, -4.60517, -4.60517, -4.60517, -4.60517, -4.60517, 
 -4.60517, 3.003749, -4.60517,
    2.045314, 2.482557, -4.60517, -4.60517, -4.60517, -4.60517, 1.592743, 
 -4.60517,
    -4.60517, 0.91328, -4.60517, -4.60517, 1.827744, 2.457795, 0.355075, 
 -4.60517, 2.39127,
    2.016987, 2.319903, 1.146683, -4.60517, -4.60517, -4.60517, 1.846162, 
 -4.60517, 2.121427, 1.973118,
    -4.60517, 2.251568, -4.60517, 2.270724, 0.70338, 0.963816, -4.60517,  
 0.023703, -4.60517,
    2.043382, 1.070586, 2.768289, 1.085169, 0.959334, -0.02428, -4.60517, 
 1.371895, 1.533227)

 zzzanova -
 structure(list(Intensity = c(t(Samp1), t(Samp2), t(Samp3), t(Samp4)),
 Group = structure(c(1,1,1,1,1,1,1,1,1,
         2,2,2,2,2,2,2,2,
         3,3,3,3,3,3,3,3,3,
         4,4,4,4,4,4,4,4,4,4,
         5,5,5,5,5,5,5,5,5,
         6,6,6,6,6,6,6,6,6), .Label = c(Group1, Group2, Group3, 
 Group4, Group5, Group6), class = factor),
    Sample = structure(c( 1, 2, 3, 4, 5, 6, 7, 8, 9,
    10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
    20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
    30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 
 40,41,42,43,44,45,46,47,48,49,50,51,52,53,54)
 ))
 , .Names = c(Intensity,
 Group, Sample), row.names =
 c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
 51, 52, 53, 54),class = data.frame)




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




-- 
Joris Meys
Statistical consultant

Ghent University
Faculty of Bioscience Engineering
Department of Applied mathematics, biometrics and process control

tel : +32 9 264 59 87
joris.m...@ugent.be
---
Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php

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

2010-07-06 Thread Amanda Peterson
Hi,

I'm trying to install R on a Solaris 10 machine.
username -a
SunOS discovery01 5.10 Generic_142900-09 sun4v sparc SUNW, T5440

I get the following error when running make install-tests. Any suggestions?

r...@discovery01:/opt/R-2.10.1 make install-tests installing specific
tests mkdir /usr/local/lib/R.test mkdir /usr/local/lib/tests/Packages
installing package tests ...
mkdir /usr/local/lib/R/library
mkdir /usr/local/lib/R/library/tools
mkdir /usr/local/lib/R/library/tools/tests
mkdir /usr/local/lib/R/library/utils
mkdir /usr/local/lib/R/libraryutils/tests
mkdir /usr/local/lib/R/library/grDevices
mkdir /usr/local/lib/R/library/grDevices/tests
mkdir /usr/local/lib/R/library/stats
mkdir /usr/local/lib/R/library/stats/tests
mkdir /usr/local/lib/R/library/grid
mkdir /usr/local/lib/R/library/grid/tests
mkdir /usr/local/lib/R/library/boot
mkdir /usr/local/lib/R/library/boot/tests
/opt/sfw/bin/install: omitting directory 'boot/tests/Examples'
*** Error code 1
The following command caused the error:
(tmp=${TMPDIR-/tmp}/R$$; mkdir ${tmp}; \  cd ../../..;
abs_top_srcdir='pwd'; \  abs_MKINSTALLDIRS='echo /bin/bash
../../../src/scripts/mkinstalldirs.in | sed
s:../../..:${abs_top_srcdir}:' \  cd ${tmp}; \  for pkg in boot cluster
codetools foreign KernSmooth lattice mgcv nlme rpart survival MASS class
nnet spatial Matrix; do \
  gzip -dc /opt/R-2.10.1/src/library/Recommended/${pkg}.tgz |
//usr/sfw/bin/gtar xf - ; \ done ; \ for pkg in boot cluster codetools
foreign KernSmoth lattice mgcv nlme rpart survival MASS class nnet
spatial Matrix; do \
  if test -d ${pkg}/tests; then \
${abs_MKINSTALLDIRS} /usr/local/lib/R/library/${pkg}/tests \
for f in ${pkg}/tests/*; do \
  /opt/sfw/bin/install -c -m 644 ${f}
/usr/local/lib/R/library/${pkg}/tests; \
done; \
  fi;  \
 done)
make: Fatal error: Command failed for target 'install-tests'
Current working directory /opt/R-2.10.1/src/library/Recommended
*** Error code 1
The folloiwng command caused the error:
(cd Recommended  make install-tests)
make: Fatal error: Command railed for target 'install-tests'
Current working directory /opt/R-2.10.1/src/library
*** Error code 1 (ignored)
The following command caused the error:
(cd src/library  make install-tests)


Thanks,
Amanda

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Rcmdr installation under Unbuntu installatiion errors

2010-07-06 Thread John Sorkin
Unbuntu 10.04
R 2.10

I am trying to install Rcmdr and receive the following messages:

The downloaded packages are in
‘/tmp/RtmpzhjDZG/downloaded_packages’
Warning messages:
1: In install.packages(Rcmdr, dependencies = TRUE) :
  installation of package 'Rmpi' had non-zero exit status
2: In install.packages(Rcmdr, dependencies = TRUE) :
  installation of package 'rpvm' had non-zero exit status
3: In install.packages(Rcmdr, dependencies = TRUE) :
  installation of package 'rsprng' had non-zero exit status
4: In install.packages(Rcmdr, dependencies = TRUE) :
  installation of package 'tkrplot' had non-zero exit status
5: In install.packages(Rcmdr, dependencies = TRUE) :
  installation of package 'rgl' had non-zero exit status
6: In install.packages(Rcmdr, dependencies = TRUE) :
  installation of package 'RODBC' had non-zero exit status
7: In install.packages(Rcmdr, dependencies = TRUE) :
  installation of package 'shapes' had non-zero exit status
8: In install.packages(Rcmdr, dependencies = TRUE) :
  installation of package 'latentnet' had non-zero exit status
9: In install.packages(Rcmdr, dependencies = TRUE) :
  installation of package 'VIM' had non-zero exit status
10: In install.packages(Rcmdr, dependencies = TRUE) :
  installation of package 'statnet' had non-zero exit status

Do I need to be concerned? Do I need to take any action?

Thank you,
John 


John Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
Baltimore VA Medical Center GRECC,
University of Maryland School of Medicine Claude D. Pepper OAIC,
University of Maryland Clinical Nutrition Research Unit, and
Baltimore VA Center Stroke of Excellence

University of Maryland School of Medicine
Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524

(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)
jsor...@grecc.umaryland.edu

Confidentiality Statement:
This email message, including any attachments, is for th...{{dropped:6}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Profiler for R ? (HFWUtils package)

2010-07-06 Thread Hadley Wickham
And the profr package for an alternative display.
Hadley

On Tuesday, July 6, 2010, Uwe Ligges lig...@statistik.tu-dortmund.de wrote:
 or just see

 ?Rprof

 and

 ?Rprofmem


 Uwe Ligges


 On 06.07.2010 01:21, Jim Callahan wrote:

 Message: 21
 Date: Mon, 5 Jul 2010 02:26:29 -0400
 From: Ralf Bralf.bie...@gmail.com
 To: r-help@r-project.orgr-help@r-project.org
 Subject: [R] Profiler for R ?

 Hi,

 is there such a thing as a profiler for R that informs about a) how
 much processing time is used by particular functions and commands and
 b) how much memory is used for creating how many objects (or types of
 data structures)?


 Haven't tried it; but stumbled across Profiling() function in the
 HFWUtils package.
 Starting at bottom of page 29-30 of HFWUtils package user manual:

 profiling
 plots tree of execution times

 Description
 determines how much time a function its and sub-functions (and
 sub-functions thereof etc) take to run (‘profiling’). Also draws
 picture of this using the interrelations of functions.


 HTH,
 Jim Callahan
 Orlando, FL


 In a way I am looking for something similar to the

 java profiler (which is started by command line and provides profiling
 information collected from the run of a particular program). Is there
 such a tool through the R command line or RGUI ? Are there profilers
 available for the Eclipse StatET or though another package or
 extension?

 Thanks,
 Ralf


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


-- 
Assistant Professor / Dobelman Family Junior Chair
Department of Statistics / Rice University
http://had.co.nz/

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


[R] Error message using mi() in mi package

2010-07-06 Thread Andrew Miles

Hello!

I get the following message when I run the mi() function from the mi  
package.


Error while imputing variable: c3 , model: mi.polr
 Error in eval(expr, envir, enclos) : could not find function  
c14ordered


Here's the situation:

I am running R v. 2.9.2 on Mac OSX v. 10.5.8.  I am trying to impute  
missing data in a data set that I've trimmed down to 302 variables.   
I've created an mi.info() object on the data, and I've updated the  
type of variable where necessary so that the mi() imputing function  
uses the most appropriate type of models to estimate imputed values.   
The data have no collinearlity.  I then run the mi function like this:


imp = mi(imp.data, info=info2, n.iter=10)

where imp.data is my data set of 302 variables and info2 is my  
modified mi.info object.  I get the message as posted above.  The  
message only occurs when working on a variable I have labeled as  
ordered-categorical.  But the mi() function processes most variables  
labeled as ordered-categorical just fine.  In fact, if shrink my  
data set (say, to just 5 variables) I can get mi() to process a  
problematic variable just fine as well.


I'm not sure what the function c14ordered is that the error message  
refers to.  My first thought is maybe it is referring to one of my  
variables in my data?  Variables names in my data follow a basic  
letter-number pattern (i.e. a1, a2, etc.), but there is no c14, rather  
c14a1, c14a2, etc.  So I'm not sure if the variable has anything to do  
with the problem, but I thought I'd mention it just in case someone  
wiser in this matter than I can see something I cannot.


I cannot post code that reproduces the problem due to the nature of  
the code and data involved.


Any help would be appreciated, as I am not sure what is happening, and  
can't see why I can sometimes impute a variable labeled as ordered- 
categorical and sometimes cannot.


Thanks!

Andrew Miles

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

2010-07-06 Thread beloitstudent

Thanks for the advice!  It has worked for the most part.  However, I am
still coming up with an error message when placing my break line in the axis
that I'm not sure what it means.  If you could help me out, that would be
fantastic...otherwise I might just have to see if I can add it on
powerpoint.  Here is the code you gave me and my script that doesn't work.
*
*#your script
*axis.break-function(axis=1,*
*breakpos=NULL,pos=NA,bgcol=**white,breakcol=black,
 style=slash,brw=0.02) {

 # get the coordinates of the outside of the plot
 figxy-par(usr)
 # flag if either axis is logarithmic
 xaxl-par(xlog)
 yaxl-par(ylog)
 # calculate the x and y offsets for the break
 xw-(figxy[2]-figxy[1])*brw
 yw-(figxy[4]-figxy[3])*brw
 if(!is.na(pos)) figxy-rep(pos,4)
 # if no break position was given, put it just off the plot origin
 if(is.null(breakpos))
  breakpos-ifelse(axis%%2,**figxy[1]+xw*2,figxy[3]+yw*2)
 if(xaxl  (axis == 1 || axis == 3)) breakpos-log10(breakpos)
 if(yaxl  (axis == 2 || axis == 4)) breakpos-log10(breakpos)
 # set up the blank rectangle (left, bottom, right, top)
 switch(axis,
  br-c(breakpos-xw/2,figxy[3]-**yw/2,breakpos+xw/2,figxy[3]+**yw/2),
  br-c(figxy[1]-xw/2,breakpos-**yw/2,figxy[1]+xw/2,breakpos+**yw/2),
  br-c(breakpos-xw/2,figxy[4]-**yw/2,breakpos+xw/2,figxy[4]+**yw/2),
  br-c(figxy[2]-xw/2,breakpos-**yw/2,figxy[2]+xw/2,breakpos+**yw/2),
  stop(Improper axis specification.))
 # get the current setting of xpd
 old.xpd-par(xpd)
 # don't cut the break off at the edge of the plot
 par(xpd=TRUE)
 # correct for logarithmic axes
 if(xaxl) br[c(1,3)]-10^br[c(1,3)]
 if(yaxl) br[c(2,4)]-10^br[c(2,4)]
 if(style == gap) {
  if(xaxl) {
   figxy[1]-10^figxy[1]
   figxy[2]-10^figxy[2]
  }
  if(yaxl) {
   figxy[3]-10^figxy[3]
   figxy[4]-10^figxy[4]
  }
  # blank out the gap area and calculate the line segments
  if(axis == 1 || axis == 3) {
   rect(breakpos,figxy[3],**breakpos+xw,figxy[4],col=**bgcol,border=bgcol)
   xbegin-c(breakpos,breakpos+**xw)
   ybegin-c(figxy[3],figxy[3])
   xend-c(breakpos,breakpos+xw)
   yend-c(figxy[4],figxy[4])
   if(xaxl) {
xbegin-10^xbegin
xend-10^xend
   }
  }
  else {
   rect(figxy[1],breakpos,figxy[**2],breakpos+yw,col=bgcol,**border=bgcol)
   xbegin-c(figxy[1],figxy[1])
   ybegin-c(breakpos,breakpos+**yw)
   xend-c(figxy[2],figxy[2])
   yend-c(breakpos,breakpos+yw)
   if(xaxl) {
xbegin-10^xbegin
xend-10^xend
   }
  }
  # clip the lines
  par(xpd=TRUE)
 }
 else {
  # draw the blank rectangle
  rect(br[1],br[2],br[3],br[4],**col=bgcol,border=bgcol)
  if(style == slash) {
   # calculate the slash ends
   if(axis == 1 || axis == 3) {
xbegin-c(breakpos-xw,**breakpos)
xend-c(breakpos,breakpos+xw)
ybegin-c(br[2],br[2])
yend-c(br[4],br[4])
if(xaxl) {
 xbegin-10^xbegin
 xend-10^xend
}
   }
   else {
xbegin-c(br[1],br[1])
xend-c(br[3],br[3])
ybegin-c(breakpos-yw,**breakpos)
yend-c(breakpos,breakpos+yw)
if(yaxl) {
 ybegin-10^ybegin
 yend-10^yend
}
   }
  }
  else {
   # calculate the zigzag ends
   if(axis == 1 || axis == 3) {
xbegin-c(breakpos-xw/2,**breakpos-xw/4,breakpos+xw/4)
xend-c(breakpos-xw/4,**breakpos+xw/4,breakpos+xw/2)

ybegin-c(ifelse(yaxl,10^**figxy[3+(axis==3)],figxy[3+(**axis==3)]),br[4],br[2])


yend-c(br[4],br[2],ifelse(**yaxl,10^figxy[3+(axis==3)],**figxy[3+(axis==3)]))

if(xaxl) {
 xbegin-10^xbegin
 xend-10^xend
}
   }
   else {

xbegin-c(ifelse(xaxl,10^**figxy[1+(axis==4)],figxy[1+(**axis==4)]),br[1],br[3])


xend-c(br[1],br[3],ifelse(**xaxl,10^figxy[1+(axis==4)],**figxy[1+(axis==4)]))

ybegin-c(breakpos-yw/2,**breakpos-yw/4,breakpos+yw/4)
yend-c(breakpos-yw/4,**breakpos+yw/4,breakpos+yw/2)
if(yaxl) {
 ybegin-10^ybegin
 yend-10^yend
}
   }
  }
 }
 # draw the segments
 segments(xbegin,ybegin,xend,**yend,col=breakcol,lty=1)
 # restore xpd
 par(xpd=FALSE)
} *

#my script
*library(plotrix)

par(mar=c(6,8,4,4))
par(xpd=TRUE)

Saline - structure(list(Time = c(-20L, 0, 30L, 45L, 60L, 80L, 110L,
140L,200L, 260L, 320L), Average =
c(119.250,118.750,117.500,132.75,151.875,159.75,142.75,160,168,167.125,143),SEM
=
c(2.211,2.569,2.665,5.435146394,6.208741369,8.363550657,8.51349469,14.30284687,15.93865792,16.76541326,13.796)),
.Names = c(Time (min),
Arterial Plasma Glucose (µg/mL), SEM), class = data.frame, row.names =
c(1, 2,
3, 4, 5, 6, 7, 8, 9, 10, 11))

plotCI(x=Saline [,1],y=Saline [,2]-80, uiw=Saline [,3], err=y,
  pt.bg=par(bg),pch=19, cex=2.5 ,gap=0, sfrac=0.005,
 xlim=c(-20,340),xaxp=c(-20,320,12), xlab=Time (min), ylim=c(0,100),
  yaxp=c(0,100,10), ylab=Arterial Plasma\nGlucose (µg/mL), las=1,
axes=FALSE,
 font.lab=2.2,cex.lab=1.6)

 Ex - structure(list(Time = c(-20L, 0, 30L, 45L, 60L, 80L, 110L, 140L,
200L, 260L, 320L), Average =
c(117.500,117.625,117.375,134.5,166.25,173.5,164.25,162.5,160.375,150.25,139.875),
SEM =

Re: [R] help with predict.lda

2010-07-06 Thread Changbin Du
Thanks all so much for your help! I went out for 2 days vacation and could
not reply your guys email. Yes, the CV=False works.

Thanks again!




On Sun, Jul 4, 2010 at 2:47 AM, Peter Ehlers ehl...@ucalgary.ca wrote:

 On 2010-07-03 21:33, Changbin Du wrote:

 HI, Dear community,

 I am using the linear discriminant analysis to build model and make new
 predictions:

  dim(train)  #training data

 [1] 1272   22

 dim(valid)  # validation data

 [1] 140  22


 lda.fit- lda(out ~ ., data=train, na.action=na.omit, CV=TRUE) # model
 fitting of linear discriminant analysis on training data

  predict(lda.fit, valid)   # make prediction on validation data

 Error in UseMethod(predict) :
   no applicable method for 'predict' applied to an object of class list

 Can anyone help with this?


 Suggestions:

 0. lda() is not a function in the base installation of R;
   I'll assume that you mean MASS::lda.

 1. ask yourself what you expect from setting CV = TRUE.

 2. *carefully* (re)read the help pages for lda (especially
   the Value section) and for predict.lda (the 'object'
   definition).

 3. use CV = FALSE

 4. whenever problems arise, *first* use str() on your
   objects to see what you've got.

 5. finally, do provide *reproducible* code; here, you
   could have used the example on the help page.

  -Peter Ehlers


 Thanks so much!




-- 
Sincerely,
Changbin
--

[[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] ctree ordering nodes

2010-07-06 Thread Paras Sharma
Hello,

 

When using the ctree function, from library (party) what is the syntax to
order the 

Variables in the nodes in a specific way?

 

For example, how would I specify to make a binary come first, then a
continuous variable?

 

Also is there a way to force ctree to show variables which are not
significant?

 

Thanks,

Paras

 

 

 


[[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] Linux-Windows problem

2010-07-06 Thread Duncan Murdoch

On 06/07/2010 11:54 AM, Bos, Roger wrote:
Its not an R problem, but an editor problem (therefore slightly off-topic I admit).  I do most of the coding in Windows and somestimes the Rkward editor (in linux) will not open the file and it tells me it cannot determine the encoding.  If I try to source that file into R (in the linux version) I sometimes get an unexpected character on line such and such, but I don't see any problem and the code runs fine in windows.  The fix, I discovered is to open the file in notepad (in windows) and resave it.  Then Rkward will open it just fine.  


I am sure there is an easy setting I can change to make the two work more 
smoothly together, but I just haven't found it yet.
  


If you avoid non-ascii characters, you shouldn't have problems.  (Ascii 
is a subset of most of the encodings used on Windows and Linux.)  The 
only likely problem would come if you save in Unicode (UCS-2), which 
your Linux system may have trouble recognizing.


Duncan Murdoch

Thanks,

Roger


-Original Message-
From: henrik.bengts...@gmail.com [mailto:henrik.bengts...@gmail.com] On Behalf 
Of Henrik Bengtsson
Sent: Tuesday, July 06, 2010 9:24 AM
To: Bos, Roger
Cc: Uwe Ligges; Ildiko Varga; r-help
Subject: Re: [R] Linux-Windows problem

Sorry for asking the obvious, but have you confirmed that you are running the 
same version of R on both systems?

/H

On Tue, Jul 6, 2010 at 2:11 PM, Bos, Roger roger@rothschild.com wrote:
 Uwe,

 I suspect I might be having a similar problem as the R code I generate 
 in Windows editor (Tinn-R) doesn't always open properly in my Linux 
 editor (Rkward), but I don't know how to change the default encoding 
 in either one.


 In R on both machines, getOption(encoding) returns native.enc, so 
 the problem is not with R, but rather at the OS or editor level.  If 
 you or anyone could help me out that would be appreciated.


 Thanks,

 Roger

 -Original Message-
 From: r-help-boun...@r-project.org 
 [mailto:r-help-boun...@r-project.org]

 On Behalf Of Uwe Ligges
 Sent: Monday, July 05, 2010 10:04 AM
 To: Ildiko Varga
 Cc: r-help@r-project.org
 Subject: Re: [R] Linux-Windows problem



 On 05.07.2010 14:31, Ildiko Varga wrote:
 Dear All,

 I faced the following problem. With the same data.frame the results 
 are different under Linux and Windows.

 Could you help on this topic?

 I guess you read in the data differently since you have different 
 default encodings on both platforms (e.g. latin1 vs. UTF-8) and you 
 data is probably not plain ASCII.


 Best,
 Uwe Ligges

 Thanks in advance,
 Ildiko

 Linux:
 d = read.csv(CRP.csv)
 d$drugCode = as.numeric(d$drug)
 cor(d, use=pairwise.complete.obs)
 PATIENT BL.CRP   X24HR.CRP   X48HR.CRP drug   
 drugCode PATIENTNA NA  NA  NA   NA 
 NA BL.CRP NA  1.000  0.84324880 -0.05699590   NA 
 -0.3367147 X24HR.CRP  NA  0.8432488  1. -0.06162383   NA 
 -0.3557316 X48HR.CRP  NA -0.0569959 -0.06162383  1.   NA  
 0.1553356 drug   NA NA  NA  NA   NA 
 NA drugCode   NA -0.3367147 -0.35573159  0.15533562   NA  
 1.000 Warning message:

 In cor(d, use = pairwise.complete.obs) : NAs introduced by coercion
 str(d)
 'data.frame':   41 obs. of  6 variables:
$ PATIENT  : Factor w/ 41 levels RV13,RV14,..: 2 3 4 6 7 12 13
 14
 15 17 ...
$ BL.CRP   : num  7.3 31.2 4.2 6.7 1.6 7.7 5.3 38.9 1 7.3 ...
$ X24HR.CRP: num  6.1 24.9 11.1 4.9 1 5 3.7 18 1 7.3 ...
$ X48HR.CRP: num  121.5 40 28.4 34.5 33.3 ...
$ drug : Factor w/ 2 levels active,placebo: 1 1 1 1 1 1 1 
 1

 1
 1 ...
$ drugCode : num  1 1 1 1 1 1 1 1 1 1 ...

   Windows:
 d = read.csv(CRP.csv)
 d$drugCode = as.numeric(d$drug)
 cor(d, use=pairwise.complete.obs) Error in cor(d, use =
 pairwise.complete.obs) : 'x' must be numeric
 str(d)
 'data.frame':   41 obs. of  6 variables:
$ PATIENT  : Factor w/ 41 levels RV13,RV14,..: 2 3 4 6 7 12 13
 14
 15 17 ...
$ BL.CRP   : num  7.3 31.2 4.2 6.7 1.6 7.7 5.3 38.9 1 7.3 ...
$ X24HR.CRP: num  6.1 24.9 11.1 4.9 1 5 3.7 18 1 7.3 ...
$ X48HR.CRP: num  121.5 40 28.4 34.5 33.3 ...
$ drug : Factor w/ 2 levels active,placebo: 1 1 1 1 1 1 1 
 1

 1
 1 ...
$ drugCode : num  1 1 1 1 1 1 1 1 1 1 ...

   [[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] Help With ANOVA

2010-07-06 Thread Joris Meys
I still can't reproduce your example. The aov output gives me the following :

 anova(aov(Intensity ~ Group, data = zzzanova))
Analysis of Variance Table

Response: Intensity
  Df Sum Sq Mean Sq F value  Pr(F)
Group  5  98.85  19.771  2.1469 0.07576 .
Residuals 48 442.03   9.209
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Next to that I noticed you have truncated data, which has implications
for the analysis as well. If you use a Kruskal-Wallis test, the
p-value becomes larger :

 kruskal.test(Intensity ~ Group, data = zzzanova)

Kruskal-Wallis rank sum test

data:  Intensity by Group
Kruskal-Wallis chi-squared = 6.6955, df = 5, p-value = 0.2443

Which is to be expected, as you have almost 50% truncated data. So a
p-value of 0.005 seems very wrong to me.

Cheers
Joris

On Tue, Jul 6, 2010 at 7:11 PM, Amit Patel amitrh...@yahoo.co.uk wrote:
 Hi Joris

 Sorry i had a misprint in the appendix code in the last email

 datalist - c(-4.60517, -4.60517, -4.60517, -4.60517, -4.60517, -4.60517, 
 -4.60517, 3.003749, -4.60517,
    2.045314, 2.482557, -4.60517, -4.60517, -4.60517, -4.60517, 1.592743, 
 -4.60517,
    -4.60517, 0.91328, -4.60517, -4.60517, 1.827744, 2.457795, 0.355075, 
 -4.60517, 2.39127,
    2.016987, 2.319903, 1.146683, -4.60517, -4.60517, -4.60517, 1.846162, 
 -4.60517, 2.121427, 1.973118,
    -4.60517, 2.251568, -4.60517, 2.270724, 0.70338, 0.963816, -4.60517,  
 0.023703, -4.60517,
    2.043382, 1.070586, 2.768289, 1.085169, 0.959334, -0.02428, -4.60517, 
 1.371895, 1.533227)

 zzzanova -
 structure(list(Intensity = datalist,
 Group = structure(c(1,1,1,1,1,1,1,1,1,
         2,2,2,2,2,2,2,2,
         3,3,3,3,3,3,3,3,3,
         4,4,4,4,4,4,4,4,4,4,
         5,5,5,5,5,5,5,5,5,
         6,6,6,6,6,6,6,6,6), .Label = c(Group1, Group2, Group3, 
 Group4, Group5, Group6), class = factor),
    Sample = structure(c( 1, 2, 3, 4, 5, 6, 7, 8, 9,
    10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
    20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
    30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 
 40,41,42,43,44,45,46,47,48,49,50,51,52,53,54)
 ))
 , .Names = c(Intensity,
 Group, Sample), row.names =
 c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
 51, 52, 53, 54),class = data.frame)


 Thanks for your reply





 - Original Message 
 From: Joris Meys jorism...@gmail.com
 To: Amit Patel amitrh...@yahoo.co.uk
 Cc: r-help@r-project.org
 Sent: Tue, 6 July, 2010 17:04:40
 Subject: Re: [R] Help With ANOVA

 We're missing the samp1 etc. in order to be able to test the code.
 Where did you get the other p-value?
 Cheers
 Joris

 On Tue, Jul 6, 2010 at 3:08 PM, Amit Patel amitrh...@yahoo.co.uk wrote:
 Hi I needed some help with ANOVA

 I have a problem with My ANOVA
 analysis. I have a dataset with a known ANOVA p-value, however I can
 not seem to re-create it in R.

 I have created a list (zzzanova) which contains
 1)Intensity Values
 2)Group Number (6 Different Groups)
 3)Sample Number (54 different samples)
 this is created by the script in Appendix 1

 I then conduct ANOVA with the command
 zzz.aov - aov(Intensity ~ Group, data = zzzanova)

 I get a p-value of
 Pr(F)1
 0.9483218

 The
 expected p-value is 0.00490 so I feel I maybe using ANOVA incorrectly
 or have put in a wrong formula. I am trying to do an ANOVA analysis
 across all 6 Groups. Is there something wrong with my formula. But I think I
 have made a mistake in the formula rather than anything else.




 APPENDIX 1

 datalist - c(-4.60517, -4.60517, -4.60517, -4.60517, -4.60517, -4.60517, 
 -4.60517, 3.003749, -4.60517,
    2.045314, 2.482557, -4.60517, -4.60517, -4.60517, -4.60517, 1.592743, 
 -4.60517,
    -4.60517, 0.91328, -4.60517, -4.60517, 1.827744, 2.457795, 0.355075, 
 -4.60517, 2.39127,
    2.016987, 2.319903, 1.146683, -4.60517, -4.60517, -4.60517, 1.846162, 
 -4.60517, 2.121427, 1.973118,
    -4.60517, 2.251568, -4.60517, 2.270724, 0.70338, 0.963816, -4.60517,  
 0.023703, -4.60517,
    2.043382, 1.070586, 2.768289, 1.085169, 0.959334, -0.02428, -4.60517, 
 1.371895, 1.533227)

 zzzanova -
 structure(list(Intensity = c(t(Samp1), t(Samp2), t(Samp3), t(Samp4)),
 Group = structure(c(1,1,1,1,1,1,1,1,1,
         2,2,2,2,2,2,2,2,
         3,3,3,3,3,3,3,3,3,
         4,4,4,4,4,4,4,4,4,4,
         5,5,5,5,5,5,5,5,5,
         6,6,6,6,6,6,6,6,6), .Label = c(Group1, Group2, Group3, 
 Group4, Group5, Group6), class = factor),
    Sample = structure(c( 1, 2, 3, 4, 5, 6, 7, 8, 9,
    10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
    20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
    30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 
 40,41,42,43,44,45,46,47,48,49,50,51,52,53,54)
 ))
 , .Names = c(Intensity,
 Group, Sample), row.names =
 c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
 41, 42, 

[R] Add1 w/ coef estimates?

2010-07-06 Thread darckeen

I was wondering if there is anyway to have Add1() display the coefficient
estimates for each candidate predictor along with the F test.  This is for
lm() btw.

Thanks
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Add1-w-coef-estimates-tp2279662p2279662.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] plotmath vector problem; full program enclosed

2010-07-06 Thread Paul Johnson
Here's another example of my plotmath whipping boy, the Normal distribution.

A colleague asks for a Normal plotted above a series of axes that
represent various other distributions (T, etc).

I want to use vectors of equations in plotmath to do this, but have
run into trouble.  Now I've isolated the problem down to a relatively
small piece of working example code (below).  If you would please run
this, I think you will see the problem.  When plotmath meets one
vector of expressions, it converts all but one to math, so in the
figure output i get, in LaTeX speak

b1 $\mu-1.0 \sigma$$\mu$

All values except the first come out correctly.

This happens only when I try to use bquote or substitute to get
variables to fill in where the 1.96, 1.0, and so forth should be.  In
the figure output, you should see a second axis where all of the
symbols are resolved correctly.

As usual, thanks in advance for your help, sorry if I've made an
obvious mistake or overlooked a manual.

### Filename: plotMathProblem.R
### Paul Johnson July 5, 2010
### email me paulj...@ku.edu

  sigma - 10.0
  mu - 4.0

  myx - seq( mu - 3.5*sigma,  mu+ 3.5*sigma, length.out=500)

  myDensity - dnorm(myx,mean=mu,sd=sigma)

  ### xpd needed to allow writing outside strict box of graph
  ### Need big bottom margin to add several x axes
  par(xpd=TRUE, ps=10, mar=c(18,2,2,2))

  plot(myx, myDensity, type=l, xlab=, ylab=Probability Density ,
main=myTitle1, axes=FALSE)
  axis(2, pos= mu - 3.6*sigma)
  axis(1, pos=0)

  lines(c(myx[1],myx[length(myx)]),c(0,0)) ### closes off axes


  addInteriorLine - function(x, m, sd){
for (i in 1:(length(x))){
  lines( c(x[i],x[i]), c(0, dnorm(x[i],m=m,sd=sd)), lty= 14, lwd=.2)
}
  }


  dividers - c(qnorm(0.025), -1, 0, 1, qnorm(0.975))
  addInteriorLine(mu+sigma*dividers, mu,sigma)

  # bquote creates an expression that text plotters can use
  t1 -  bquote( mu== .(mu))
  mtext(bquote( mu == .(mu)), 1, at=mu, line=-1)


  addInteriorLabel - function(pos1, pos2, m, s){
area - abs(100*( pnorm(m+pos1*s,m,s)-pnorm(m+pos2*s, m,s)))
mid - m+0.5*(pos1+pos2)*s
text(mid, 0.5*dnorm(mid,m,s),label=paste(round(area,2),%))
  }


  addInteriorLabel(dividers[1],dividers[2],  mu, sigma)
  addInteriorLabel(dividers[2],dividers[3],  mu, sigma)
  addInteriorLabel(dividers[3],dividers[4],  mu, sigma)
  addInteriorLabel(dividers[4],dividers[5],  mu, sigma)




### Following is problem point: axis will
### end up with correct labels, except for first point,
### where we end up with b1 instead of mu - 1.96*sigma.
   b1 - substitute( mu - d*sigma, list(d=round(dividers[1],2)) )
   b2 - substitute( mu - sigma )
   b3 - substitute( mu )
   b4 - substitute( mu + sigma )
   b5 - substitute( mu + d*sigma, list(d=round(dividers[5],2)) )
##   plot(-20:50,-20:50,type=n,axes=F)
   axis(1, line=4,at=mu+dividers*sigma,
labels=c(expression(b1),b2,b3,b4,b5), padj=-1)




### This gets right result but have to hard code the dividers
  b1 - expression( mu - 1.96*sigma )
  b2 - expression( mu - sigma )
  b3 - expression( mu )
  b4 - expression( mu + sigma )
  b5 - expression( mu + 1.96*sigma )
  axis(1, line=8,at=mu+dividers*sigma, labels=c(b1,b2,b3,b4,b5), padj=-1)




-- 
Paul E. Johnson
Professor, Political Science
1541 Lilac Lane, Room 504
University of Kansas

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

2010-07-06 Thread Tim Smith
Hi,

I was trying to install the package 
'lme4'. Here is the code and the sessionInfo() that I am using:

 install.packages(lme4,dependencies=T)
Warning in 
install.packages(lme4, dependencies = T) :
  argument 'lib' is 
missing: using '/Users/ts/Library/R/2.11/library'
Warning message:
In getDependencies(pkgs, dependencies, available, lib) :
  package 
‘lme4’ is not available



 sessionInfo()
R version 
2.11.1 (2010-05-31) 
x86_64-apple-darwin9.8.0 

attached base packages:
[1] stats graphics  grDevices utils datasets  
methods   base 

loaded via a namespace (and not attached):
[1] grid_2.11.1   nlme_3.1-96   
stats4_2.11.1 tools_2.11.1 
 

Should I be passing in some 
other parameters?

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.


[R] how to calculate summary statistics for each factor

2010-07-06 Thread karena

I have a dataset like the following:

subject   class value
123110
1241 12
125112
223223
224   2 18
225 219
3233 21
324   3  10
325   3   19
326   3   20

how to calculate the summary value for each factor?
thanks

karena
-- 
View this message in context: 
http://r.789695.n4.nabble.com/how-to-calculate-summary-statistics-for-each-factor-tp2279777p2279777.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] information reduction-database management question

2010-07-06 Thread Paul Chatfield

I don't think the approach would change much with text.  You would have
to write a function which picks the 'min' or whatever that means to you
with text and then it should work ok,

 

Paul

 

From: Brad Patrick Schneid [via R]
[mailto:ml-node+2279677-1095983982-120...@n4.nabble.com] 
Sent: 06 July 2010 15:48
To: Paul Chatfield
Subject: Re: information reduction-database management question

 

Thanks Paul, 
I appreciate your time and this is an interesting approach.
Unfortunately, I need it to work for all types of information, including
character data (i.e. text).   
Again.. thanks for your suggestion! 
Brad 

Paul Chatfield wrote:

If you redefine your NAs as below to be detected as some arbitrary large
number, then the code should work through.  Any 5's left in your dataset
can be replaced just as easily by NAs again.  Not elegant, but
effective. 

site - c(s1, s1, s1, s2,s2, s2) 
pref - c(1, 2, 3, 1, 2, 3) 
R1 - c(NA, NA, 1, NA,NA,NA) 
R2 - c(NA, 0, 1, 1, NA, 1) 
R3 - c(NA, 1, 1, NA, 1, 1) 
R4 - c(0, NA, 0, 1, NA, 0) 
R5 - c(NA, 0, 1, NA, 1, 1) 
datum - data.frame(site, pref, R1, R2, R3, R4, R5) 

## For 1 column; 
datum$R1[is.na(datum$R1)==T]-5 
tapply(datum$R1, datum$site, min, na.rm=T) 

## Can loop this over all columns; 

new-matrix(0,5,2) 
for (i in 3:7) 
{datum[,i][is.na(datum[,i])==T]-5 
new[i-2,]-tapply(datum[,i], datum$site, min, na.rm=T)} 

 



View message @
http://r.789695.n4.nabble.com/information-reduction-database-management-
question-tp2278863p2279677.html 
To unsubscribe from Re: information reduction-database management
question, click here
 (link removed) 
NoYXRmaWVsZEByZWFkaW5nLmFjLnVrfDIyNzkzODV8LTE4MjM2NDg5MTM= . 

 


-- 
View this message in context: 
http://r.789695.n4.nabble.com/information-reduction-database-management-question-tp2278863p2279688.html
Sent from the R help mailing list archive at Nabble.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] Differencing with auto.arima and xreg

2010-07-06 Thread Kim H

I am having some issues with differencing using auto.arima when also
specifying an xreg dataframe.

The xreg dataframe contains dummy variables that specify time periods that
had a promotion running.

When I model diff(y) with order (1,0,1), the coefficients for these dummy
variables are very different than when I model y with order=(1,1,1).  I
think when modeling diff(y) the coefficients represent the % change from
week to week (since the response is first logged) and I may be having issues
trying to run a dummy variable over consecutive periods since the promotion
is probably not causing increases in % changes at each time period.  But
when modeling y with order (1,1,1), it allows me to use the dummy variable
over consecutive periods and appears meaningful.

Does anyone have any insight as to why this is happening?
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Differencing-with-auto-arima-and-xreg-tp2279874p2279874.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] package error

2010-07-06 Thread Tim Smith
Hi,

I was trying to install the 'lme4' package on a mac, but keep getting 
installation errors. 

I'm trying to post the actual code that I use, and the output, but it keeps on 
getting bounced off the message board as spam (Moderator approval required). Is 
there something that I should be using/not using in the subject line/body of 
the email?

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] how to save summary(lm) and anova (lm) in format?

2010-07-06 Thread Paul Johnson
There are R packages that can make nice R regression tables in LaTeX
documents. I've used memisc and its good, there is also apsrtable
and the old standby xtable.  Also I use my own function outreg, but
that's just a 'not invented here' attitude.

Your problem is that you need this to go into Word, in which case I
think a reasonable strategy is to create html output in R and then in
Word you can use paste special HTML and it will bring in the html as
a Word table.

I recently made a presentation about this, you might scan down to the
end where I have the html example for the poor-pitiful Word users of
the world:

http://pj.freefaculty.org/SummerCamp2010/regression2.pdf

Look down around slide 75

pj


On Fri, Jul 2, 2010 at 12:34 PM, Yi liuyi.fe...@gmail.com wrote:
 Hi, folks,

 I would like to copy the output of summary(lm) and anova (lm) in R to my
 word file. But the output will be a mess if I just copy after I call summary
 and anova.

 #
 x=rnorm(10)
 y=rnorm(10,mean=3)
 lm=lm(y~x)
 summary(lm)

 Call:
 lm(formula = y ~ x)
 Residuals:
      Min        1Q    Median        3Q       Max
 -1.278567 -0.312017  0.001938  0.297578  1.310113
 Coefficients:
            Estimate Std. Error t value Pr(|t|)
 (Intercept)   2.5221     0.2272  11.103 3.87e-06 ***
 x            -0.5988     0.2731  -2.192   0.0597 .
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 Residual standard error: 0.7181 on 8 degrees of freedom
 Multiple R-squared: 0.3753,     Adjusted R-squared: 0.2972
 F-statistic: 4.807 on 1 and 8 DF,  p-value: 0.0597
 

 How can I get the exact ouput as shown in R but not as the above?


 Thanks.

 Yi

        [[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.





-- 
Paul E. Johnson
Professor, Political Science
1541 Lilac Lane, Room 504
University of Kansas

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

2010-07-06 Thread Alex Foley
Hi,

I was trying to install lme4 package, but got the following errors:

 install.packages(lme4)
Warning in install.packages(lme4) :
  argument 'lib' is missing: using '/Users/xx/Library/R/2.11/library'
Warning message:
In getDependencies(pkgs, dependencies, available, lib) :
  package ‘lme4’ is not available



The session info is:


 sessionInfo()
R version 2.11.1 (2010-05-31) 
x86_64-apple-darwin9.8.0 

locale:
[1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base 

other attached packages:
[1] Matrix_0.999375-39 lattice_0.18-8

loaded via a namespace (and not attached):
[1] grid_2.11.1   nlme_3.1-96   stats4_2.11.1 tools_2.11.1 

Should I be installing this in a different manner?


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] Error message using mi() in mi package

2010-07-06 Thread Peter Ehlers

On 2010-07-06 10:37, Andrew Miles wrote:

Hello!

I get the following message when I run the mi() function from the mi
package.

Error while imputing variable: c3 , model: mi.polr
Error in eval(expr, envir, enclos) : could not find function c14ordered

Here's the situation:

I am running R v. 2.9.2 on Mac OSX v. 10.5.8. I am trying to impute
missing data in a data set that I've trimmed down to 302 variables. I've
created an mi.info() object on the data, and I've updated the type of
variable where necessary so that the mi() imputing function uses the
most appropriate type of models to estimate imputed values. The data
have no collinearlity. I then run the mi function like this:

imp = mi(imp.data, info=info2, n.iter=10)

where imp.data is my data set of 302 variables and info2 is my modified
mi.info object. I get the message as posted above. The message only
occurs when working on a variable I have labeled as
ordered-categorical. But the mi() function processes most variables
labeled as ordered-categorical just fine. In fact, if shrink my data
set (say, to just 5 variables) I can get mi() to process a problematic
variable just fine as well.

I'm not sure what the function c14ordered is that the error message
refers to. My first thought is maybe it is referring to one of my
variables in my data? Variables names in my data follow a basic
letter-number pattern (i.e. a1, a2, etc.), but there is no c14, rather
c14a1, c14a2, etc. So I'm not sure if the variable has anything to do
with the problem, but I thought I'd mention it just in case someone
wiser in this matter than I can see something I cannot.

I cannot post code that reproduces the problem due to the nature of the
code and data involved.

Any help would be appreciated, as I am not sure what is happening, and
can't see why I can sometimes impute a variable labeled as
ordered-categorical and sometimes cannot.

Thanks!

Andrew Miles



This looks suspiciously like a syntax problem.
I would get my text editor to search for 'c14ordered'
in the code. You might have missed some punctuation.

  -Peter Ehlers

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 calculate summary statistics for each factor

2010-07-06 Thread Erik Iverson

?tapply is one way.

karena wrote:

I have a dataset like the following:

subject   class value
123110
1241 12
125112
223223
224   2 18
225 219
3233 21
324   3  10
325   3   19
326   3   20

how to calculate the summary value for each factor?
thanks

karena


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] plotmath vector problem; full program enclosed

2010-07-06 Thread Duncan Murdoch

On 06/07/2010 10:54 AM, Paul Johnson wrote:

Here's another example of my plotmath whipping boy, the Normal distribution.

A colleague asks for a Normal plotted above a series of axes that
represent various other distributions (T, etc).

I want to use vectors of equations in plotmath to do this, but have
run into trouble.  Now I've isolated the problem down to a relatively
small piece of working example code (below).  If you would please run
this, I think you will see the problem.  When plotmath meets one
vector of expressions, it converts all but one to math, so in the
figure output i get, in LaTeX speak

b1 $\mu-1.0 \sigma$$\mu$

All values except the first come out correctly.

This happens only when I try to use bquote or substitute to get
variables to fill in where the 1.96, 1.0, and so forth should be.  In
the figure output, you should see a second axis where all of the
symbols are resolved correctly.

As usual, thanks in advance for your help, sorry if I've made an
obvious mistake or overlooked a manual.

### Filename: plotMathProblem.R
### Paul Johnson July 5, 2010
### email me paulj...@ku.edu

  sigma - 10.0
  mu - 4.0

  myx - seq( mu - 3.5*sigma,  mu+ 3.5*sigma, length.out=500)

  myDensity - dnorm(myx,mean=mu,sd=sigma)

  ### xpd needed to allow writing outside strict box of graph
  ### Need big bottom margin to add several x axes
  par(xpd=TRUE, ps=10, mar=c(18,2,2,2))

  plot(myx, myDensity, type=l, xlab=, ylab=Probability Density ,
main=myTitle1, axes=FALSE)
  axis(2, pos= mu - 3.6*sigma)
  axis(1, pos=0)

  lines(c(myx[1],myx[length(myx)]),c(0,0)) ### closes off axes


  addInteriorLine - function(x, m, sd){
for (i in 1:(length(x))){
  lines( c(x[i],x[i]), c(0, dnorm(x[i],m=m,sd=sd)), lty= 14, lwd=.2)
}
  }


  dividers - c(qnorm(0.025), -1, 0, 1, qnorm(0.975))
  addInteriorLine(mu+sigma*dividers, mu,sigma)

  # bquote creates an expression that text plotters can use
  t1 -  bquote( mu== .(mu))
  mtext(bquote( mu == .(mu)), 1, at=mu, line=-1)


  addInteriorLabel - function(pos1, pos2, m, s){
area - abs(100*( pnorm(m+pos1*s,m,s)-pnorm(m+pos2*s, m,s)))
mid - m+0.5*(pos1+pos2)*s
text(mid, 0.5*dnorm(mid,m,s),label=paste(round(area,2),%))
  }


  addInteriorLabel(dividers[1],dividers[2],  mu, sigma)
  addInteriorLabel(dividers[2],dividers[3],  mu, sigma)
  addInteriorLabel(dividers[3],dividers[4],  mu, sigma)
  addInteriorLabel(dividers[4],dividers[5],  mu, sigma)




### Following is problem point: axis will
### end up with correct labels, except for first point,
### where we end up with b1 instead of mu - 1.96*sigma.
   b1 - substitute( mu - d*sigma, list(d=round(dividers[1],2)) )
   b2 - substitute( mu - sigma )
   b3 - substitute( mu )
   b4 - substitute( mu + sigma )
   b5 - substitute( mu + d*sigma, list(d=round(dividers[5],2)) )
##   plot(-20:50,-20:50,type=n,axes=F)
   axis(1, line=4,at=mu+dividers*sigma,
labels=c(expression(b1),b2,b3,b4,b5), padj=-1)
  


You want as.expression(b1), not expression(b1).  The latter means 
the expression consisting of the symbol b1.  The former means take 
the object stored in b1, and convert it to an expression..


It's not perfect, because you'll end up with mu - -1.96sigma (i.e. two 
minus signs), but it's closer than what you had.


Duncan Murdoch





### This gets right result but have to hard code the dividers
  b1 - expression( mu - 1.96*sigma )
  b2 - expression( mu - sigma )
  b3 - expression( mu )
  b4 - expression( mu + sigma )
  b5 - expression( mu + 1.96*sigma )
  axis(1, line=8,at=mu+dividers*sigma, labels=c(b1,b2,b3,b4,b5), padj=-1)







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


Re: [R] lme4

2010-07-06 Thread Douglas Bates
On OS X you need to install the source package for lme4 as the binary
package fails one of the tests.  We have been unable to reproduce this
failure under other operating systems, which makes it hard to debug.


On Tue, Jul 6, 2010 at 11:09 AM, Alex Foley alex_foley_...@yahoo.com wrote:
 Hi,

 I was trying to install lme4 package, but got the following errors:

 install.packages(lme4)
 Warning in install.packages(lme4) :
  argument 'lib' is missing: using '/Users/xx/Library/R/2.11/library'
 Warning message:
 In getDependencies(pkgs, dependencies, available, lib) :
  package ‘lme4’ is not available



 The session info is:


 sessionInfo()
 R version 2.11.1 (2010-05-31)
 x86_64-apple-darwin9.8.0

 locale:
 [1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8

 attached base packages:
 [1] stats     graphics  grDevices utils     datasets  methods   base

 other attached packages:
 [1] Matrix_0.999375-39 lattice_0.18-8

 loaded via a namespace (and not attached):
 [1] grid_2.11.1   nlme_3.1-96   stats4_2.11.1 tools_2.11.1

 Should I be installing this in a different manner?


 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.



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] plotmath vector problem; full program enclosed

2010-07-06 Thread David Winsemius


On Jul 6, 2010, at 1:41 PM, Duncan Murdoch wrote:


On 06/07/2010 10:54 AM, Paul Johnson wrote:
Here's another example of my plotmath whipping boy, the Normal  
distribution.


A colleague asks for a Normal plotted above a series of axes that
represent various other distributions (T, etc).

I want to use vectors of equations in plotmath to do this, but have
run into trouble.  Now I've isolated the problem down to a relatively
small piece of working example code (below).  If you would please run
this, I think you will see the problem.  When plotmath meets one
vector of expressions, it converts all but one to math, so in the
figure output i get, in LaTeX speak

b1 $\mu-1.0 \sigma$$\mu$

All values except the first come out correctly.

This happens only when I try to use bquote or substitute to get
variables to fill in where the 1.96, 1.0, and so forth should be.  In
the figure output, you should see a second axis where all of the
symbols are resolved correctly.

As usual, thanks in advance for your help, sorry if I've made an
obvious mistake or overlooked a manual.

### Filename: plotMathProblem.R
### Paul Johnson July 5, 2010
### email me paulj...@ku.edu

 sigma - 10.0
 mu - 4.0

 myx - seq( mu - 3.5*sigma,  mu+ 3.5*sigma, length.out=500)

 myDensity - dnorm(myx,mean=mu,sd=sigma)

 ### xpd needed to allow writing outside strict box of graph
 ### Need big bottom margin to add several x axes
 par(xpd=TRUE, ps=10, mar=c(18,2,2,2))

 plot(myx, myDensity, type=l, xlab=, ylab=Probability Density ,
main=myTitle1, axes=FALSE)
 axis(2, pos= mu - 3.6*sigma)
 axis(1, pos=0)

 lines(c(myx[1],myx[length(myx)]),c(0,0)) ### closes off axes


 addInteriorLine - function(x, m, sd){
   for (i in 1:(length(x))){
 lines( c(x[i],x[i]), c(0, dnorm(x[i],m=m,sd=sd)), lty= 14,  
lwd=.2)

   }
 }


 dividers - c(qnorm(0.025), -1, 0, 1, qnorm(0.975))
 addInteriorLine(mu+sigma*dividers, mu,sigma)

 # bquote creates an expression that text plotters can use
 t1 -  bquote( mu== .(mu))
 mtext(bquote( mu == .(mu)), 1, at=mu, line=-1)


 addInteriorLabel - function(pos1, pos2, m, s){
   area - abs(100*( pnorm(m+pos1*s,m,s)-pnorm(m+pos2*s, m,s)))
   mid - m+0.5*(pos1+pos2)*s
   text(mid, 0.5*dnorm(mid,m,s),label=paste(round(area,2),%))
 }


 addInteriorLabel(dividers[1],dividers[2],  mu, sigma)
 addInteriorLabel(dividers[2],dividers[3],  mu, sigma)
 addInteriorLabel(dividers[3],dividers[4],  mu, sigma)
 addInteriorLabel(dividers[4],dividers[5],  mu, sigma)




### Following is problem point: axis will
### end up with correct labels, except for first point,
### where we end up with b1 instead of mu - 1.96*sigma.
  b1 - substitute( mu - d*sigma, list(d=round(dividers[1],2)) )
  b2 - substitute( mu - sigma )
  b3 - substitute( mu )
  b4 - substitute( mu + sigma )
  b5 - substitute( mu + d*sigma, list(d=round(dividers[5],2)) )
##   plot(-20:50,-20:50,type=n,axes=F)
  axis(1, line=4,at=mu+dividers*sigma,
labels=c(expression(b1),b2,b3,b4,b5), padj=-1)



You want as.expression(b1), not expression(b1).  The latter  
means the expression consisting of the symbol b1.  The former  
means take the object stored in b1, and convert it to an  
expression..


It's not perfect, because you'll end up with mu - -1.96sigma (i.e.  
two minus signs), but it's closer than what you had.


Easily addressed in this case with ~ instead of -. The value of  
d provides the minus:


b1 - substitute( mu ~ d*sigma, list(d=round(dividers[1],2)) )



Duncan Murdoch





### This gets right result but have to hard code the dividers
 b1 - expression( mu - 1.96*sigma )
 b2 - expression( mu - sigma )
 b3 - expression( mu )
 b4 - expression( mu + sigma )
 b5 - expression( mu + 1.96*sigma )
 axis(1, line=8,at=mu+dividers*sigma, labels=c(b1,b2,b3,b4,b5),  
padj=-1)








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


David Winsemius, MD
West Hartford, CT

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


Re: [R] Rcmdr installation under Unbuntu installatiion errors

2010-07-06 Thread Vojtěch Zeisek
Hello,
1) Try Rkward, it is much more better,
2) You can install Rcmdr directly via r-cran-rcmdr (generaly, search r-cran-*; 
You can also add special repository, see http://cran.r-
project.org/bin/linux/ubuntu/) and/or
3) if You wish to install normal packages from source (that is IMHO what 
install.packages() does in Linux), You need development packages for R as well 
as kernel header, compiler for Fortran, C and so on. I do not know how are 
those packages named in Ubuntu, but it should be no problem to find it via 
Synaptics. But I'd prefer solution 2. :-)
Best regards,
Vojtěch Zeisek

Dne Út 6. července 2010 18:29:47 John Sorkin napsal(a):
 Unbuntu 10.04
 R 2.10
 
 I am trying to install Rcmdr and receive the following messages:
 
 The downloaded packages are in
 ‘/tmp/RtmpzhjDZG/downloaded_packages’
 Warning messages:
 1: In install.packages(Rcmdr, dependencies = TRUE) :
   installation of package 'Rmpi' had non-zero exit status
 2: In install.packages(Rcmdr, dependencies = TRUE) :
   installation of package 'rpvm' had non-zero exit status
 3: In install.packages(Rcmdr, dependencies = TRUE) :
   installation of package 'rsprng' had non-zero exit status
 4: In install.packages(Rcmdr, dependencies = TRUE) :
   installation of package 'tkrplot' had non-zero exit status
 5: In install.packages(Rcmdr, dependencies = TRUE) :
   installation of package 'rgl' had non-zero exit status
 6: In install.packages(Rcmdr, dependencies = TRUE) :
   installation of package 'RODBC' had non-zero exit status
 7: In install.packages(Rcmdr, dependencies = TRUE) :
   installation of package 'shapes' had non-zero exit status
 8: In install.packages(Rcmdr, dependencies = TRUE) :
   installation of package 'latentnet' had non-zero exit status
 9: In install.packages(Rcmdr, dependencies = TRUE) :
   installation of package 'VIM' had non-zero exit status
 10: In install.packages(Rcmdr, dependencies = TRUE) :
   installation of package 'statnet' had non-zero exit status
 
 Do I need to be concerned? Do I need to take any action?
 
 Thank you,
 John
 
 
 John Sorkin M.D., Ph.D.
 Chief, Biostatistics and Informatics
 Baltimore VA Medical Center GRECC,
 University of Maryland School of Medicine Claude D. Pepper OAIC,
 University of Maryland Clinical Nutrition Research Unit, and
 Baltimore VA Center Stroke of Excellence
 
 University of Maryland School of Medicine
 Division of Gerontology
 Baltimore VA Medical Center
 10 North Greene Street
 GRECC (BT/18/GR)
 Baltimore, MD 21201-1524
 
 (Phone) 410-605-7119
 (Fax) 410-605-7913 (Please call phone number above prior to faxing)
 jsor...@grecc.umaryland.edu
 
 Confidentiality Statement:
 This email message, including any attachments, is for th...{{dropped:6}}
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
  http://www.R-project.org/posting-guide.html and provide commented,
  minimal, self-contained, reproducible code.
-- 
Vojtěch Zeisek

Department of Botany, Faculty of Science, Charles Uni., Prague, CZ
Institute of Botany, Academy of Science, Czech Republic
Community of the openSUSE GNU/Linux

https://www.natur.cuni.cz/faculty-en?set_language=en
http://www.ibot.cas.cz/?p=indexamp;site=en
http://www.opensuse.org/
http://web.natur.cuni.cz/~zeisek/


signature.asc
Description: This is a digitally signed message part.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Pseudo F statistics with index.G1

2010-07-06 Thread Allan Engelhardt
I guess that in R you have to be explicit about what you want to do.  
You can't just drop them, so you'll have to assign them some (other) 
value.  Try which(table(C)==1) to give you the values you need to change 
and then decide what to change them to.  The SAS documentation may tell 
you what it does; otherwise,  I'd suggest you look at the source, but I 
guess that is not an option for you.  Welcome to R.


Allan

On 06/07/10 15:48, Henrik Aldberg wrote:

Thank you Allan,

As i understand it the index.G1 function does not work if one of the 
clusters in the partition only contains one object. Is there a way to 
get around this in R? In SAS the PSF function seems to ignore the 
presence of singleton clusters.



Sincerely

Henrik





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


Re: [R] Error message using mi() in mi package

2010-07-06 Thread Andrew Miles

On Jul 6, 2010, at 1:30 PM, Peter Ehlers wrote:


On 2010-07-06 10:37, Andrew Miles wrote:

Hello!

I get the following message when I run the mi() function from the mi
package.

Error while imputing variable: c3 , model: mi.polr
Error in eval(expr, envir, enclos) : could not find function  
c14ordered


Here's the situation:

I am running R v. 2.9.2 on Mac OSX v. 10.5.8. I am trying to impute
missing data in a data set that I've trimmed down to 302 variables.  
I've
created an mi.info() object on the data, and I've updated the  
type of

variable where necessary so that the mi() imputing function uses the
most appropriate type of models to estimate imputed values. The data
have no collinearlity. I then run the mi function like this:

imp = mi(imp.data, info=info2, n.iter=10)

where imp.data is my data set of 302 variables and info2 is my  
modified

mi.info object. I get the message as posted above. The message only
occurs when working on a variable I have labeled as
ordered-categorical. But the mi() function processes most variables
labeled as ordered-categorical just fine. In fact, if shrink my  
data
set (say, to just 5 variables) I can get mi() to process a  
problematic

variable just fine as well.

I'm not sure what the function c14ordered is that the error message
refers to. My first thought is maybe it is referring to one of my
variables in my data? Variables names in my data follow a basic
letter-number pattern (i.e. a1, a2, etc.), but there is no c14,  
rather

c14a1, c14a2, etc. So I'm not sure if the variable has anything to do
with the problem, but I thought I'd mention it just in case someone
wiser in this matter than I can see something I cannot.

I cannot post code that reproduces the problem due to the nature of  
the

code and data involved.

Any help would be appreciated, as I am not sure what is happening,  
and

can't see why I can sometimes impute a variable labeled as
ordered-categorical and sometimes cannot.

Thanks!

Andrew Miles



This looks suspiciously like a syntax problem.
I would get my text editor to search for 'c14ordered'
in the code. You might have missed some punctuation.

 -Peter Ehlers


A good thought.  I checked my own code (the stuff coding the data and  
using the mi package) and, for good measure, the code of the mi.polr  
function in the mi package, but the string c14ordered never appears.


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


Re: [R] Error message using mi() in mi package

2010-07-06 Thread Erik Iverson





This looks suspiciously like a syntax problem.
I would get my text editor to search for 'c14ordered'
in the code. You might have missed some punctuation.

 -Peter Ehlers


A good thought.  I checked my own code (the stuff coding the data and 
using the mi package) and, for good measure, the code of the mi.polr 
function in the mi package, but the string c14ordered never appears.


And what does

traceback()

tell 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] Help With ANOVA (corrected please ignore last email)

2010-07-06 Thread Joshua Wiley
Hello,

Are you saying that the -4.60517 values are supposed to be treated as
missing?  Unless you set them to NA in R, they will be treated as real
values.  This would make a huge difference.

I can tell you that your formula: aov(Intensity ~ Group, data =
zzzanova) is treating the variable 'Intensity' as the outcome (or
dependent variable) and the variable 'Group' as the predictor (or
independent variable).  But of course, whether you are using an ANOVA
correctly depends on a lot more.  Oh, and yes, Pr(F) is the p-value.

However, I may have found the source of your discrepancy.  Looking at
the data you provided:

 str(zzzanova)
'data.frame':   54 obs. of  3 variables:
 $ Intensity: num  -4.61 -4.61 -4.61 -4.61 -4.61 ...
 $ Group: Factor w/ 6 levels Group1,Group2,..: 1 1 1 1 1 1 1 1 1 2 ...
 $ Sample   : num  1 2 3 4 5 6 7 8 9 10 ...

Notice that 'Group' is a factor with 6 levels.  That means that R is
not treating this variable's scale as interval or ratio variable.  In
fact, it is being treated as nominal.  Looking at the summary of the
ANOVA:

 summary(aov(Intensity ~ Group, data = zzzanova))
Df Sum Sq Mean Sq F value  Pr(F)
Group5  98.85  19.771  2.1469 0.07576 .
Residuals   48 442.03   9.209
---

indicates that the variable 'Group' has 5 degrees of freedom, which
makes sense given that it has 6 levels.  Now look what happens when I
convert the variable 'Group' from a factor to a numeric variable
(using as.numeric() on the variable in the formula).

 summary(aov(Intensity ~ as.numeric(zzzanova$Group), data = zzzanova))
   Df Sum Sq Mean Sq F value   Pr(F)
as.numeric(zzzanova$Group)  1  80.18  80.182  9.0503 0.004042 **
Residuals  52 460.70   8.860

Now Group only has 1 degree of freedom (because it is being treated as
a single predictor now instead of having multiple levels).  Also the
new p-value seems within rounding error of what you expected.


All that said, my guess is that you really do want to treat the
variable 'Group' as a factor.  Treated as a numeric variable what you
are basically testing is the relationship between Group and Intensity
as group increases from 1 to 6.  This is very different from testing
the relationship between 6 different groups and Intensity.

As a disclaimer, that is not really the best explanation of what is
going on, but every other way I tried to say it seemed worse, so if
you are unclear you might consult a statistics text or perhaps others
can give a more erudite explanation.


So to summarize, take a look at the original data for Group that you
pulled into R.  If you do not want it to be a factor you should either
convert it somehow, or read the data in differently.  This may not be
the actual issue, but it was the closest I could come to replicating
your p-value.

Best regards,

Josh


On Tue, Jul 6, 2010 at 10:09 AM, Amit Patel amitrh...@yahoo.co.uk wrote:
 Hi Joshua

 Thanks very much for your help. I will take your advice and work on my 
 problem. I am getting my expected p-values from GeneSpring software. I 
 thought the problem lied in the treatment of the NA (-4.60517) values and 
 maybe how they are treated in the analysis (i.e ignored or -4.60517). Im not 
 sure what your ANOVA background is but I just wanted to check that I'm using 
 ANOVA correctly. Im using the assumption that the Pr(F)1 value is the actual 
 p-value.

 Thanks again for your help



 - Original Message 
 From: Joshua Wiley jwiley.ps...@gmail.com
 To: Amit Patel amitrh...@yahoo.co.uk
 Cc: r-help@r-project.org
 Sent: Tue, 6 July, 2010 16:46:30
 Subject: Re: [R] Help With ANOVA (corrected please ignore last email)

 Hi Amit,

 When I copy in your data and run

 aov(Intensity ~ Group, data = zzzanova)

 I get neither the p-value you showed nor the one you expected.  My
 suggestions at things to look at would be

 1) Where/How did you get the expected p-value?  Another statistics
 program (e.g., SPSS or SAS)?  It helps to know where the expected
 value came from.  If it was from another program, I would also
 recommend searching the R-help archives (for instance for SPSS and
 ANOVA)

 2) Perhaps something is happening with the unbalanced design (e.g.,
 aov() in R handles it differently than you expected)?

 3) You only mention that the p-values do not match.  What about other
 aspects?  Are the df the same?  I would try to be certain that the
 data in R is the same as what was used to calculate the expected
 p-value.  summary(zzzanova) will return some nice summary statistics
 for each column of your dataframe.

 Cheers,


 Josh



 On Tue, Jul 6, 2010 at 6:12 AM, Amit Patel amitrh...@yahoo.co.uk wrote:
 Sorry i had a misprint in the appendix code in the last email


 Hi I needed some help with ANOVA

 I have a problem with My ANOVA
 analysis. I have a dataset with a known ANOVA p-value, however I can
 not seem to re-create it in R.

 I have created a list (zzzanova) which contains
 1)Intensity Values
 2)Group Number 

Re: [R] Error message using mi() in mi package

2010-07-06 Thread Andrew Miles



On Jul 6, 2010, at 2:15 PM, Erik Iverson wrote:






This looks suspiciously like a syntax problem.
I would get my text editor to search for 'c14ordered'
in the code. You might have missed some punctuation.

-Peter Ehlers
A good thought.  I checked my own code (the stuff coding the data  
and using the mi package) and, for good measure, the code of the  
mi.polr function in the mi package, but the string c14ordered  
never appears.


And what does

traceback()

tell you?



Here is the original error message:

Error while imputing variable: a8 , model: mi.polr
 Error in eval(expr, envir, enclos) : could not find function  
c14ordered


I confess I'm not good at reading traceback() output, but here is what  
it says:


18: eval(expr, envir, enclos)
17: eval(predvars, data, env)
16: model.frame.default(formula = form, data = data,  
drop.unused.levels = FALSE)

15: model.frame(formula = form, data = data, drop.unused.levels = FALSE)
14: model.frame(formula = form, data = data, drop.unused.levels = FALSE)
13: eval(expr, envir, enclos)
12: eval(expr, p)
11: eval.parent(m)
10: bayespolr(formula = form, data = data, start = 0, method =  
c(logistic),

drop.unused.levels = FALSE, n.iter = n.iter)
9: mi.polr(formula = a8 ~ rural_now + a3 + a7 . . . [more arguments  
here including the rest of the variables in the data set and a list of  
all the actual data], start = NULL,

   n.iter = 100, draw.from.beta = TRUE)
8: do.call(model.type, args = c(list(formula = info[[CurrentVar]] 
$imp.formula,

   data = dat, start = if (!is.null(start.val[[i]][[jj]])) {
   start.val[[i]][[jj]]
   } else {
   NULL
   }), info[[CurrentVar]]$params))
7: eval(expr, envir, enclos)
6: eval(substitute(expr), data, enclos = parent.frame())
5: with.default(data = dat, do.call(model.type, args = c(list(formula  
= info[[CurrentVar]]$imp.formula,

   data = dat, start = if (!is.null(start.val[[i]][[jj]])) {
   start.val[[i]][[jj]]
   } else {
   NULL
   }), info[[CurrentVar]]$params)))
4: with(data = dat, do.call(model.type, args = c(list(formula =  
info[[CurrentVar]]$imp.formula,

   data = dat, start = if (!is.null(start.val[[i]][[jj]])) {
   start.val[[i]][[jj]]
   } else {
   NULL
   }), info[[CurrentVar]]$params)))
3: .local(object, ...)
2: mi(imp.data, info = info2, n.iter = 10)
1: mi(imp.data, info = info2, n.iter = 10)

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

2010-07-06 Thread Marino Taussig De Bodonia, Agnese
Hello,

I am currently analyzing responses to questionnaires about general attitudes. I 
have performed a PCA on my data, and have retained two Principal Components. 
Now I would like to use the scores of both the principal comonents in a 
multiple regression. I would like to know if it makes sense to use the scores 
of one principal component to explain the variance in the scores of another 
principal component:

lm(scores of principal component 1~scores of principal component 2+ age, 
gender, etc..)

Please could you let me know if this is statistically sound?

Thank you in advance for you time,

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

2010-07-06 Thread Xiaoxi Gao

Hello R users,

I have two quick questions while using lpSolve package for linear 
programming. (1) the result contains both characters and numbers, e.g., 
Success: the objective function is 40.5, but I only need the number, can I only 
store the number? (2) How to set boundaries for variables? e.g., all variable 
are positive.

Thanks a lot!

Xiaoxi 

  
_
Hotmail has tools for the New Busy. Search, chat and e-mail from your inbox.

N:WL:en-US:WM_HMP:042010_1
[[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] Get the indices of non-zero entries of a sparse matrix in R

2010-07-06 Thread G FANG
Hi,

I am trying to get the indices of non-zero entries of a sparse matrix in R

 sr d
1 1089 3772 1
2 1109  190 1
3 1109 2460 1
4 1109 3071 1
5 1109 3618 1
6 1109   38 1

I found that the following can create a sparse matrix,

library(Matrix)
Y - sparseMatrix(s,r,x=d)

but have not idea and did not find online how to convert a sparse
matrix back to three columns efficiently, i.e. get the indices of
non-zero entries of a sparse matrix.

In matlab, the find function can do this efficiently, I am wondering
is there a similar one in R?

Please advice. Thanks!

--gang

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] plotmath vector problem; full program enclosed

2010-07-06 Thread Allan Engelhardt



On 06/07/10 18:51, David Winsemius wrote:


Easily addressed in this case with ~ instead of -. The value of 
d provides the minus:


b1 - substitute( mu ~ d*sigma, list(d=round(dividers[1],2)) )


Neat trick!  But it gives a slightly different minus sign in the 
display, so perhaps simply


b1- substitute( mu - d*sigma, list(d=-round(dividers[1],2)) )


instead?

Allan

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

2010-07-06 Thread David Winsemius


On Jul 6, 2010, at 3:32 PM, Xiaoxi Gao wrote:



Hello R users,

I have two quick questions while using lpSolve package for linear  
programming. (1) the result contains both characters and numbers,  
e.g., Success: the objective function is 40.5, but I only need the  
number, can I only store the number?


?lp
?lp.object   # or just use str()

(2) How to set boundaries for variables? e.g., all variable are  
positive.


It appears you need to take time to re-read the help page and work  
through the examples.





--

David Winsemius, MD
West Hartford, CT

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


Re: [R] plotmath vector problem; full program enclosed

2010-07-06 Thread Charles C. Berry

On Tue, 6 Jul 2010, David Winsemius wrote:



On Jul 6, 2010, at 1:41 PM, Duncan Murdoch wrote:


On 06/07/2010 10:54 AM, Paul Johnson wrote:
 Here's another example of my plotmath whipping boy, the Normal 
 distribution.
 
 A colleague asks for a Normal plotted above a series of axes that

 represent various other distributions (T, etc).
 
 I want to use vectors of equations in plotmath to do this, but have

 run into trouble.  Now I've isolated the problem down to a relatively
 small piece of working example code (below).  If you would please run
 this, I think you will see the problem.  When plotmath meets one
 vector of expressions, it converts all but one to math, so in the
 figure output i get, in LaTeX speak
 
 b1 $\mu-1.0 \sigma$$\mu$
 
 All values except the first come out correctly.
 
 This happens only when I try to use bquote or substitute to get

 variables to fill in where the 1.96, 1.0, and so forth should be.  In
 the figure output, you should see a second axis where all of the
 symbols are resolved correctly.
 
 As usual, thanks in advance for your help, sorry if I've made an

 obvious mistake or overlooked a manual.
 
 ### Filename: plotMathProblem.R

 ### Paul Johnson July 5, 2010
 ### email me paulj...@ku.edu
 
  sigma - 10.0

  mu - 4.0
 
  myx - seq( mu - 3.5*sigma,  mu+ 3.5*sigma, length.out=500)
 
  myDensity - dnorm(myx,mean=mu,sd=sigma)
 
  ### xpd needed to allow writing outside strict box of graph

  ### Need big bottom margin to add several x axes
  par(xpd=TRUE, ps=10, mar=c(18,2,2,2))
 
 plot(myx, myDensity, type=l, xlab=, ylab=Probability Density ,

 main=myTitle1, axes=FALSE)
  axis(2, pos= mu - 3.6*sigma)
  axis(1, pos=0)
 
  lines(c(myx[1],myx[length(myx)]),c(0,0)) ### closes off axes
 
 
  addInteriorLine - function(x, m, sd){

for (i in 1:(length(x))){
  lines( c(x[i],x[i]), c(0, dnorm(x[i],m=m,sd=sd)), lty= 14, lwd=.2)
}

}
 
 
  dividers - c(qnorm(0.025), -1, 0, 1, qnorm(0.975))

  addInteriorLine(mu+sigma*dividers, mu,sigma)
 
  # bquote creates an expression that text plotters can use

  t1 -  bquote( mu== .(mu))
  mtext(bquote( mu == .(mu)), 1, at=mu, line=-1)
 
 
  addInteriorLabel - function(pos1, pos2, m, s){

area - abs(100*( pnorm(m+pos1*s,m,s)-pnorm(m+pos2*s, m,s)))
mid - m+0.5*(pos1+pos2)*s
text(mid, 0.5*dnorm(mid,m,s),label=paste(round(area,2),%))

}
 
 
  addInteriorLabel(dividers[1],dividers[2],  mu, sigma)

  addInteriorLabel(dividers[2],dividers[3],  mu, sigma)
  addInteriorLabel(dividers[3],dividers[4],  mu, sigma)
  addInteriorLabel(dividers[4],dividers[5],  mu, sigma)
 
 
 
 
 ### Following is problem point: axis will

 ### end up with correct labels, except for first point,
 ### where we end up with b1 instead of mu - 1.96*sigma.
   b1 - substitute( mu - d*sigma, list(d=round(dividers[1],2)) )
   b2 - substitute( mu - sigma )
   b3 - substitute( mu )
   b4 - substitute( mu + sigma )
   b5 - substitute( mu + d*sigma, list(d=round(dividers[5],2)) )
 ##   plot(-20:50,-20:50,type=n,axes=F)
   axis(1, line=4,at=mu+dividers*sigma,
 labels=c(expression(b1),b2,b3,b4,b5), padj=-1)
 

You want as.expression(b1), not expression(b1).  The latter means the 
expression consisting of the symbol b1.  The former means take the object 
stored in b1, and convert it to an expression..


It's not perfect, because you'll end up with mu - -1.96sigma (i.e. two 
minus signs), but it's closer than what you had.


Easily addressed in this case with ~ instead of -. The value of d 
provides the minus:


b1 - substitute( mu ~ d*sigma, list(d=round(dividers[1],2)) )



But if d = 0, there will be no sign so maybe use

b1 - substitute( mu ~ d*sigma, list(d=sprintf(%+.2f, dividers[1]))) 
b5 - substitute( mu ~ d*sigma, list(d=sprintf(%+.2f, dividers[5])))


etc.

Chuck



Duncan Murdoch

 
 
 
 ### This gets right result but have to hard code the dividers

  b1 - expression( mu - 1.96*sigma )
  b2 - expression( mu - sigma )
  b3 - expression( mu )
  b4 - expression( mu + sigma )
  b5 - expression( mu + 1.96*sigma )
  axis(1, line=8,at=mu+dividers*sigma, labels=c(b1,b2,b3,b4,b5), padj=-1)
 
 
 
 
 


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

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


David Winsemius, MD
West Hartford, CT

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



Charles C. Berry(858) 534-2098
Dept of Family/Preventive Medicine
E mailto:cbe...@tajo.ucsd.edu   UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901


[R] How to plot confidence bands for nls

2010-07-06 Thread Claudia Penaloza
I adjusted an exponential regression to the following data and wish to plot
confidence bands as well. Is this possible?

Any help greatly appreciated.
Claudia

x - c(1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,
2002,2003,2004,2005,2006,2007,2008,2009)
y - c(987,937,810,749,1087,807,1050,1294,1455,1022,927,403,698,1191,1078,
1176,1125,936,1263,647,868)
dat - data.frame(x,y)

f - function(x,a,b){a*exp(b*x)}
fems - nls(y ~ f(x,a,b), data = dat, start = c(a=100, b=0))

plot(y~x, xlim=c(1989,2010), ylab=Nesting Females, xlab=Year)

curve(f(x, a=co[1], b=co[2]), add=T)

[[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.


  1   2   >