Re: [R] file copy to password protected network drive

2013-02-15 Thread Jeff Newmiller
I have no idea what TACC is, or what your OS is, or what file networking 
scheme your system is using, and these issues are all outside the topic area 
for this list. You should go talk to your network admin or local help desk 
about how to accomplish this task at the command line, and then if it involves 
anything more than ordinary file system access then you should Google for how 
to invoke shell commands on your OS (shell/system/system2).
---
Jeff NewmillerThe .   .  Go Live...
DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live Go...
  Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
--- 
Sent from my phone. Please excuse my brevity.

Kumar Mainali kpmain...@gmail.com wrote:

I am trying to copy files to a password protected drive which is
Ranch at
TACC from another network drive. I am logged in to the source drive and
can
run R there. The following code does not even find the destination
folder.

file.copy(sourcedrive/file.tar, 
usern...@ranch.tacc.utexas.edu/uniqueID/destinationfolder/file.tar,
overwrite = FALSE)

Thanks in advance,
Kumar

On Thu, Feb 14, 2013 at 4:41 AM, e-letter inp...@gmail.com wrote:

 Readers,

 For this data set:

 testvalues-c(10,20,30,40)

 How to amend the plot instruction:

 plot(testvalues,ann=FALSE,type='l',yaxt='n',xaxt='n')

 so that x axis ticks labels can be added to existing graph with
 arbitrary value such as 0,100,200,300)?

 Thanks in advance.

 --
 r2151

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] plot custom x axis ticks values

2013-02-15 Thread e-letter
On 14/02/2013, MacQueen, Don macque...@llnl.gov wrote:
 plot(testvalues,ann=FALSE,type='l',yaxt='n',xaxt='n')
 par()$usr
 [1]  0.88  4.12  8.80 41.20

 The x axis range is from 0.88 to 4.12, so tick labels at 0, 100, 200, 300
 makes no sense.


True per se, but the purpose of the tick labels is to indicate
approximate correlation with another data set.

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


Re: [R] plot custom x axis ticks values

2013-02-15 Thread e-letter
On 15/02/2013, Jim Lemon j...@bitwrit.com.au wrote:
 On 02/14/2013 09:41 PM, e-letter wrote:
 Readers,

 For this data set:

 testvalues-c(10,20,30,40)

 How to amend the plot instruction:

 plot(testvalues,ann=FALSE,type='l',yaxt='n',xaxt='n')

 so that x axis ticks labels can be added to existing graph with
 arbitrary value such as 0,100,200,300)?

 Hi r2151,
 If you want the labels to fit on the axis you will have to include this
 information in the call to plot:

 plot(testvalues,ann=FALSE,type='l',yaxt='n',xaxt='n',ylim=c(0,300))
 axis(2,at=c(0,100,200,300))


plot(testvalues,ann=FALSE,type='l',yaxt='n',xaxt='n',xlim=c(0,300))
axis(1,at=c(0,100,200,300))

The use of 'xlim' makes the graph unacceptable. Instead, it would be
better to plot:

plot(testvalues,ann=FALSE,type='l',yaxt='n')

The x-axis shows tick marks '1.0', '1.5', ...

It is required to replace these values with custom values (e.g. 10, 20... etc.)

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


Re: [R] Troubleshooting underidentification issues in structural equation modelling (SEM)

2013-02-15 Thread Ruijie
Thanks Prof Fox for your guidance. My purpose in fitting this model is to
contrast it with another model that I am proposing which I believe will be
a better fit.

On the point of some of the items being close to invariant, I had a close
look at my data and indeed that is the case I am aware of it. However, I am
not sure what to do with these items. Do I remove them? If I do, what
threshold of variance do I set for removal? How do I decide on that
threshold?

I've combed a number of textbooks for answers but sadly have not found
much. Hope you could offer some advice, thanks!

Regards,
Ruijie (RJ)


He who has a why can endure any how.

~ Friedrich Nietzsche


On 10 February 2013 00:38, John Fox j...@mcmaster.ca wrote:

 Dear Ruijie,

 Your model is underidentified by virtue of two of the factors having only
 one observed indicator each. No SEM software can magically estimate this
 model as it stands. Beyond that, I won't comment on the wisdom of what
 you're doing, such as computing covariances between ordinal variables --
 but
 see what I discovered below.

 Removing these two variables and the associated factors produces the
 following model:

 - snip 

  model - cfa(reference.indicators=FALSE)
 1: F01: I01, I02, I03
 2: F02: I04, I05, I06, I07, I08, I09, I10, I11, I12, I13
 3: F03: I14, I15, I16, I17, I18, I19, I20, I21, I22, I23, I24, I25, I26
 4: F04: I27, I28, I29, I30, I31, I32, I33, I34
 5: F05: I35, I36, I37, I38, I39, I40, I41, I42, I43
 6: F07: I46, I47, I48, I49, I50, I51
 7: F08: I54, I55, I56, I57, I58, I59, I60, I61, I62, I63, I64
 8: F09: I65, I66, I67
 9: F11: I69, I70, I71
 10:
 Read 9 items
 NOTE: adding 66 variances to the model
 
  cfa.output - sem(model, cov.mat, N = 900)

 - snip 

 sem() ran out of iterations, but the summary output is revealing:

 - snip 

  summary(cfa.output)

  Model Chisquare =  5677.1   Df =  2043 Pr(Chisq) = 0
  AIC =  6013.1
  BIC =  -8220.193

  Normalized Residuals
Min. 1st Qu.  MedianMean 3rd Qu.Max.
 -3.9910 -0.5887 -0.1486  0.2588  0.8092 17.2900

  R-square for Endogenous Variables
 I01 I02 I03 I04 I05 I06 I07 I08 I09
 I10
  0.0953  0.1263  0.  0.1131  0.4039  0.2519  0.1168  0.0468  0.0005
 0.0059
 I11 I12 I13 I14 I15 I16 I17 I18 I19
 I20
  0.0479  0.0228  0.1150  0.2813  0.0001  0.0388  0.2106  0.0001  0.0913
 0.0063
 I21 I22 I23 I24 I25 I26 I27 I28 I29
 I30
  0.0041  0.0077  0.0022  0.  0.0299  0.0067  0.0019  0.0011  0.0010
 0.
 I31 I32 I33 I34 I35 I36 I37 I38 I39
 I40
  0.0005  0.0117  0.0270  0.0001  0.0084  0.0001  0.0256  0.4969  0.0613
 0.0515
 I41 I42 I43 I46 I47 I48 I49 I50 I51
 I54
  0.0005  0.0052  0.0307  0.0003  0.1131  0.0014  0.  0.1276  0.9728
 0.0520
 I55 I56 I57 I58 I59 I60 I61 I62 I63
 I64
  0.2930  0.0127  0.0543  0.0500  0.0378  0.0001  0.3048  0.0002  0.0304
 0.0001
 I65 I66 I67 I69 I70 I71
 56.7264  0.  0.0002  0.2220  0.2342  0.2240

  Parameter Estimates
  Estimate  Std Errorz value  Pr(|z|)

 lam[I01:F01]  3.023074e-02 5.133785e-03  5.888586224  3.895133e-09 I01 ---
 F01
 lam[I02:F01]  3.283192e-02 5.291069e-03  6.205157975  5.464199e-10 I02 ---
 F01
 lam[I03:F01]  1.123398e-04 2.695713e-03  0.041673509  9.667590e-01 I03 ---
 F01
 lam[I04:F02]  1.365329e-01 1.555023e-02  8.780124358  1.632940e-18 I04 ---
 F02
 lam[I05:F02]  9.525580e-02 5.517838e-03 17.263245517  8.896692e-67 I05 ---
 F02
 lam[I06:F02]  1.720147e-01 1.277593e-02 13.463962882  2.548717e-41 I06 ---
 F02
 lam[I07:F02]  3.164280e-02 3.543421e-03  8.930015663  4.259485e-19 I07 ---
 F02
 lam[I08:F02]  5.685988e-02 1.021854e-02  5.564386503  2.630763e-08 I08 ---
 F02
 lam[I09:F02]  1.234516e-03 2.228298e-03  0.554017268  5.795670e-01 I09 ---
 F02
 lam[I10:F02]  1.656005e-02 8.458411e-03  1.957820181  5.025112e-02 I10 ---
 F02
 lam[I11:F02]  8.785114e-02 1.560646e-02  5.629151062  1.810987e-08 I11 ---
 F02
 lam[I12:F02]  3.022114e-02 7.815459e-03  3.866842129  1.102537e-04 I12 ---
 F02
 lam[I13:F02]  5.075487e-02 5.732307e-03  8.854177302  8.430329e-19 I13 ---
 F02
 lam[I14:F03]  2.587670e-01 2.308125e-02 11.211137448  3.595430e-29 I14 ---
 F03
 lam[I15:F03] -2.999816e-04 1.469667e-03 -0.204115351  8.382634e-01 I15 ---
 F03
 lam[I16:F03]  2.314973e-02 5.256310e-03  4.404179628  1.061849e-05 I16 ---
 F03
 lam[I17:F03]  9.333201e-02 9.301123e-03 10.034488472  1.075152e-23 I17 ---
 F03
 lam[I18:F03] -3.389770e-04 1.469665e-03 -0.230649144  8.175874e-01 I18 ---
 F03
 lam[I19:F03]  6.783532e-02 1.005099e-02  6.749117110  1.487475e-11 I19 ---
 F03
 lam[I20:F03]  3.916003e-02 2.208166e-02  1.773418523  7.615938e-02 I20 ---
 F03
 lam[I21:F03]  7.260062e-03 5.059696e-03  1.434881038  1.513210e-01 I21 ---
 F03
 

Re: [R] plot custom x axis ticks values

2013-02-15 Thread Duncan Murdoch

On 13-02-15 3:28 AM, e-letter wrote:

On 15/02/2013, Jim Lemon j...@bitwrit.com.au wrote:

On 02/14/2013 09:41 PM, e-letter wrote:

Readers,

For this data set:

testvalues-c(10,20,30,40)

How to amend the plot instruction:

plot(testvalues,ann=FALSE,type='l',yaxt='n',xaxt='n')

so that x axis ticks labels can be added to existing graph with
arbitrary value such as 0,100,200,300)?


Hi r2151,
If you want the labels to fit on the axis you will have to include this
information in the call to plot:

plot(testvalues,ann=FALSE,type='l',yaxt='n',xaxt='n',ylim=c(0,300))
axis(2,at=c(0,100,200,300))



plot(testvalues,ann=FALSE,type='l',yaxt='n',xaxt='n',xlim=c(0,300))
axis(1,at=c(0,100,200,300))

The use of 'xlim' makes the graph unacceptable. Instead, it would be
better to plot:

plot(testvalues,ann=FALSE,type='l',yaxt='n')

The x-axis shows tick marks '1.0', '1.5', ...

It is required to replace these values with custom values (e.g. 10, 20... etc.)


Why not rescale the values before plotting, and use the automatic ticks? 
 You can lie about the user scale, but it doesn't always give a helpful 
plot, e.g.


axis(1, at=1:4, labels=c(100, 300, 400, 200))

is kind of hard to interpret.

Duncan Murdoch

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


Re: [R] An extended Hodgkin-Huxley model that doesn't want to work.

2013-02-15 Thread Jannetta Steyn
Sorry that was supposed to be Berend and not Ben


On 15 February 2013 09:28, Jannetta Steyn janne...@henning.org wrote:

 Hi Ben

 Thank you so much for spending the time to look at my code. I really
 appreciate it.

 The unnecessary parameters were artefacts from me working on the model and
 putting things in and taking things out.

 The model isn't quite producing what I expect yet, but it is definitely
 starting to look more interesting :-)

 I am now getting this error:

  plot(ode)
 Error in is.null(times) : 'times' is missing

 What on earth does it mean because I cant see any 'times' being missing.
 To me it seems that all is there?!?

 Kind Regards
 Jannetta


 On 14 February 2013 17:41, Berend Hasselman b...@xs4all.nl wrote:


 Forgot Reply to All.

 On 13-02-2013, at 23:30, Jannetta Steyn janne...@henning.org wrote:

  Hi All
 
  I have been struggling with this model for some time now and I just
 can't
  get it to work correctly. The messages I get when running the code is:
 
  DLSODA-  Warning..Internal T (=R1) and H (=R2) are
  such that in the machine, T + H = T on the next step
 (H = step size). Solver will continue anyway.
  In above message, R =
  [1] 0 0
  DINTDY-  T (=R1) illegal
  In above message, R =
  [1] 0.1
  T not in interval TCUR - HU (= R1) to TCUR (=R2)
  In above message, R =
  [1] 0 0
  DINTDY-  T (=R1) illegal
  In above message, R =
  [1] 0.2
  T not in interval TCUR - HU (= R1) to TCUR (=R2)
  In above message, R =
  [1] 0 0
  DLSODA-  Trouble in DINTDY.  ITASK = I1, TOUT = R1
  In above message, I =
  [1] 1
  In above message, R =
  [1] 0.2
  Error in lsoda(y, times, func, parms, ...) :
  illegal input detected before taking any integration steps - see written
  message
 
 
 
  I'll first paste the formulae and then I'll paste my code. If anyone can
  spot something wrong with my implementation it would really make my day.
 
  (1)
  dV/dt = (I_ext - I_int-I_coup)/C
  I_ext = injected current
  I_int = Sum of all ion currents
  I_coup = coupling current (but we're not using it here )
 
  (2)
  I_i = g_i * m_i^pi * h_i^pi(V-E)
  i identifies the ion, thus I_K would be Potassium current.
 
  (3)
  dm/dt = (m_inf*V - m)/tau_m
 
  (4)
  dh/dt = (h_inf*V-h)/tau_h
 
  (5)
  The Nernst equation is used to calculate reversal potential for Ca:
  Eca = 12.2396 * log(13000/Ca2+)
 
  (6)
  d[Ca_2+]/dt = (F*I_Ca - [Ca2+] + C0)/Tau_Ca
 
 
  tau_m, tau_h, m_inf and h_inf are all calculated according to formulae
  provided in a paper. In my code these are calculated for the different
  channels into the following variables:
 
  CaTminf, CaThinf, CaTtaum, CaTtauh, CaSminf, CaStaum, Napminf, Naphinf,
  taumna, tauhna, hminf, htaum, Kminf and Ktaum
 
  The E (reversal potential) values for all the channels are given, except
  for CaT and CaS which uses Eca as calculated in (5).
 
  Current for Ca is calculated by summing the CaT and CaS currents, hence
  CaI = gCaT*CaTm^3*CaTh*(v-Eca(v)) + gCaS*CaSm^3(v-ECa(v)
 
 
  Here is the code:
 
  library(simecol)
  ## Hodkin-Huxley model
  HH_soma -  function(time, init, parms) {
  with(as.list(c(init, parms)),{
 
# Na only used in Axon
#Naminf -1/(1+exp(-(v+24.7)/5.29));
#Nataum - function(v) 1.32 - (1.26/(1+exp(-(v+120)/25)));
#Nahinf -1/(1+exp((v+489)/5.18));
#Natauh -(0.67/(1+exp(-(v+62.9)/10))) *
 (1.5+(1/(1+exp((v+34.9)/36;
 
#PD
# mca10
CaTminf - function(v) 1/(1+exp(-(v+25)/7.2));
# hca10
CaThinf - function(v) 1/(1+exp(v+36)/7);
# taumca1
CaTtaum - function(v) 55- (49.5/(1+exp(-v+58)/17));
# tauhca1
CaTtauh - function(v) 350 - (300/(1+exp(-v+50)/16.9));
 
#mca20
CaSminf - function(v) 1/(1+exp(-(v+22)/8.5));
#taumca2
CaStaum - function(v) 16-(13.1/(1+exp(-(v+25.1)/26.4)));
 
 
# mna0
Napminf - function(v) 1/(1+exp(-(v+26.8)/8.2));
# hna0
Naphinf - function(v) 1/(1+exp(-(v+48.5)/5.18));
 
taumna - function(v) 19.8-(10.7/(1+exp(-(v+26.5)/8.6)));
tauhna - function(v) 666-(379/(1+exp(-(v+33.6)/11.7)));
 
# mh0
hminf - function(v)  1/(1+exp(v+70)/6);
# taumh
htaum - function(v)  272+(1499/(1+exp(-(v+42.2)/8.73)));
 
Kminf - function(v)  1/(1+exp(-(v+14.2)/11.8));
Ktaum - function(v)  7.2-(6.4/(1+exp(-(v+28.3)/19.2)));
 
# Reversal potential of intracellular calcium concentration
# Nernst Equation using  extracellular concentration of Ca = 13mM
# eca
ECa - function(Ca2) 12.2396*log(13000/(Ca2));
#ECa - function(CaI) 12.2396*log(13000/(CaI));
 
 
#Sum of all the Ca
# function(v) CaTminf(v) + CaSminf(v);
CaI - gCaT*CaTm^3*CaTh*(v-ECa(CaI)) + gCaS*CaSm^3*(v-ECa(CaI))
 
#AB
#dCa2 - (((-F*Caminf(v))-Caminf(v) + C0)/TauCa)
dCa2 - (((-F*CaI) - Ca2 + C0)/TauCa)
 
# mk20
KCaminf - function(v, Ca2) (Ca2/(Ca2+30))*(1/(1+exp(-(v+51)/8)));
# taumk
KCataum - function(v) 90.3 - ((75.09/(1+exp(-(v+46)/22.7;
 
#AB
Aminf - function(v) 

Re: [R] An extended Hodgkin-Huxley model that doesn't want to work.

2013-02-15 Thread Berend Hasselman

On 15-02-2013, at 10:28, Jannetta Steyn janne...@henning.org wrote:

 Hi Ben
 
 Thank you so much for spending the time to look at my code. I really
 appreciate it.
 
 The unnecessary parameters were artefacts from me working on the model and
 putting things in and taking things out.
 
 The model isn't quite producing what I expect yet, but it is definitely
 starting to look more interesting :-)
 
 I am now getting this error:
 
 plot(ode)
 Error in is.null(times) : 'times' is missing
 
 What on earth does it mean because I cant see any 'times' being missing. To
 me it seems that all is there?!?
 

Difficult to tell as you haven't shown how ode was defined.

I just did this:

out-ode(y=init, times=times, func=HH_soma, parms=parms)
plot(out)

and got the plots of all state variables.

Berend

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


Re: [R] An extended Hodgkin-Huxley model that doesn't want to work.

2013-02-15 Thread Jannetta Steyn
Hi Berend

This is the code. Pretty much just changed to what you suggested which is
CaI=1, removing the unnecessary variables and using deSolve:

rm(list = ls())
library(deSolve)
## Hodkin-Huxley model
HH_soma -  function(times, init, parms) {
  with(as.list(c(init, parms)),{

# Na only used in Axon
#Naminf -1/(1+exp(-(v+24.7)/5.29));
#Nataum - function(v) 1.32 - (1.26/(1+exp(-(v+120)/25)));
#Nahinf -1/(1+exp((v+489)/5.18));
#Natauh -(0.67/(1+exp(-(v+62.9)/10))) * (1.5+(1/(1+exp((v+34.9)/36;

#PD
# mca10
CaTminf - function(v) 1/(1+exp(-(v+25)/7.2));
# hca10
CaThinf - function(v) 1/(1+exp(v+36)/7);
# taumca1
CaTtaum - function(v) 55- (49.5/(1+exp(-v+58)/17));
# tauhca1
CaTtauh - function(v) 350 - (300/(1+exp(-v+50)/16.9));

#mca20
CaSminf - function(v) 1/(1+exp(-(v+22)/8.5));
#taumca2
CaStaum - function(v) 16-(13.1/(1+exp(-(v+25.1)/26.4)));


# mna0
Napminf - function(v) 1/(1+exp(-(v+26.8)/8.2));
# hna0
Naphinf - function(v) 1/(1+exp(-(v+48.5)/5.18));

taumna - function(v) 19.8-(10.7/(1+exp(-(v+26.5)/8.6)));
tauhna - function(v) 666-(379/(1+exp(-(v+33.6)/11.7)));

# mh0
hminf - function(v)  1/(1+exp(v+70)/6);
# taumh
htaum - function(v)  272+(1499/(1+exp(-(v+42.2)/8.73)));

Kminf - function(v)  1/(1+exp(-(v+14.2)/11.8));
Ktaum - function(v)  7.2-(6.4/(1+exp(-(v+28.3)/19.2)));

# Reversal potential of intracellular calcium concentration
# Nernst Equation using  extracellular concentration of Ca = 13mM
# eca
ECa - function(Ca2) 12.2396*log(13000/(Ca2));
#ECa - function(CaI) 12.2396*log(13000/(CaI));


#Sum of all the Ca
# function(v) CaTminf(v) + CaSminf(v);
CaI - gCaT*CaTm^3*CaTh*(v-ECa(CaI)) + gCaS*CaSm^3*(v-ECa(CaI))

#AB
#dCa2 - (((-F*Caminf(v))-Caminf(v) + C0)/TauCa)
dCa2 - (((-F*CaI) - Ca2 + C0)/TauCa)

# mk20
KCaminf - function(v, Ca2) (Ca2/(Ca2+30))*(1/(1+exp(-(v+51)/8)));
# taumk
KCataum - function(v) 90.3 - ((75.09/(1+exp(-(v+46)/22.7;

#AB
Aminf - function(v) 1/(1+exp(-(v+27)/8.7));
Ahinf - function(v) 1/(1+exp((v+56.9)/4.9));
Ataum - function(v) 11.6-(10.4/(1+exp(-(v+32.9)/15.2)));
Atauh - function(v) 38.6-(29.2*(1+exp(-(v+38.9)/26.5)));

#proc
#mp0
procminf - function(v) 1/(1+exp((v+56.9)/4));
#taump
proctaum - function(v) 0.5;


 dv - (-1*(I
   + CaI
   + gNap*Napm^3*Naph*(v-ENap)
   + gh*hm*(v-Eh)
   + gK*Km^4*(v-EK)
   + gKCa * KCam^4*(v-EKCa)
   + gA*Am^4*Ah*(v-EA)
   + gL*(v-EL))
   / C);


dCaTm - (CaTminf(v) - CaTm)/CaTtaum(v);
dCaTh - (CaThinf(v) - CaTh)/CaTtauh(v);

dCaSm - (CaSminf(v) - CaSm)/CaStaum(v);

dNapm - (Napminf(v) - Napm)/taumna(v);
dNaph - (Napminf(v) - Naph)/tauhna(v);

dhm - (hminf(v) - hm)/htaum(v);

dKm - (Kminf(v) - Km)/Ktaum(v);

dKCam - (KCaminf(v, Ca2) - KCam)/KCataum(v);

dAm - (Aminf(v) - Am)/Ataum(v);
dAh - (Ahinf(v) - Ah)/Atauh(v);


list(c(dv,
   dCaTm, dCaTh,
   dCaSm,
   dNapm, dNaph,
   dhm,
   dKm,
   dKCam,
   dAm, dAh))
  })
}
parms = c(gCaT=22.5, gCaS=60, gNap=4.38, gh=0.219,
  gK=1576.8, gKCa=251.85, gA=39.42, gL=0.105,
  ENap=50, Ca2=0.52,
  Eh=-20, EK=-80, EL=-55, EKCa=-80, EA=-80,
  C=1/12, I=6.5, F=0.418, TauCa=303, C0=0.5, CaI=1);
time = seq(from=0, to=1000, by=0.1);
init = c(v=-65, CaTm=0.52 , CaTh=0.52, CaSm=0.52,
 Napm=0.52, Naph=0.52, hm=0.52, Km=0.52,
 KCam=0.52, Am=0.52, Ah=0.52);


out-ode(y=init, times=time, func=HH_soma, parms=parms);
plot(ode)
o-data.frame(out);
plot(o$time, o$v, type='l');


Regards
Jannetta


On 15 February 2013 09:46, Berend Hasselman b...@xs4all.nl wrote:


 On 15-02-2013, at 10:28, Jannetta Steyn janne...@henning.org wrote:

  Hi Ben
 
  Thank you so much for spending the time to look at my code. I really
  appreciate it.
 
  The unnecessary parameters were artefacts from me working on the model
 and
  putting things in and taking things out.
 
  The model isn't quite producing what I expect yet, but it is definitely
  starting to look more interesting :-)
 
  I am now getting this error:
 
  plot(ode)
  Error in is.null(times) : 'times' is missing
 
  What on earth does it mean because I cant see any 'times' being missing.
 To
  me it seems that all is there?!?
 

 Difficult to tell as you haven't shown how ode was defined.

 I just did this:

 out-ode(y=init, times=times, func=HH_soma, parms=parms)
 plot(out)

 and got the plots of all state variables.

 Berend





-- 

===
Web site: http://www.jannetta.com
Email: janne...@henning.org
===

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list

Re: [R] Troubleshooting underidentification issues in structural equation modelling (SEM)

2013-02-15 Thread Bert Gunter
These are statistical, not R issues, so please do not post further
here. You are clearly out of your depth statistically. You need to get
local statistical help, or you can try posting on a statistical list
like stats.stackexchange.com if you care to take advice from unknown
sources who don't understand the details of your situation.

-- Bert

On Fri, Feb 15, 2013 at 1:11 AM, Ruijie breakaw...@gmail.com wrote:
 Thanks Prof Fox for your guidance. My purpose in fitting this model is to
 contrast it with another model that I am proposing which I believe will be
 a better fit.

 On the point of some of the items being close to invariant, I had a close
 look at my data and indeed that is the case I am aware of it. However, I am
 not sure what to do with these items. Do I remove them? If I do, what
 threshold of variance do I set for removal? How do I decide on that
 threshold?

 I've combed a number of textbooks for answers but sadly have not found
 much. Hope you could offer some advice, thanks!

 Regards,
 Ruijie (RJ)

 
 He who has a why can endure any how.

 ~ Friedrich Nietzsche


 On 10 February 2013 00:38, John Fox j...@mcmaster.ca wrote:

 Dear Ruijie,

 Your model is underidentified by virtue of two of the factors having only
 one observed indicator each. No SEM software can magically estimate this
 model as it stands. Beyond that, I won't comment on the wisdom of what
 you're doing, such as computing covariances between ordinal variables --
 but
 see what I discovered below.

 Removing these two variables and the associated factors produces the
 following model:

 - snip 

  model - cfa(reference.indicators=FALSE)
 1: F01: I01, I02, I03
 2: F02: I04, I05, I06, I07, I08, I09, I10, I11, I12, I13
 3: F03: I14, I15, I16, I17, I18, I19, I20, I21, I22, I23, I24, I25, I26
 4: F04: I27, I28, I29, I30, I31, I32, I33, I34
 5: F05: I35, I36, I37, I38, I39, I40, I41, I42, I43
 6: F07: I46, I47, I48, I49, I50, I51
 7: F08: I54, I55, I56, I57, I58, I59, I60, I61, I62, I63, I64
 8: F09: I65, I66, I67
 9: F11: I69, I70, I71
 10:
 Read 9 items
 NOTE: adding 66 variances to the model
 
  cfa.output - sem(model, cov.mat, N = 900)

 - snip 

 sem() ran out of iterations, but the summary output is revealing:

 - snip 

  summary(cfa.output)

  Model Chisquare =  5677.1   Df =  2043 Pr(Chisq) = 0
  AIC =  6013.1
  BIC =  -8220.193

  Normalized Residuals
Min. 1st Qu.  MedianMean 3rd Qu.Max.
 -3.9910 -0.5887 -0.1486  0.2588  0.8092 17.2900

  R-square for Endogenous Variables
 I01 I02 I03 I04 I05 I06 I07 I08 I09
 I10
  0.0953  0.1263  0.  0.1131  0.4039  0.2519  0.1168  0.0468  0.0005
 0.0059
 I11 I12 I13 I14 I15 I16 I17 I18 I19
 I20
  0.0479  0.0228  0.1150  0.2813  0.0001  0.0388  0.2106  0.0001  0.0913
 0.0063
 I21 I22 I23 I24 I25 I26 I27 I28 I29
 I30
  0.0041  0.0077  0.0022  0.  0.0299  0.0067  0.0019  0.0011  0.0010
 0.
 I31 I32 I33 I34 I35 I36 I37 I38 I39
 I40
  0.0005  0.0117  0.0270  0.0001  0.0084  0.0001  0.0256  0.4969  0.0613
 0.0515
 I41 I42 I43 I46 I47 I48 I49 I50 I51
 I54
  0.0005  0.0052  0.0307  0.0003  0.1131  0.0014  0.  0.1276  0.9728
 0.0520
 I55 I56 I57 I58 I59 I60 I61 I62 I63
 I64
  0.2930  0.0127  0.0543  0.0500  0.0378  0.0001  0.3048  0.0002  0.0304
 0.0001
 I65 I66 I67 I69 I70 I71
 56.7264  0.  0.0002  0.2220  0.2342  0.2240

  Parameter Estimates
  Estimate  Std Errorz value  Pr(|z|)

 lam[I01:F01]  3.023074e-02 5.133785e-03  5.888586224  3.895133e-09 I01 ---
 F01
 lam[I02:F01]  3.283192e-02 5.291069e-03  6.205157975  5.464199e-10 I02 ---
 F01
 lam[I03:F01]  1.123398e-04 2.695713e-03  0.041673509  9.667590e-01 I03 ---
 F01
 lam[I04:F02]  1.365329e-01 1.555023e-02  8.780124358  1.632940e-18 I04 ---
 F02
 lam[I05:F02]  9.525580e-02 5.517838e-03 17.263245517  8.896692e-67 I05 ---
 F02
 lam[I06:F02]  1.720147e-01 1.277593e-02 13.463962882  2.548717e-41 I06 ---
 F02
 lam[I07:F02]  3.164280e-02 3.543421e-03  8.930015663  4.259485e-19 I07 ---
 F02
 lam[I08:F02]  5.685988e-02 1.021854e-02  5.564386503  2.630763e-08 I08 ---
 F02
 lam[I09:F02]  1.234516e-03 2.228298e-03  0.554017268  5.795670e-01 I09 ---
 F02
 lam[I10:F02]  1.656005e-02 8.458411e-03  1.957820181  5.025112e-02 I10 ---
 F02
 lam[I11:F02]  8.785114e-02 1.560646e-02  5.629151062  1.810987e-08 I11 ---
 F02
 lam[I12:F02]  3.022114e-02 7.815459e-03  3.866842129  1.102537e-04 I12 ---
 F02
 lam[I13:F02]  5.075487e-02 5.732307e-03  8.854177302  8.430329e-19 I13 ---
 F02
 lam[I14:F03]  2.587670e-01 2.308125e-02 11.211137448  3.595430e-29 I14 ---
 F03
 lam[I15:F03] -2.999816e-04 1.469667e-03 -0.204115351  8.382634e-01 I15 ---
 F03
 lam[I16:F03]  2.314973e-02 5.256310e-03  4.404179628  

Re: [R] An extended Hodgkin-Huxley model that doesn't want to work.

2013-02-15 Thread Jannetta Steyn
Hi Ben

Thank you so much for spending the time to look at my code. I really
appreciate it.

The unnecessary parameters were artefacts from me working on the model and
putting things in and taking things out.

The model isn't quite producing what I expect yet, but it is definitely
starting to look more interesting :-)

I am now getting this error:

 plot(ode)
Error in is.null(times) : 'times' is missing

What on earth does it mean because I cant see any 'times' being missing. To
me it seems that all is there?!?

Kind Regards
Jannetta


On 14 February 2013 17:41, Berend Hasselman b...@xs4all.nl wrote:


 Forgot Reply to All.

 On 13-02-2013, at 23:30, Jannetta Steyn janne...@henning.org wrote:

  Hi All
 
  I have been struggling with this model for some time now and I just can't
  get it to work correctly. The messages I get when running the code is:
 
  DLSODA-  Warning..Internal T (=R1) and H (=R2) are
  such that in the machine, T + H = T on the next step
 (H = step size). Solver will continue anyway.
  In above message, R =
  [1] 0 0
  DINTDY-  T (=R1) illegal
  In above message, R =
  [1] 0.1
  T not in interval TCUR - HU (= R1) to TCUR (=R2)
  In above message, R =
  [1] 0 0
  DINTDY-  T (=R1) illegal
  In above message, R =
  [1] 0.2
  T not in interval TCUR - HU (= R1) to TCUR (=R2)
  In above message, R =
  [1] 0 0
  DLSODA-  Trouble in DINTDY.  ITASK = I1, TOUT = R1
  In above message, I =
  [1] 1
  In above message, R =
  [1] 0.2
  Error in lsoda(y, times, func, parms, ...) :
  illegal input detected before taking any integration steps - see written
  message
 
 
 
  I'll first paste the formulae and then I'll paste my code. If anyone can
  spot something wrong with my implementation it would really make my day.
 
  (1)
  dV/dt = (I_ext - I_int-I_coup)/C
  I_ext = injected current
  I_int = Sum of all ion currents
  I_coup = coupling current (but we're not using it here )
 
  (2)
  I_i = g_i * m_i^pi * h_i^pi(V-E)
  i identifies the ion, thus I_K would be Potassium current.
 
  (3)
  dm/dt = (m_inf*V - m)/tau_m
 
  (4)
  dh/dt = (h_inf*V-h)/tau_h
 
  (5)
  The Nernst equation is used to calculate reversal potential for Ca:
  Eca = 12.2396 * log(13000/Ca2+)
 
  (6)
  d[Ca_2+]/dt = (F*I_Ca - [Ca2+] + C0)/Tau_Ca
 
 
  tau_m, tau_h, m_inf and h_inf are all calculated according to formulae
  provided in a paper. In my code these are calculated for the different
  channels into the following variables:
 
  CaTminf, CaThinf, CaTtaum, CaTtauh, CaSminf, CaStaum, Napminf, Naphinf,
  taumna, tauhna, hminf, htaum, Kminf and Ktaum
 
  The E (reversal potential) values for all the channels are given, except
  for CaT and CaS which uses Eca as calculated in (5).
 
  Current for Ca is calculated by summing the CaT and CaS currents, hence
  CaI = gCaT*CaTm^3*CaTh*(v-Eca(v)) + gCaS*CaSm^3(v-ECa(v)
 
 
  Here is the code:
 
  library(simecol)
  ## Hodkin-Huxley model
  HH_soma -  function(time, init, parms) {
  with(as.list(c(init, parms)),{
 
# Na only used in Axon
#Naminf -1/(1+exp(-(v+24.7)/5.29));
#Nataum - function(v) 1.32 - (1.26/(1+exp(-(v+120)/25)));
#Nahinf -1/(1+exp((v+489)/5.18));
#Natauh -(0.67/(1+exp(-(v+62.9)/10))) *
 (1.5+(1/(1+exp((v+34.9)/36;
 
#PD
# mca10
CaTminf - function(v) 1/(1+exp(-(v+25)/7.2));
# hca10
CaThinf - function(v) 1/(1+exp(v+36)/7);
# taumca1
CaTtaum - function(v) 55- (49.5/(1+exp(-v+58)/17));
# tauhca1
CaTtauh - function(v) 350 - (300/(1+exp(-v+50)/16.9));
 
#mca20
CaSminf - function(v) 1/(1+exp(-(v+22)/8.5));
#taumca2
CaStaum - function(v) 16-(13.1/(1+exp(-(v+25.1)/26.4)));
 
 
# mna0
Napminf - function(v) 1/(1+exp(-(v+26.8)/8.2));
# hna0
Naphinf - function(v) 1/(1+exp(-(v+48.5)/5.18));
 
taumna - function(v) 19.8-(10.7/(1+exp(-(v+26.5)/8.6)));
tauhna - function(v) 666-(379/(1+exp(-(v+33.6)/11.7)));
 
# mh0
hminf - function(v)  1/(1+exp(v+70)/6);
# taumh
htaum - function(v)  272+(1499/(1+exp(-(v+42.2)/8.73)));
 
Kminf - function(v)  1/(1+exp(-(v+14.2)/11.8));
Ktaum - function(v)  7.2-(6.4/(1+exp(-(v+28.3)/19.2)));
 
# Reversal potential of intracellular calcium concentration
# Nernst Equation using  extracellular concentration of Ca = 13mM
# eca
ECa - function(Ca2) 12.2396*log(13000/(Ca2));
#ECa - function(CaI) 12.2396*log(13000/(CaI));
 
 
#Sum of all the Ca
# function(v) CaTminf(v) + CaSminf(v);
CaI - gCaT*CaTm^3*CaTh*(v-ECa(CaI)) + gCaS*CaSm^3*(v-ECa(CaI))
 
#AB
#dCa2 - (((-F*Caminf(v))-Caminf(v) + C0)/TauCa)
dCa2 - (((-F*CaI) - Ca2 + C0)/TauCa)
 
# mk20
KCaminf - function(v, Ca2) (Ca2/(Ca2+30))*(1/(1+exp(-(v+51)/8)));
# taumk
KCataum - function(v) 90.3 - ((75.09/(1+exp(-(v+46)/22.7;
 
#AB
Aminf - function(v) 1/(1+exp(-(v+27)/8.7));
Ahinf - function(v) 1/(1+exp((v+56.9)/4.9));
Ataum - function(v) 11.6-(10.4/(1+exp(-(v+32.9)/15.2)));
Atauh - function(v) 

Re: [R] plot custom x axis ticks values

2013-02-15 Thread Jim Lemon

On 02/15/2013 07:28 PM, e-letter wrote:

...
plot(testvalues,ann=FALSE,type='l',yaxt='n',xaxt='n',xlim=c(0,300))
axis(1,at=c(0,100,200,300))

The use of 'xlim' makes the graph unacceptable. Instead, it would be
better to plot:

plot(testvalues,ann=FALSE,type='l',yaxt='n')

The x-axis shows tick marks '1.0', '1.5', ...

It is required to replace these values with custom values (e.g. 10, 20... etc.)


Okay, I think I see what you want:

plot(testvalues,ann=FALSE,type='l',yaxt='n')
axis(2,at=c(1,1.5,2,2.5,3),labels=c(10,20,30,40,50))

Jim

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


Re: [R] An extended Hodgkin-Huxley model that doesn't want to work.

2013-02-15 Thread Berend Hasselman

On 15-02-2013, at 11:22, Jannetta Steyn janne...@henning.org wrote:

 Hi Berend
 
 This is the code. Pretty much just changed to what you suggested which is
 CaI=1, removing the unnecessary variables and using deSolve:
 
 rm(list = ls())
 library(deSolve)
 ## Hodkin-Huxley model
 HH_soma -  function(times, init, parms) {
  with(as.list(c(init, parms)),{
 
# Na only used in Axon
#Naminf -1/(1+exp(-(v+24.7)/5.29));
#Nataum - function(v) 1.32 - (1.26/(1+exp(-(v+120)/25)));
#Nahinf -1/(1+exp((v+489)/5.18));
#Natauh -(0.67/(1+exp(-(v+62.9)/10))) * (1.5+(1/(1+exp((v+34.9)/36;
 
#PD
# mca10
CaTminf - function(v) 1/(1+exp(-(v+25)/7.2));
# hca10
CaThinf - function(v) 1/(1+exp(v+36)/7);
# taumca1
CaTtaum - function(v) 55- (49.5/(1+exp(-v+58)/17));
# tauhca1
CaTtauh - function(v) 350 - (300/(1+exp(-v+50)/16.9));
 
#mca20
CaSminf - function(v) 1/(1+exp(-(v+22)/8.5));
#taumca2
CaStaum - function(v) 16-(13.1/(1+exp(-(v+25.1)/26.4)));
 
 
# mna0
Napminf - function(v) 1/(1+exp(-(v+26.8)/8.2));
# hna0
Naphinf - function(v) 1/(1+exp(-(v+48.5)/5.18));
 
taumna - function(v) 19.8-(10.7/(1+exp(-(v+26.5)/8.6)));
tauhna - function(v) 666-(379/(1+exp(-(v+33.6)/11.7)));
 
# mh0
hminf - function(v)  1/(1+exp(v+70)/6);
# taumh
htaum - function(v)  272+(1499/(1+exp(-(v+42.2)/8.73)));
 
Kminf - function(v)  1/(1+exp(-(v+14.2)/11.8));
Ktaum - function(v)  7.2-(6.4/(1+exp(-(v+28.3)/19.2)));
 
# Reversal potential of intracellular calcium concentration
# Nernst Equation using  extracellular concentration of Ca = 13mM
# eca
ECa - function(Ca2) 12.2396*log(13000/(Ca2));
#ECa - function(CaI) 12.2396*log(13000/(CaI));
 
 
#Sum of all the Ca
# function(v) CaTminf(v) + CaSminf(v);
CaI - gCaT*CaTm^3*CaTh*(v-ECa(CaI)) + gCaS*CaSm^3*(v-ECa(CaI))
 
#AB
#dCa2 - (((-F*Caminf(v))-Caminf(v) + C0)/TauCa)
dCa2 - (((-F*CaI) - Ca2 + C0)/TauCa)
 
# mk20
KCaminf - function(v, Ca2) (Ca2/(Ca2+30))*(1/(1+exp(-(v+51)/8)));
# taumk
KCataum - function(v) 90.3 - ((75.09/(1+exp(-(v+46)/22.7;
 
#AB
Aminf - function(v) 1/(1+exp(-(v+27)/8.7));
Ahinf - function(v) 1/(1+exp((v+56.9)/4.9));
Ataum - function(v) 11.6-(10.4/(1+exp(-(v+32.9)/15.2)));
Atauh - function(v) 38.6-(29.2*(1+exp(-(v+38.9)/26.5)));
 
#proc
#mp0
procminf - function(v) 1/(1+exp((v+56.9)/4));
#taump
proctaum - function(v) 0.5;
 
 
 dv - (-1*(I
   + CaI
   + gNap*Napm^3*Naph*(v-ENap)
   + gh*hm*(v-Eh)
   + gK*Km^4*(v-EK)
   + gKCa * KCam^4*(v-EKCa)
   + gA*Am^4*Ah*(v-EA)
   + gL*(v-EL))
   / C);
 
 
dCaTm - (CaTminf(v) - CaTm)/CaTtaum(v);
dCaTh - (CaThinf(v) - CaTh)/CaTtauh(v);
 
dCaSm - (CaSminf(v) - CaSm)/CaStaum(v);
 
dNapm - (Napminf(v) - Napm)/taumna(v);
dNaph - (Napminf(v) - Naph)/tauhna(v);
 
dhm - (hminf(v) - hm)/htaum(v);
 
dKm - (Kminf(v) - Km)/Ktaum(v);
 
dKCam - (KCaminf(v, Ca2) - KCam)/KCataum(v);
 
dAm - (Aminf(v) - Am)/Ataum(v);
dAh - (Ahinf(v) - Ah)/Atauh(v);
 
 
list(c(dv,
   dCaTm, dCaTh,
   dCaSm,
   dNapm, dNaph,
   dhm,
   dKm,
   dKCam,
   dAm, dAh))
  })
 }
 parms = c(gCaT=22.5, gCaS=60, gNap=4.38, gh=0.219,
  gK=1576.8, gKCa=251.85, gA=39.42, gL=0.105,
  ENap=50, Ca2=0.52,
  Eh=-20, EK=-80, EL=-55, EKCa=-80, EA=-80,
  C=1/12, I=6.5, F=0.418, TauCa=303, C0=0.5, CaI=1);
 time = seq(from=0, to=1000, by=0.1);
 init = c(v=-65, CaTm=0.52 , CaTh=0.52, CaSm=0.52,
 Napm=0.52, Naph=0.52, hm=0.52, Km=0.52,
 KCam=0.52, Am=0.52, Ah=0.52);
 
 
 out-ode(y=init, times=time, func=HH_soma, parms=parms);
 plot(ode)
 o-data.frame(out);
 plot(o$time, o$v, type='l');


You are not doing what I told you to do.

You should plot the result of ode not ode itself (that's the function).
So do

plot(out)

and you will get the plots.

Berend

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


Re: [R] plot custom x axis ticks values

2013-02-15 Thread e-letter
On 15/02/2013, Duncan Murdoch murdoch.dun...@gmail.com wrote:
 On 13-02-15 3:28 AM, e-letter wrote:
 On 15/02/2013, Jim Lemon j...@bitwrit.com.au wrote:
 On 02/14/2013 09:41 PM, e-letter wrote:
 Readers,

 For this data set:

 testvalues-c(10,20,30,40)

 How to amend the plot instruction:

 plot(testvalues,ann=FALSE,type='l',yaxt='n',xaxt='n')

 so that x axis ticks labels can be added to existing graph with
 arbitrary value such as 0,100,200,300)?

 Hi r2151,
 If you want the labels to fit on the axis you will have to include this
 information in the call to plot:

 plot(testvalues,ann=FALSE,type='l',yaxt='n',xaxt='n',ylim=c(0,300))
 axis(2,at=c(0,100,200,300))


 plot(testvalues,ann=FALSE,type='l',yaxt='n',xaxt='n',xlim=c(0,300))
 axis(1,at=c(0,100,200,300))

 The use of 'xlim' makes the graph unacceptable. Instead, it would be
 better to plot:

 plot(testvalues,ann=FALSE,type='l',yaxt='n')

 The x-axis shows tick marks '1.0', '1.5', ...

 It is required to replace these values with custom values (e.g. 10, 20...
 etc.)

 Why not rescale the values before plotting, and use the automatic ticks?
   You can lie about the user scale, but it doesn't always give a helpful
 plot, e.g.

 axis(1, at=1:4, labels=c(100, 300, 400, 200))


Tried, but realised that this works when the vectors for 'at' and
'labels' are equal size. For a data set of e.g. 10 values, syntax
error returns.

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


Re: [R] lm regression query

2013-02-15 Thread Bert Gunter
Smells like homework to me. If so, we don't do homework on this list.

-- Bert

On Thu, Feb 14, 2013 at 3:55 PM, email email8...@gmail.com wrote:
 Hello:

 I have a 4-column dataset: Crime, Education, Urbanization, Age. I want to
 construct a multiple linear regression to find the effect of Education,
 Urbanization, and Age on Crime

 lm(Crime ~ Education + Urbanization + Age)

 If I use + in above statement, does it mean it will build a model to find
 the relationship between Crime and Education when Urbanization and Age are
 held constant?

 What would be the difference if I drop the term Urbanization + Age ?

 lm(Crime ~ Education)

 Regards:
 John

 [[alternative HTML version deleted]]

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



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

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

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

2013-02-15 Thread Giovanna Ottaviani
Hello everybody,
Anyone familiar with the package frontier? I have some general questions on how 
to approach the model design.
Thanks in advance

Giovanna
Giovanna Ottaviani Aalmo
Stipendiat/Ph.D. Student
---
Norsk institutt for skog og landskap
Pb 115, NO-1431 Ås
T (+47) 64 94 9094
M(+47) 980 30 422
F(+47) 64 94  90 80
---
www.skogoglandskap.nohttp://www.skogoglandskap.no/
---


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

2013-02-15 Thread Arne Henningsen
Dear Giovanna

On 15 February 2013 11:50, Giovanna Ottaviani g...@skogoglandskap.no wrote:
 Anyone familiar with the package frontier? I have some general
 questions on how to approach the model design.

The R package frontier is based on FRONTIER 4.1 [1] but it has
some improvements [2]. The models that can be estimated by the R
package frontier are the same as the models that can be estimated by
FRONTIER 4.1. These models are explained in the file Front41.pdf
that is included in the archive FRONT41-xp1.zip, which is available
at [1]. If you have any further questions regarding the frontier
package, please use a forum at the package's R-Forge site [3].

[1] http://www.uq.edu.au/economics/cepa/frontier.php

[2] http://frontier.r-forge.r-project.org/

[3] http://r-forge.r-project.org/projects/frontier/

Best wishes from Copenhagen,
Arne

--
Arne Henningsen
http://www.arne-henningsen.name

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


Re: [R] 2 setGeneric's, same name, different method signatures

2013-02-15 Thread Greg Minshall
Martin,

fantastic.  thank you *very* much!  that clears lots of things up for
me.

(for the record: i think that setGeneric overwriting a previous is more
surprising -- thus violating the principle of least surprise -- than one
function overwriting a previous, in that we think of (or, one way to
think of) OOP is that, after finagling to have no clashes on *class*
names, we should not worry about other than intra-class name clashes.
as i said, for the record.)

thanks again.

cheers, Greg

Martin Morgan mtmor...@fhcrc.org wrote:
 
 Hi Greg --
 
 this setGeneric is over-writing the first, as would
 
 f = function() first
 f = function() second
 f() # second
 
 If you'd like to dispatch on a single argument, then
 
 setGeneric(one, function(x, ...) standardGeneric(one))
 setMethod(one, A, function(x, ...) A-method)
 setMetohd(one, B, function(x, y, ...) B-method)
 
 The '...' in the generic allow you to add arguments that are 'picked
 off' by methods. The user could provide any value for y, not only an
 object of class B.
 
 If you'd like to dispatch sometimes on two arguments then
 
 setGeneric(two, function(x, y, ...) standardGeneric(two))
 setMethod(two, c(A, ANY), function(x, y, ...) A,ANY-method)
 setMethod(two, c(B, B), function(x, y, ...) B,B-method)
 
 then two(new(A)), two(new(A), new(A)) and two(new(A),
 new(B)) end up in A,ANY,two-method while two(new(B), new(B))
 ends up in B,B,two-method. Other combinations are errors. One might
 instead not define A,ANY but instead
 
 setMethod(two, c(A, missing), function(x, y, ...) A,missing-method)
 
 and then two(new(A), new(A)) would be an error.
 
 Multiple dispatch is complicated, and perhaps best to avoid if possible.
 
 It's possible to write a generic and methods that dispatch on '...',
 with the requirement that all classes are the same; this in the spirit
 of comparing two B's and returning the smaller; see ?dotsMethods
 though again this is not a trivial use case.
 
 Hope that helps enough.
 
 Martin

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

2013-02-15 Thread Amir Kasaeian
Dear all,
Good day!

I have a question in my codes, would you please help me how to rectify it?
these are my coded but at the last line I received the error!

mac_30 = read.dta(MAC results4.dta)
mac_30
map_30 = read.dta(MAP results4.dta)
map_30
mac_30$weight = 1/nrow(mac_30)
mac_30
map_30$weight = 1/nrow(map_30)
map_30
input.data = merge(mac_30, map_30, all = TRUE)
input.data.no15to19 = subset(input.data, agegroup != 15-19 | is.na(agegroup))
input.data.no15to19
fit = loess(v5q0 ~ year, input.data.no15to19, weights = weight, span = 0.5)

years.predict = seq(min(input.data.no15to19$year), 
max(input.data.no15to19$year))

combined.method = data.frame(iso3 = unique(input.data$iso3), svdate = 
unique(input.data$svdate), method = Combined method, year = years.predict)
the error is:

Error in data.frame(iso3 = unique(input.data3$iso3), svdate 
=unique(input.data3$svdate),  :
  arguments imply differing number of rows: 1, 4, 45


I know there is only one iso3 and 4 years but what I need is forty five 1 and 
forty five rows of different 4 svdate exactly equivalent to forty five rows of 
different years.
 Now, how can I have the same number of rows for the different variable to 
continue my code for rest of analysis?

Thank you so much.

Kind regards,
Amir


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


Re: [R] 3D-plots of 2D-grids

2013-02-15 Thread jas
Hello David,
thanks again for your reply.

Two things remain unclear. That the data is disjointed is ok, as there are
only values in hectares, where there are actually buildings and stuff,
forest/nature is NA. 

The presp-scan that I get from persp(grd) has no similarity to the image in
2d. My guess is that the order and regularity in the data somehow gets lost
in the process of making the matrix?! 

The grid-extraction that I uploaded contains 40x20 cells. if its a bigger
grid with 7million cells, does it still work the same way or can the vectors
only be a maximum of 100 cells? 
 
Thank you for your help
jas



--
View this message in context: 
http://r.789695.n4.nabble.com/3D-plots-of-2D-grids-tp4658517p4658639.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] post hoc comparisons on random effects

2013-02-15 Thread Vasillis Papathanasiou

Dear all,
is it possible to run a post hoc comparison using the glht function on random 
factors from a mixed model?
For example, in the model lie(Yield~1+Meadow, random=1|Site/Area/Quadrate).
thank you,

Vasillis Papathanasiou,
Marine SciencesPalaion Fokon 35-37Thessaloniki, 54454Greecetel:+306945874855
  
[[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] Using paste on results from tapply?

2013-02-15 Thread arun
HI,

You could also use these:
sapply(groups,paste,collapse= )
#  1   3   6 
 #   6   4 4 5 3 2 
#or
gsub(^ | $,,gsub(\\D+, ,paste(groups)))
#[1] 6 4 4   5 3 2
A.K.





- Original Message -
From: Nikola Janevski njane...@gmail.com
To: r-help@r-project.org
Cc: 
Sent: Thursday, February 14, 2013 8:28 PM
Subject: [R] Using paste on results from tapply?

Hi,

This is what I'm trying to do:
1. I have a vector split.indexes - c(1, 3, 3, 6, 6, 6)
2. I have another vector with something like data - c(6, 4, 4, 5, 3, 2)
3. I use groups - tapply(data, split.indexes, c) to create a list of
vectors based on the levels in split.indexes. As a result groups has the
following values:
$`1`
[1] 6

$`3`
[1] 4 4

$`6`
[1] 5 3 2

This is the desired result.

4. The next thing I want to do is call paste on groups so I can get a
string representation that I can write to a file. However, when I call
paste(groups) it returns the following:
[1] 6          c(4, 4)    c(5, 3, 2)

I do not understand why the paste function is adding the c() part. I
would like the result to be:
[1] 6          4 4    5 3 2

Do you know any efficient way to get rid of the c()?

Note: I'm working with a very large data set 1-10 millions of points and I
would rather not use for loops since the performance are really bad.

Sincerely,
Nikola

    [[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] match a task No. to a few person-IDs

2013-02-15 Thread Mat
hello together,

i have a task No. in a data.frame
   Task_No   Team 
A49397   1 
B49396   1  

and now i want to match my Person-IDs to each Task_No. My Person Ids look
like this one:
IDTeam
A 5019  1
B 472   1
C 140   1
D 5013  1

The solution should be this:

   Task_No   Team ID
A49397   1 5019
A49397   1 472
A49397   1 140
A49397   1 5013
B49396   1 5019
B49396   1472
B49396   1140
B49396   15013

perhabs anyone can help me.

Thank you

Mat



--
View this message in context: 
http://r.789695.n4.nabble.com/match-a-task-No-to-a-few-person-IDs-tp4658646.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] Fitting pareto distribution / plotting observed fitted dists

2013-02-15 Thread Stefan Sobernig
Some background: I have some data on structural dependencies in a base 
of code artifacts. The dependency structure is reflected in terms of 
relative node degrees, with each node representing some code unit (just 
as an example).


This gives me real data of the following form (sorry for the longish 
posting):


dat1 -
c(0.00245098039215686, 0, 0, 0, 0, 0, 0, 0, 0.0563725490196078,
0, 0, 0, 0.071078431372549, 0, 0, 0, 0, 0, 0, 0, 0.00490196078431373,
0.00245098039215686, 0, 0, 0, 0.00245098039215686, 0.0196078431372549,
0, 0.0588235294117647, 0.00490196078431373, 0.00980392156862745,
0, 0, 0.00245098039215686, 0, 0, 0, 0.0220588235294118, 0, 0,
0.00245098039215686, 0, 0, 0.00245098039215686, 0, 0.00245098039215686,
0.00980392156862745, 0.0294117647058824, 0.0637254901960784,
0, 0, 0, 0.00245098039215686, 0.0392156862745098, 0, 0.0147058823529412,
0, 0.017156862745098, 0.00245098039215686, 0, 0.00980392156862745,
0, 0.00735294117647059, 0.00490196078431373, 0.0514705882352941,
0.00245098039215686, 0, 0, 0, 0, 0.00245098039215686, 0, 
0.00245098039215686,

0, 0.00245098039215686, 0, 0, 0.0147058823529412, 0.0367647058823529,
0.0269607843137255, 0.0269607843137255, 0.00735294117647059,
0.0441176470588235, 0, 0, 0.0196078431372549, 0, 0.00490196078431373,
0.0245098039215686, 0.00490196078431373, 0.00490196078431373,
0.0196078431372549, 0, 0.0318627450980392, 0.0245098039215686,
0, 0.00245098039215686, 0, 0.0417, 0, 0, 0.00490196078431373,
0.00490196078431373, 0.00245098039215686, 0.0122549019607843,
0.00490196078431373, 0.00490196078431373, 0.071078431372549,
0, 0, 0, 0, 0.00245098039215686, 0.00245098039215686, 0.00490196078431373,
0.0343137254901961, 0.00980392156862745, 0.00245098039215686,
0.053921568627451, 0, 0.0245098039215686, 0.00245098039215686,
0.0245098039215686, 0, 0.0294117647058824, 0.00490196078431373,
0.00980392156862745, 0.0367647058823529, 0, 0, 0.017156862745098,
0, 0.0245098039215686, 0, 0.071078431372549, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00980392156862745, 0.00980392156862745,
0.00490196078431373, 0.0343137254901961, 0.0147058823529412,
0.0122549019607843, 0.00735294117647059, 0, 0.00245098039215686,
0, 0.00245098039215686, 0.110294117647059, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0.00490196078431373, 0.00735294117647059, 0, 0.00245098039215686,
0, 0, 0.00735294117647059, 0.0735294117647059, 0, 0, 0, 0, 0,
0.00490196078431373, 0, 0.00245098039215686, 0.105392156862745,
0, 0, 0, 0, 0, 0.00735294117647059, 0.00980392156862745, 
0.00245098039215686,

0.0147058823529412, 0, 0, 0.00245098039215686, 0.00490196078431373,
0.00245098039215686, 0.00735294117647059, 0.0563725490196078,
0, 0, 0, 0, 0.0318627450980392, 0, 0, 0.00490196078431373, 
0.0465686274509804,

0, 0.0147058823529412, 0.017156862745098, 0.00735294117647059,
0.0245098039215686, 0.017156862745098, 0, 0, 0, 0.071078431372549,
0, 0, 0, 0, 0, 0, 0, 0.00980392156862745, 0.0122549019607843,
0.00980392156862745, 0.0196078431372549, 0, 0, 0, 0.00980392156862745,
0.00490196078431373, 0.00735294117647059, 0, 0.0196078431372549,
0.0220588235294118, 0, 0, 0.00490196078431373, 0.0661764705882353,
0, 0, 0, 0, 0, 0, 0, 0, 0.00490196078431373, 0.0955882352941176,
0, 0, 0, 0, 0, 0.00490196078431373, 0.00490196078431373, 0, 
0.00245098039215686,

0.0588235294117647, 0, 0, 0, 0, 0.0857843137254902, 0, 0, 0,
0, 0.017156862745098, 0, 0, 0.0294117647058824, 0, 0, 0.0122549019607843,
0, 0, 0.0147058823529412, 0, 0, 0.0318627450980392, 0, 0, 
0.00245098039215686,

0.00245098039215686, 0.00980392156862745, 0.00245098039215686,
0.00245098039215686, 0.00245098039215686, 0.0784313725490196,
0, 0, 0.0784313725490196, 0, 0, 0, 0.00735294117647059, 
0.00245098039215686,

0.00490196078431373, 0.00245098039215686, 0, 0.00245098039215686,
0, 0.0122549019607843, 0, 0.0857843137254902, 0, 0, 0, 0, 
0.0196078431372549,

0, 0.00980392156862745, 0.00245098039215686, 0, 0, 0.00490196078431373,
0.00735294117647059, 0.00735294117647059, 0.00735294117647059,
0.00735294117647059, 0.00735294117647059, 0.00735294117647059,
0.00735294117647059, 0, 0, 0.00735294117647059, 0.00980392156862745,
0.0122549019607843, 0.0245098039215686, 0.00980392156862745,
0.00245098039215686, 0.00490196078431373, 0.00490196078431373,
0.00245098039215686, 0.00245098039215686, 0.00735294117647059,
0.00490196078431373, 0, 0.00245098039215686, 0.017156862745098,
0.00490196078431373, 0.00490196078431373, 0.00245098039215686,
0.00245098039215686, 0.0269607843137255, 0, 0.00490196078431373,
0.00490196078431373, 0.00490196078431373, 0.00980392156862745,
0.00735294117647059, 0.00980392156862745, 0.0245098039215686,
0, 0.0563725490196078, 0, 0, 0, 0, 0, 0, 0.053921568627451, 0,
0, 0.0220588235294118, 0, 0, 0, 0.00245098039215686, 0, 0, 0,
0, 0.00490196078431373, 0.00735294117647059, 0.00735294117647059,
0.0122549019607843, 0.0147058823529412, 0, 0.00980392156862745,
0, 0, 0.0196078431372549, 0, 0, 0.0122549019607843, 0, 0, 
0.0147058823529412,

0, 0.00490196078431373, 

Re: [R] Troubleshooting underidentification issues in structural equation modelling (SEM)

2013-02-15 Thread John Fox
Dear Ruijie and Bert,

I agree with Bert that it's very difficult to do effective statistical 
consulting long-distance by email. I think that you'd be much better served by 
getting competent statistical help locally, as I've already suggested. 

On the other hand, I'm surprised that your reading didn't suggest that a CFA 
model with latent variables that have just one observed indicator each, and in 
which both the factor loadings and error variances for these variables are free 
parameters, is underidentified. If you haven't already read it, I recommend 
Bollen's Structural Equations with Latent Variables (Wiley, 1989), despite its 
age.

Best,
 John

On Fri, 15 Feb 2013 02:11:01 -0800
 Bert Gunter gunter.ber...@gene.com wrote:
 These are statistical, not R issues, so please do not post further
 here. You are clearly out of your depth statistically. You need to get
 local statistical help, or you can try posting on a statistical list
 like stats.stackexchange.com if you care to take advice from unknown
 sources who don't understand the details of your situation.
 
 -- Bert
 
 On Fri, Feb 15, 2013 at 1:11 AM, Ruijie breakaw...@gmail.com wrote:
  Thanks Prof Fox for your guidance. My purpose in fitting this model is to
  contrast it with another model that I am proposing which I believe will be
  a better fit.
 
  On the point of some of the items being close to invariant, I had a close
  look at my data and indeed that is the case I am aware of it. However, I am
  not sure what to do with these items. Do I remove them? If I do, what
  threshold of variance do I set for removal? How do I decide on that
  threshold?
 
  I've combed a number of textbooks for answers but sadly have not found
  much. Hope you could offer some advice, thanks!
 
  Regards,
  Ruijie (RJ)
 
  
  He who has a why can endure any how.
 
  ~ Friedrich Nietzsche
 
 
  On 10 February 2013 00:38, John Fox j...@mcmaster.ca wrote:
 
  Dear Ruijie,
 
  Your model is underidentified by virtue of two of the factors having only
  one observed indicator each. No SEM software can magically estimate this
  model as it stands. Beyond that, I won't comment on the wisdom of what
  you're doing, such as computing covariances between ordinal variables --
  but
  see what I discovered below.
 
  Removing these two variables and the associated factors produces the
  following model:
 
  - snip 
 
   model - cfa(reference.indicators=FALSE)
  1: F01: I01, I02, I03
  2: F02: I04, I05, I06, I07, I08, I09, I10, I11, I12, I13
  3: F03: I14, I15, I16, I17, I18, I19, I20, I21, I22, I23, I24, I25, I26
  4: F04: I27, I28, I29, I30, I31, I32, I33, I34
  5: F05: I35, I36, I37, I38, I39, I40, I41, I42, I43
  6: F07: I46, I47, I48, I49, I50, I51
  7: F08: I54, I55, I56, I57, I58, I59, I60, I61, I62, I63, I64
  8: F09: I65, I66, I67
  9: F11: I69, I70, I71
  10:
  Read 9 items
  NOTE: adding 66 variances to the model
  
   cfa.output - sem(model, cov.mat, N = 900)
 
  - snip 
 
  sem() ran out of iterations, but the summary output is revealing:
 
  - snip 
 
   summary(cfa.output)
 
   Model Chisquare =  5677.1   Df =  2043 Pr(Chisq) = 0
   AIC =  6013.1
   BIC =  -8220.193
 
   Normalized Residuals
 Min. 1st Qu.  MedianMean 3rd Qu.Max.
  -3.9910 -0.5887 -0.1486  0.2588  0.8092 17.2900
 
   R-square for Endogenous Variables
  I01 I02 I03 I04 I05 I06 I07 I08 I09
  I10
   0.0953  0.1263  0.  0.1131  0.4039  0.2519  0.1168  0.0468  0.0005
  0.0059
  I11 I12 I13 I14 I15 I16 I17 I18 I19
  I20
   0.0479  0.0228  0.1150  0.2813  0.0001  0.0388  0.2106  0.0001  0.0913
  0.0063
  I21 I22 I23 I24 I25 I26 I27 I28 I29
  I30
   0.0041  0.0077  0.0022  0.  0.0299  0.0067  0.0019  0.0011  0.0010
  0.
  I31 I32 I33 I34 I35 I36 I37 I38 I39
  I40
   0.0005  0.0117  0.0270  0.0001  0.0084  0.0001  0.0256  0.4969  0.0613
  0.0515
  I41 I42 I43 I46 I47 I48 I49 I50 I51
  I54
   0.0005  0.0052  0.0307  0.0003  0.1131  0.0014  0.  0.1276  0.9728
  0.0520
  I55 I56 I57 I58 I59 I60 I61 I62 I63
  I64
   0.2930  0.0127  0.0543  0.0500  0.0378  0.0001  0.3048  0.0002  0.0304
  0.0001
  I65 I66 I67 I69 I70 I71
  56.7264  0.  0.0002  0.2220  0.2342  0.2240
 
   Parameter Estimates
   Estimate  Std Errorz value  Pr(|z|)
 
  lam[I01:F01]  3.023074e-02 5.133785e-03  5.888586224  3.895133e-09 I01 ---
  F01
  lam[I02:F01]  3.283192e-02 5.291069e-03  6.205157975  5.464199e-10 I02 ---
  F01
  lam[I03:F01]  1.123398e-04 2.695713e-03  0.041673509  9.667590e-01 I03 ---
  F01
  lam[I04:F02]  1.365329e-01 1.555023e-02  8.780124358  1.632940e-18 I04 ---
  F02
  lam[I05:F02]  9.525580e-02 5.517838e-03 17.263245517  8.896692e-67 I05 ---
  F02
  

Re: [R] Troubleshooting underidentification issues in structural equation modelling (SEM)

2013-02-15 Thread Ruijie
Much thanks to both for your advice. I'll take note to post relevant
questions going forward.

Regards,
Ruijie (RJ)


He who has a why can endure any how.

~ Friedrich Nietzsche


On 15 February 2013 21:12, John Fox j...@mcmaster.ca wrote:

 Dear Ruijie and Bert,

 I agree with Bert that it's very difficult to do effective statistical
 consulting long-distance by email. I think that you'd be much better served
 by getting competent statistical help locally, as I've already suggested.

 On the other hand, I'm surprised that your reading didn't suggest that a
 CFA model with latent variables that have just one observed indicator each,
 and in which both the factor loadings and error variances for these
 variables are free parameters, is underidentified. If you haven't already
 read it, I recommend Bollen's Structural Equations with Latent Variables
 (Wiley, 1989), despite its age.

 Best,
  John

 On Fri, 15 Feb 2013 02:11:01 -0800
  Bert Gunter gunter.ber...@gene.com wrote:
  These are statistical, not R issues, so please do not post further
  here. You are clearly out of your depth statistically. You need to get
  local statistical help, or you can try posting on a statistical list
  like stats.stackexchange.com if you care to take advice from unknown
  sources who don't understand the details of your situation.
 
  -- Bert
 
  On Fri, Feb 15, 2013 at 1:11 AM, Ruijie breakaw...@gmail.com wrote:
   Thanks Prof Fox for your guidance. My purpose in fitting this model is
 to
   contrast it with another model that I am proposing which I believe
 will be
   a better fit.
  
   On the point of some of the items being close to invariant, I had a
 close
   look at my data and indeed that is the case I am aware of it. However,
 I am
   not sure what to do with these items. Do I remove them? If I do, what
   threshold of variance do I set for removal? How do I decide on that
   threshold?
  
   I've combed a number of textbooks for answers but sadly have not found
   much. Hope you could offer some advice, thanks!
  
   Regards,
   Ruijie (RJ)
  
   
   He who has a why can endure any how.
  
   ~ Friedrich Nietzsche
  
  
   On 10 February 2013 00:38, John Fox j...@mcmaster.ca wrote:
  
   Dear Ruijie,
  
   Your model is underidentified by virtue of two of the factors having
 only
   one observed indicator each. No SEM software can magically estimate
 this
   model as it stands. Beyond that, I won't comment on the wisdom of what
   you're doing, such as computing covariances between ordinal variables
 --
   but
   see what I discovered below.
  
   Removing these two variables and the associated factors produces the
   following model:
  
   - snip 
  
model - cfa(reference.indicators=FALSE)
   1: F01: I01, I02, I03
   2: F02: I04, I05, I06, I07, I08, I09, I10, I11, I12, I13
   3: F03: I14, I15, I16, I17, I18, I19, I20, I21, I22, I23, I24, I25,
 I26
   4: F04: I27, I28, I29, I30, I31, I32, I33, I34
   5: F05: I35, I36, I37, I38, I39, I40, I41, I42, I43
   6: F07: I46, I47, I48, I49, I50, I51
   7: F08: I54, I55, I56, I57, I58, I59, I60, I61, I62, I63, I64
   8: F09: I65, I66, I67
   9: F11: I69, I70, I71
   10:
   Read 9 items
   NOTE: adding 66 variances to the model
   
cfa.output - sem(model, cov.mat, N = 900)
  
   - snip 
  
   sem() ran out of iterations, but the summary output is revealing:
  
   - snip 
  
summary(cfa.output)
  
Model Chisquare =  5677.1   Df =  2043 Pr(Chisq) = 0
AIC =  6013.1
BIC =  -8220.193
  
Normalized Residuals
  Min. 1st Qu.  MedianMean 3rd Qu.Max.
   -3.9910 -0.5887 -0.1486  0.2588  0.8092 17.2900
  
R-square for Endogenous Variables
   I01 I02 I03 I04 I05 I06 I07 I08
 I09
   I10
0.0953  0.1263  0.  0.1131  0.4039  0.2519  0.1168  0.0468
  0.0005
   0.0059
   I11 I12 I13 I14 I15 I16 I17 I18
 I19
   I20
0.0479  0.0228  0.1150  0.2813  0.0001  0.0388  0.2106  0.0001
  0.0913
   0.0063
   I21 I22 I23 I24 I25 I26 I27 I28
 I29
   I30
0.0041  0.0077  0.0022  0.  0.0299  0.0067  0.0019  0.0011
  0.0010
   0.
   I31 I32 I33 I34 I35 I36 I37 I38
 I39
   I40
0.0005  0.0117  0.0270  0.0001  0.0084  0.0001  0.0256  0.4969
  0.0613
   0.0515
   I41 I42 I43 I46 I47 I48 I49 I50
 I51
   I54
0.0005  0.0052  0.0307  0.0003  0.1131  0.0014  0.  0.1276
  0.9728
   0.0520
   I55 I56 I57 I58 I59 I60 I61 I62
 I63
   I64
0.2930  0.0127  0.0543  0.0500  0.0378  0.0001  0.3048  0.0002
  0.0304
   0.0001
   I65 I66 I67 I69 I70 I71
   56.7264  0.  0.0002  0.2220  0.2342  0.2240
  
Parameter Estimates
Estimate  Std Errorz value  Pr(|z|)
  
   lam[I01:F01]  3.023074e-02 5.133785e-03  5.888586224  

Re: [R] match a task No. to a few person-IDs

2013-02-15 Thread PIKAL Petr
Hi

?merge

Regards
Petr

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
 project.org] On Behalf Of Mat
 Sent: Friday, February 15, 2013 10:45 AM
 To: r-help@r-project.org
 Subject: [R] match a task No. to a few person-IDs
 
 hello together,
 
 i have a task No. in a data.frame
Task_No   Team
 A49397 1
 B49396 1
 
 and now i want to match my Person-IDs to each Task_No. My Person Ids
 look like this one:
 IDTeam
 A 50191
 B 472 1
 C 140 1
 D 50131
 
 The solution should be this:
 
Task_No   Team ID
 A49397 1 5019
 A49397 1 472
 A49397 1 140
 A49397 1 5013
 B49396 1 5019
 B49396 1472
 B49396 1140
 B49396 15013
 
 perhabs anyone can help me.
 
 Thank you
 
 Mat
 
 
 
 --
 View this message in context: http://r.789695.n4.nabble.com/match-a-
 task-No-to-a-few-person-IDs-tp4658646.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] plot custom x axis ticks values

2013-02-15 Thread PIKAL Petr
Hi

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
 project.org] On Behalf Of e-letter
 Sent: Friday, February 15, 2013 11:37 AM
 To: Duncan Murdoch
 Cc: r-help@r-project.org
 Subject: Re: [R] plot custom x axis ticks values
 
 On 15/02/2013, Duncan Murdoch murdoch.dun...@gmail.com wrote:
  On 13-02-15 3:28 AM, e-letter wrote:
  On 15/02/2013, Jim Lemon j...@bitwrit.com.au wrote:
  On 02/14/2013 09:41 PM, e-letter wrote:
  Readers,
 
  For this data set:
 
  testvalues-c(10,20,30,40)
 
  How to amend the plot instruction:
 
  plot(testvalues,ann=FALSE,type='l',yaxt='n',xaxt='n')
 
  so that x axis ticks labels can be added to existing graph with
  arbitrary value such as 0,100,200,300)?
 
  Hi r2151,
  If you want the labels to fit on the axis you will have to include
  this information in the call to plot:
 
  plot(testvalues,ann=FALSE,type='l',yaxt='n',xaxt='n',ylim=c(0,300))
  axis(2,at=c(0,100,200,300))
 
 
  plot(testvalues,ann=FALSE,type='l',yaxt='n',xaxt='n',xlim=c(0,300))
  axis(1,at=c(0,100,200,300))
 
  The use of 'xlim' makes the graph unacceptable. Instead, it would be
  better to plot:
 
  plot(testvalues,ann=FALSE,type='l',yaxt='n')
 
  The x-axis shows tick marks '1.0', '1.5', ...
 
  It is required to replace these values with custom values (e.g. 10,
 20...
  etc.)
 
  Why not rescale the values before plotting, and use the automatic
 ticks?
You can lie about the user scale, but it doesn't always give a
  helpful plot, e.g.
 
  axis(1, at=1:4, labels=c(100, 300, 400, 200))
 
 
 Tried, but realised that this works when the vectors for 'at' and
 'labels' are equal size. For a data set of e.g. 10 values, syntax error
 returns.

It is rather twisted idea to expect that you say to R:

Please, put 10 labels to x axis at 4 positions.

and expect that the program will not complain.

Maybe you shall consult axis help page, which clearly says:

at - the points at which tick-marks are to be drawn. Non-finite (infinite, NaN 
or NA) values are omitted. By default (when NULL) tickmark locations are 
computed, see 'Details' below.
 
labels  - this can either be a logical value specifying whether (numerical) 
annotations are to be made at the tickmarks, or a character or expression 
vector of labels to be placed at the tickpoints. (Other objects are coerced by 
as.graphicsAnnot.) If this is not logical, at should also be supplied and of 
the same length. If labels is of length zero after coercion, it has the same 
effect as supplying TRUE.

Regards
Petr



 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] match a task No. to a few person-IDs

2013-02-15 Thread arun
Hi,
Try this:
dat1- read.table(text=
  Task_No  Team   
A    49397    1 
B    49396    1
,sep=,header=TRUE)

dat2- read.table(text=
    ID    Team
A 5019    1
B 472    1
C 140    1
D 5013    1
,sep=,header=TRUE)
 res1- as.matrix(merge(dat1,dat2,by=Team))
 row.names(res1)-LETTERS[match(res1[,2],dat1[,1])]  


library(plyr)
res2- as.matrix(join(dat1,dat2,by=Team))
 row.names(res2)-LETTERS[match(res2[,1],dat1[,1])]
 res2
#  Task_No Team   ID
#A   49397    1 5019
#A   49397    1  472
#A   49397    1  140
#A   49397    1 5013
#B   49396    1 5019
#B   49396    1  472
#B   49396    1  140
#B   49396    1 5013
A.K.




- Original Message -
From: Mat matthias.we...@fnt.de
To: r-help@r-project.org
Cc: 
Sent: Friday, February 15, 2013 4:45 AM
Subject: [R] match a task No. to a few person-IDs

hello together,

i have a task No. in a data.frame
       Task_No       Team    
A    49397         1          
B    49396         1      

and now i want to match my Person-IDs to each Task_No. My Person Ids look
like this one:
    ID        Team
A 5019    1
B 472    1
C 140    1
D 5013    1

The solution should be this:

       Task_No       Team     ID
A    49397         1               5019
A    49397         1               472
A    49397         1               140
A    49397         1               5013
B    49396         1         5019
B    49396         1            472
B    49396         1            140
B    49396         1            5013

perhabs anyone can help me.

Thank you

Mat



--
View this message in context: 
http://r.789695.n4.nabble.com/match-a-task-No-to-a-few-person-IDs-tp4658646.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] Why no line? (ex. from Andy Filed book)

2013-02-15 Thread PIKAL Petr
Hi

The code gives me some warnings but all seems to be OK. 4 points, dashed red 
line and error bars. Did you expect some other line?

maybe issue of version/OS? (unstated)
Petr

 sessionInfo()
R Under development (unstable) (2012-10-29 r61044)
Platform: i386-w64-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=Czech_Czech Republic.1250  LC_CTYPE=Czech_Czech Republic.1250   
[3] LC_MONETARY=Czech_Czech Republic.1250 LC_NUMERIC=C 
[5] LC_TIME=Czech_Czech Republic.1250

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

other attached packages:
[1] Hmisc_3.9-2  survival_2.36-14 reshape_0.8.4plyr_1.7.1  
[5] ggplot2_0.9.0lattice_0.20-10  fun_1.0 

loaded via a namespace (and not attached):
 [1] cluster_1.14.3 colorspace_1.1-1   digest_0.5.1   dichromat_1.2-4   
 [5] grid_2.16.0MASS_7.3-22memoise_0.1munsell_0.3   
 [9] nlme_3.1-105   proto_0.3-9.2  RColorBrewer_1.0-5 reshape2_1.2.1
[13] scales_0.2.0   stringr_0.6tools_2.16.0  







 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
 project.org] On Behalf Of Kevin Mc Inerney
 Sent: Friday, February 15, 2013 1:32 AM
 To: r-help@r-project.org
 Subject: [R] Why no line? (ex. from Andy Filed book)
 
 The following script is written by the author of a book on R--Andy
 Field:
 You can also download the small datafile, hiccups.dat, from this
 address:
 
 http://www.sagepub.com/dsur/study/articles.htm
 
 The script:
 
 hiccupsData - read.delim(Hiccups.dat,  header = TRUE)
 hiccups-stack(hiccupsData)
 names(hiccups)-c(Hiccups,Intervention)
 hiccups$Intervention_Factor-factor(hiccups$Intervention, levels =
 hiccups$Intervention)
 
 line - ggplot(hiccups, aes(Intervention_Factor, Hiccups)) line +
 stat_summary(fun.y = mean, geom = point) + stat_summary(fun.y = mean,
 geom = line, aes(group=1),colour = Red, linetype = dashed)+
 stat_summary(fun.data = mean_cl_boot, geom = errorbar, width = 0.2) +
 labs(x = Intervention, y = Mean Number of Hiccups)
 saveInImageDirectory(04 Hiccups Line.png)
 
 Why does this script not give me a line?
 
 It's driving e crazy :-|
 
   [[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] Making the plot window wider and using the predict function

2013-02-15 Thread Rasmus Hedegaard




Hello,

I am new to R and have a couple of questions. My data set contains the 
variables Bwt and Hwt, which are bodyweight and heartweight, respectively, 
of a group of cats.
With the following code, I am making two plots, both to be viewed in the same 
plot window in R:

library(MASS)
 maleData - subset(cats, Sex == M)
 linreg0 - lm(maleData$Hwt ~ maleData$Bwt)
 par(mfrow=c(1,2))
 plot(maleData$Hwt ~ maleData$Bwt)
 plot(rstandard(linreg0) ~ fitted(linreg0))

My problem is that the two plots end up all oblong and squashed, because the 
plot window doesn't have a size suited for having two plots next to one 
another. Can I tell R to adjust the plot window?

Also, I would like to do something along the lines of:

newData - data.frame(Bwt = 3.5)
predict(linreg0,newData,interval=p)

This doesn't work - my guess is that to R, Bwt is not a variable in maleData, 
but maleData$Bwt is. I could use an attach command, but is it possible to get 
this to work without doing so?

Kind regards,
Rasmus Hedegaard. 
[[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 deal with zero-inflated proportion data?

2013-02-15 Thread Stefan Sobernig

Dear List,

In my previous posting 
(https://stat.ethz.ch/pipermail/r-help/2013-February/347593.html), I 
refer to playing around with distribution fitting for a particular sort 
of data (see dat1 in the previous posting):


I collected absolute node degrees of some network structure and computed 
the relative node degrees. This is because I want to compare different 
networks at another stage.


Based on relative degrees, I attempted some distribution fitting 
procedures (see my previous posting). However, I realized that zero 
values (representing isolates in the network) in an otherwise continuous 
variable are problematic (well, not much of a surprise).


I am wondering how to best deal with those isolates/zeros? Excluding 
them entirely would be an option but isolates are relevant for my 
analysis. Are there best practises for dealing with such a hybrid variable?


I am aware of censored methods (for regression analysis etc.), but not 
when it comes down to distribution fitting, I fear.


Many thanks,
Stefan

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


Re: [R] Making the plot window wider and using the predict function

2013-02-15 Thread PIKAL Petr
Hi
 
 
 Hello,
 
 I am new to R and have a couple of questions. My data set contains the
 variables Bwt and Hwt, which are bodyweight and heartweight,
 respectively, of a group of cats.
 With the following code, I am making two plots, both to be viewed in
 the same plot window in R:
 
 library(MASS)
  maleData - subset(cats, Sex == M)
  linreg0 - lm(maleData$Hwt ~ maleData$Bwt)

The second question is answered by:

linreg0 - lm(Hwt ~ Bwt, data=maleData)

after that your predict code works.

  par(mfrow=c(1,2))
  plot(maleData$Hwt ~ maleData$Bwt)
  plot(rstandard(linreg0) ~ fitted(linreg0))
 
 My problem is that the two plots end up all oblong and squashed,
 because the plot window doesn't have a size suited for having two plots
 next to one another. Can I tell R to adjust the plot window?

Can you use mouse for resizing?

Regards 
Petr

 
 Also, I would like to do something along the lines of:
 
 newData - data.frame(Bwt = 3.5)
 predict(linreg0,newData,interval=p)
 
 This doesn't work - my guess is that to R, Bwt is not a variable in
 maleData, but maleData$Bwt is. I could use an attach command, but is it
 possible to get this to work without doing so?
 
 Kind regards,
 Rasmus Hedegaard.
   [[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] Remove site path from .libPaths

2013-02-15 Thread Rainer M Krug
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Hi

I am sure I am missing something really basic, but I can't figure it out.

I want to start R so that I can specify the location for the Library tree. In 
principel simple:

As I only want it dependent on the directory I stat R in, I put a .Rprofile 
file in the directory.

My default path is:

 .libPaths()
[1] /home/rkrug/R/i686-pc-linux-gnu-library/2.15
[2] /usr/lib/R/site-library
[3] /usr/lib/R/library
 

But I want to have it

 .libPaths()
[1] /home/rkrug/THE_DIRECTORY/library
[2] /usr/lib/R/library
 

The first part is easy:

.libPaths(/home/rkrug/THE_DIRECTORY/library)

but how can I remove the site library?

If I set

In R 15.2 under Ubuntu, started with

R --vanilla

 .libPaths()
[1] /home/rkrug/R/i686-pc-linux-gnu-library/2.15
[2] /usr/lib/R/site-library
[3] /usr/lib/R/library
 .Library.site
[1] /usr/lib/R/site-library /usr/lib/R/library
 .Library
[1] /usr/lib/R/library
 .Library.site -  .Library.site
[1] 
 .libPaths()
[1] /home/rkrug/R/i686-pc-linux-gnu-library/2.15
[2] /usr/lib/R/site-library
[3] /usr/lib/R/library
 .libPaths() .libPaths()
[1] /usr/lib/R/site-library /usr/lib/R/library
 .Library
[1] /usr/lib/R/library
 .Library.site
[1] 
 


 .libPaths() .libPaths()
[1] /usr/lib/R/site-library /usr/lib/R/library
 

even executing

 .libPaths(.libPaths())

does not change anything.

Am I missing something or is there a bug in .libPaths()?

Cheers,

Rainer



- -- 
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology, 
UCT), Dipl. Phys.
(Germany)

Centre of Excellence for Invasion Biology
Stellenbosch University
South Africa

Tel :   +33 - (0)9 53 10 27 44
Cell:   +33 - (0)6 85 62 59 98
Fax :   +33 - (0)9 58 10 27 44

Fax (D):+49 - (0)3 21 21 25 22 44

email:  rai...@krugs.de

Skype:  RMkrug
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.11 (GNU/Linux)
Comment: Using GnuPG with Thunderbird - http://www.enigmail.net/

iQEcBAEBAgAGBQJRHkX0AAoJENvXNx4PUvmCgDYH/jwOadmO1cpPAK2yNaCfldjm
gzQ4wLf6o846xM8UQj78rIczvlycLqJP84jZe8kf8iHItFKSQbekpzujazIVEpcV
TCcaiLbm9ex5qpnRVAH8VOph59acsnEjo67dztmBl+uPb30l4hmG/dwgUeJE7HrP
oy+jtJZVostGzSiKV8wXv9C4Ohu9WWQ3Rw7VEGRKsAWcDfHGl79isEA5+ONbJjdz
32erpXBbCquJPFV4OGT7DqI8D/xh8dTVNK9ew/cdUjXlCysGAHUEQXRX2z+oxcJ5
OIwNzCADjgM1/yIc5hCtTZtFIvFWo169vTgM2/0UDnpneJLHCrHFvuU94eo/6cI=
=pJ7s
-END PGP SIGNATURE-

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


Re: [R] match a task No. to a few person-IDs

2013-02-15 Thread Rui Barradas

Hello,

Use ?merge.


tasks - read.table(text = 
Task_No   Team
A49397   1  
B49396   1  
, header = TRUE)

ids - read.table(text = 
IDTeam
A 5019  1
B 472  1
C 140  1
D 5013  1
, header = TRUE)


merge(tasks, ids, all.x = TRUE)


Hope this helps,

Rui Barradas

Em 15-02-2013 09:45, Mat escreveu:

hello together,

i have a task No. in a data.frame
Task_No   Team
A49397   1  
B49396   1  

and now i want to match my Person-IDs to each Task_No. My Person Ids look
like this one:
 IDTeam
A 5019  1
B 472   1
C 140   1
D 5013  1

The solution should be this:

Task_No   Team ID
A49397   1 5019
A49397   1 472
A49397   1 140
A49397   1 5013
B49396   1 5019
B49396   1472
B49396   1140
B49396   15013

perhabs anyone can help me.

Thank you

Mat



--
View this message in context: 
http://r.789695.n4.nabble.com/match-a-task-No-to-a-few-person-IDs-tp4658646.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] Remove site path from .libPaths

2013-02-15 Thread Duncan Murdoch

On 13-02-15 9:28 AM, Rainer M Krug wrote:

-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Hi

I am sure I am missing something really basic, but I can't figure it out.

I want to start R so that I can specify the location for the Library tree. In 
principel simple:

As I only want it dependent on the directory I stat R in, I put a .Rprofile 
file in the directory.

My default path is:


.libPaths()

[1] /home/rkrug/R/i686-pc-linux-gnu-library/2.15
[2] /usr/lib/R/site-library
[3] /usr/lib/R/library




But I want to have it


.libPaths()

[1] /home/rkrug/THE_DIRECTORY/library
[2] /usr/lib/R/library




The first part is easy:

.libPaths(/home/rkrug/THE_DIRECTORY/library)

but how can I remove the site library?

If I set

In R 15.2 under Ubuntu, started with

R --vanilla


.libPaths()

[1] /home/rkrug/R/i686-pc-linux-gnu-library/2.15
[2] /usr/lib/R/site-library
[3] /usr/lib/R/library

.Library.site

[1] /usr/lib/R/site-library /usr/lib/R/library

.Library

[1] /usr/lib/R/library

.Library.site -  .Library.site

[1] 

.libPaths()

[1] /home/rkrug/R/i686-pc-linux-gnu-library/2.15
[2] /usr/lib/R/site-library
[3] /usr/lib/R/library

.libPaths() .libPaths()

[1] /usr/lib/R/site-library /usr/lib/R/library

.Library

[1] /usr/lib/R/library

.Library.site

[1] 






.libPaths() .libPaths()

[1] /usr/lib/R/site-library /usr/lib/R/library




even executing


.libPaths(.libPaths())


does not change anything.

Am I missing something or is there a bug in .libPaths()?


I don't think there is a bug.  As ?.libPaths says,

Function .libPaths always uses the values of .Library and .Library.site 
in the base namespace. .Library.site can be set by the site in 
‘Rprofile.site’, which should be followed by a call to 
.libPaths(.libPaths()) to make use of the updated value.


So you can't change the value of .Library.site after starting R.  You'll 
need to use environment variables (as described on that page) to do it.


Duncan Murdoch



Cheers,

Rainer



- --
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology, 
UCT), Dipl. Phys.
(Germany)

Centre of Excellence for Invasion Biology
Stellenbosch University
South Africa

Tel :   +33 - (0)9 53 10 27 44
Cell:   +33 - (0)6 85 62 59 98
Fax :   +33 - (0)9 58 10 27 44

Fax (D):+49 - (0)3 21 21 25 22 44

email:  rai...@krugs.de

Skype:  RMkrug
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.11 (GNU/Linux)
Comment: Using GnuPG with Thunderbird - http://www.enigmail.net/

iQEcBAEBAgAGBQJRHkX0AAoJENvXNx4PUvmCgDYH/jwOadmO1cpPAK2yNaCfldjm
gzQ4wLf6o846xM8UQj78rIczvlycLqJP84jZe8kf8iHItFKSQbekpzujazIVEpcV
TCcaiLbm9ex5qpnRVAH8VOph59acsnEjo67dztmBl+uPb30l4hmG/dwgUeJE7HrP
oy+jtJZVostGzSiKV8wXv9C4Ohu9WWQ3Rw7VEGRKsAWcDfHGl79isEA5+ONbJjdz
32erpXBbCquJPFV4OGT7DqI8D/xh8dTVNK9ew/cdUjXlCysGAHUEQXRX2z+oxcJ5
OIwNzCADjgM1/yIc5hCtTZtFIvFWo169vTgM2/0UDnpneJLHCrHFvuU94eo/6cI=
=pJ7s
-END PGP SIGNATURE-

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

2013-02-15 Thread Rainer M Krug
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

On 15/02/13 15:41, Duncan Murdoch wrote:
 On 13-02-15 9:28 AM, Rainer M Krug wrote: Hi
 
 I am sure I am missing something really basic, but I can't figure it out.
 
 I want to start R so that I can specify the location for the Library tree. In 
 principel
 simple:
 
 As I only want it dependent on the directory I stat R in, I put a .Rprofile 
 file in the
 directory.
 
 My default path is:
 
 .libPaths()
 [1] /home/rkrug/R/i686-pc-linux-gnu-library/2.15 [2] 
 /usr/lib/R/site-library [3]
 /usr/lib/R/library
 
 
 But I want to have it
 
 .libPaths()
 [1] /home/rkrug/THE_DIRECTORY/library [2] /usr/lib/R/library
 
 
 The first part is easy:
 
 .libPaths(/home/rkrug/THE_DIRECTORY/library)
 
 but how can I remove the site library?
 
 If I set
 
 In R 15.2 under Ubuntu, started with
 
 R --vanilla
 
 .libPaths()
 [1] /home/rkrug/R/i686-pc-linux-gnu-library/2.15 [2] 
 /usr/lib/R/site-library [3]
 /usr/lib/R/library
 .Library.site
 [1] /usr/lib/R/site-library /usr/lib/R/library
 .Library
 [1] /usr/lib/R/library
 .Library.site -  .Library.site
 [1] 
 .libPaths()
 [1] /home/rkrug/R/i686-pc-linux-gnu-library/2.15 [2] 
 /usr/lib/R/site-library [3]
 /usr/lib/R/library
 .libPaths() .libPaths()
 [1] /usr/lib/R/site-library /usr/lib/R/library
 .Library
 [1] /usr/lib/R/library
 .Library.site
 [1] 
 
 
 
 .libPaths() .libPaths()
 [1] /usr/lib/R/site-library /usr/lib/R/library
 
 
 even executing
 
 .libPaths(.libPaths())
 
 does not change anything.
 
 Am I missing something or is there a bug in .libPaths()?
 
 I don't think there is a bug.  As ?.libPaths says,
 
 Function .libPaths always uses the values of .Library and .Library.site in 
 the base
 namespace. .Library.site can be set by the site in ‘Rprofile.site’, which 
 should be followed
 by a call to .libPaths(.libPaths()) to make use of the updated value.
 
 So you can't change the value of .Library.site after starting R.  You'll 
 need to use
 environment variables (as described on that page) to do it.


Thanks for the clarification - the use of the word site in set by the site 
in ‘Rprofile.site’ 
was not clear to me. I thik it would be much clearer, if it states that the 
variables .Library and
.Library.site can not be changed while R is running although it looks as if 
they can be changed
(but then new ones in the top environment are created).

Now just one more question: where can I see that

.Library - /A/New/Existing/Directory

has no impact?

I expected to see this in the code of .libPaths, but I didn't:

 .libPaths
function (new)
{
if (!missing(new)) {
new - Sys.glob(path.expand(new))
paths - unique(normalizePath(c(new, .Library.site, .Library),
/))
.lib.loc - paths[file.info(paths)$isdir %in% TRUE]
}
else .lib.loc
}
bytecode: 0x88e7640
environment: 0x88e685c
 

Thanks,

Rainer



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

- -- 
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology, 
UCT), Dipl. Phys.
(Germany)

Centre of Excellence for Invasion Biology
Stellenbosch University
South Africa

Tel :   +33 - (0)9 53 10 27 44
Cell:   +33 - (0)6 85 62 59 98
Fax :   +33 - (0)9 58 10 27 44

Fax (D):+49 - (0)3 21 21 25 22 44

email:  rai...@krugs.de

Skype:  RMkrug
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.11 (GNU/Linux)
Comment: Using GnuPG with Thunderbird - http://www.enigmail.net/

iQEcBAEBAgAGBQJRHkwmAAoJENvXNx4PUvmCnEcIAKPC2wxPWSX5Y/h5jfCg5b9g
OKsGsmL/0bu7MdfoMFKMdFL+TPQ1werGEgcKlmwUuZ/QGzd21BkN+C/Dl5gXSYRB
gHCGrQWUSemUlt+V2BTW1IaoksLhbj+QQBVNJ767eO5LkoIFujIwrCZ8o7upn/cy
6W3EhP+utHWz2C3+/gUR/c5FqPM6r4+9JiNVCxRzjCU/8/wip6dmoQRpZGXgAb8d
Koa5lRJlMUV0IMItMAoIIKzr1IkXaZ91VcTI1X13UTvzT7R/v5xgOnfJPEb5DRqI
clFGa8VaWx2RfmU6ihlmj8k8q4lOyI1hNXE6q8wX7aI3ERW4gSFCuDKASshyZXM=
=XGla
-END PGP SIGNATURE-

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


Re: [R] Remove site path from .libPaths

2013-02-15 Thread Duncan Murdoch

On 13-02-15 9:54 AM, Rainer M Krug wrote:

-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

On 15/02/13 15:41, Duncan Murdoch wrote:

On 13-02-15 9:28 AM, Rainer M Krug wrote: Hi

I am sure I am missing something really basic, but I can't figure it out.

I want to start R so that I can specify the location for the Library tree. In 
principel
simple:

As I only want it dependent on the directory I stat R in, I put a .Rprofile 
file in the
directory.

My default path is:


.libPaths()

[1] /home/rkrug/R/i686-pc-linux-gnu-library/2.15 [2] 
/usr/lib/R/site-library [3]
/usr/lib/R/library




But I want to have it


.libPaths()

[1] /home/rkrug/THE_DIRECTORY/library [2] /usr/lib/R/library




The first part is easy:

.libPaths(/home/rkrug/THE_DIRECTORY/library)

but how can I remove the site library?

If I set

In R 15.2 under Ubuntu, started with

R --vanilla


.libPaths()

[1] /home/rkrug/R/i686-pc-linux-gnu-library/2.15 [2] 
/usr/lib/R/site-library [3]
/usr/lib/R/library

.Library.site

[1] /usr/lib/R/site-library /usr/lib/R/library

.Library

[1] /usr/lib/R/library

.Library.site -  .Library.site

[1] 

.libPaths()

[1] /home/rkrug/R/i686-pc-linux-gnu-library/2.15 [2] 
/usr/lib/R/site-library [3]
/usr/lib/R/library

.libPaths() .libPaths()

[1] /usr/lib/R/site-library /usr/lib/R/library

.Library

[1] /usr/lib/R/library

.Library.site

[1] 






.libPaths() .libPaths()

[1] /usr/lib/R/site-library /usr/lib/R/library




even executing


.libPaths(.libPaths())


does not change anything.

Am I missing something or is there a bug in .libPaths()?


I don't think there is a bug.  As ?.libPaths says,



Function .libPaths always uses the values of .Library and .Library.site in the 
base
namespace. .Library.site can be set by the site in ‘Rprofile.site’, which 
should be followed
by a call to .libPaths(.libPaths()) to make use of the updated value.



So you can't change the value of .Library.site after starting R.  You'll need 
to use
environment variables (as described on that page) to do it.



Thanks for the clarification - the use of the word site in set by the site in 
‘Rprofile.site’ 
was not clear to me. I thik it would be much clearer, if it states that the 
variables .Library and
.Library.site can not be changed while R is running although it looks as if 
they can be changed
(but then new ones in the top environment are created).

Now just one more question: where can I see that

.Library - /A/New/Existing/Directory

has no impact?


If you look at environment(.libPaths) you'll see that it is not 
.GlobalEnv, so the .lib.loc assignment might not happen where you think. 
 Indeed,


ls(environment(.libPaths), all=TRUE)

will show that it contains one object, .lib.loc.  You can modify that 
variable, but it's quite a risky (in the sense of being undocumented and 
unsupported) thing to do.


Duncan Murdoch



I expected to see this in the code of .libPaths, but I didn't:


.libPaths

function (new)
{
 if (!missing(new)) {
 new - Sys.glob(path.expand(new))
 paths - unique(normalizePath(c(new, .Library.site, .Library),
 /))
 .lib.loc - paths[file.info(paths)$isdir %in% TRUE]
 }
 else .lib.loc
}
bytecode: 0x88e7640
environment: 0x88e685c




Thanks,

Rainer






Duncan Murdoch



Cheers,

Rainer





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





- --
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology, 
UCT), Dipl. Phys.
(Germany)

Centre of Excellence for Invasion Biology
Stellenbosch University
South Africa

Tel :   +33 - (0)9 53 10 27 44
Cell:   +33 - (0)6 85 62 59 98
Fax :   +33 - (0)9 58 10 27 44

Fax (D):+49 - (0)3 21 21 25 22 44

email:  rai...@krugs.de

Skype:  RMkrug
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.11 (GNU/Linux)
Comment: Using GnuPG with Thunderbird - http://www.enigmail.net/

iQEcBAEBAgAGBQJRHkwmAAoJENvXNx4PUvmCnEcIAKPC2wxPWSX5Y/h5jfCg5b9g
OKsGsmL/0bu7MdfoMFKMdFL+TPQ1werGEgcKlmwUuZ/QGzd21BkN+C/Dl5gXSYRB
gHCGrQWUSemUlt+V2BTW1IaoksLhbj+QQBVNJ767eO5LkoIFujIwrCZ8o7upn/cy
6W3EhP+utHWz2C3+/gUR/c5FqPM6r4+9JiNVCxRzjCU/8/wip6dmoQRpZGXgAb8d
Koa5lRJlMUV0IMItMAoIIKzr1IkXaZ91VcTI1X13UTvzT7R/v5xgOnfJPEb5DRqI
clFGa8VaWx2RfmU6ihlmj8k8q4lOyI1hNXE6q8wX7aI3ERW4gSFCuDKASshyZXM=
=XGla
-END PGP SIGNATURE-



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 stack row vector on top of each other?

2013-02-15 Thread Giovanni Petris
How about 

c(a, b)

?

HTH,
Giovanni

From: r-help-boun...@r-project.org [r-help-boun...@r-project.org] on behalf of 
Nordlund, Dan (DSHS/RDA) [nord...@dshs.wa.gov]
Sent: Thursday, February 14, 2013 7:19 PM
To: r-help
Subject: Re: [R] How to stack row vector on top of each other?

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
 project.org] On Behalf Of C W
 Sent: Thursday, February 14, 2013 5:08 PM
 To: r-help
 Subject: [R] How to stack row vector on top of each other?

 Hi list,
 How do you actually stack a vector on top of each other?  Say, I want
 everything in a row vector.  Neither rbind(), nor cbind() will do the
 job.
  It gives me 2 dimension.

 Here's my reproducible example:
  a - rnorm(10)
  b - rnorm(10)
  c - cbind(a,b)
  dim(c)
 [1] 10  2

  d - rbind(a,b)
  dim(d)
 [1]  2 10

 Thanks,
 Mike

I guess I don't know what you mean by actually stack a vector on top of each 
other.  Given the vectors

a - 1:3
b - 4:6

What result do you want from stacking a and b?


Dan

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


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

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


Re: [R] Remove site path from .libPaths

2013-02-15 Thread Rainer M Krug
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

On 15/02/13 16:06, Duncan Murdoch wrote:
 On 13-02-15 9:54 AM, Rainer M Krug wrote: On 15/02/13 15:41, Duncan Murdoch 
 wrote:
 On 13-02-15 9:28 AM, Rainer M Krug wrote: Hi
 
 I am sure I am missing something really basic, but I can't figure it out.
 
 I want to start R so that I can specify the location for the Library tree. 
 In principel 
 simple:
 
 As I only want it dependent on the directory I stat R in, I put a 
 .Rprofile file in the 
 directory.
 
 My default path is:
 
 .libPaths()
 [1] /home/rkrug/R/i686-pc-linux-gnu-library/2.15 [2] 
 /usr/lib/R/site-library [3] 
 /usr/lib/R/library
 
 
 But I want to have it
 
 .libPaths()
 [1] /home/rkrug/THE_DIRECTORY/library [2] /usr/lib/R/library
 
 
 The first part is easy:
 
 .libPaths(/home/rkrug/THE_DIRECTORY/library)
 
 but how can I remove the site library?
 
 If I set
 
 In R 15.2 under Ubuntu, started with
 
 R --vanilla
 
 .libPaths()
 [1] /home/rkrug/R/i686-pc-linux-gnu-library/2.15 [2] 
 /usr/lib/R/site-library [3] 
 /usr/lib/R/library
 .Library.site
 [1] /usr/lib/R/site-library /usr/lib/R/library
 .Library
 [1] /usr/lib/R/library
 .Library.site -  .Library.site
 [1] 
 .libPaths()
 [1] /home/rkrug/R/i686-pc-linux-gnu-library/2.15 [2] 
 /usr/lib/R/site-library [3] 
 /usr/lib/R/library
 .libPaths() .libPaths()
 [1] /usr/lib/R/site-library /usr/lib/R/library
 .Library
 [1] /usr/lib/R/library
 .Library.site
 [1] 
 
 
 
 .libPaths() .libPaths()
 [1] /usr/lib/R/site-library /usr/lib/R/library
 
 
 even executing
 
 .libPaths(.libPaths())
 
 does not change anything.
 
 Am I missing something or is there a bug in .libPaths()?
 
 I don't think there is a bug.  As ?.libPaths says,
 
 Function .libPaths always uses the values of .Library and .Library.site 
 in the base 
 namespace. .Library.site can be set by the site in ‘Rprofile.site’, which 
 should be
 followed by a call to .libPaths(.libPaths()) to make use of the updated 
 value.
 
 So you can't change the value of .Library.site after starting R.  You'll 
 need to use 
 environment variables (as described on that page) to do it.
 
 
 Thanks for the clarification - the use of the word site in set by the site 
 in
 ‘Rprofile.site’  was not clear to me. I thik it would be much clearer, if it 
 states that the
 variables .Library and .Library.site can not be changed while R is running 
 although it looks as
 if they can be changed (but then new ones in the top environment are created).
 
 Now just one more question: where can I see that
 
 .Library - /A/New/Existing/Directory
 
 has no impact?
 
 If you look at environment(.libPaths) you'll see that it is not .GlobalEnv, 
 so the .lib.loc 
 assignment might not happen where you think.  Indeed,
 
 ls(environment(.libPaths), all=TRUE)
 
 will show that it contains one object, .lib.loc.  You can modify that 
 variable, but it's
 quite a risky (in the sense of being undocumented and unsupported) thing to 
 do.

Ok, thanks for the clarification, and agreed that it is not advisable. I will 
then go a different
route and start R with given Environmental variables.
Thanks  lot,

Rainer


 
 Duncan Murdoch
 
 
 I expected to see this in the code of .libPaths, but I didn't:
 
 .libPaths
 function (new) { if (!missing(new)) { new - Sys.glob(path.expand(new)) paths 
 -
 unique(normalizePath(c(new, .Library.site, .Library), /)) .lib.loc -
 paths[file.info(paths)$isdir %in% TRUE] } else .lib.loc } bytecode: 
 0x88e7640 environment:
 0x88e685c
 
 
 Thanks,
 
 Rainer
 
 
 
 
 Duncan Murdoch
 
 
 Cheers,
 
 Rainer
 
 
 
 
 __ R-help@r-project.org 
 mailing list 
 https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting 
 guide 
 http://www.R-project.org/posting-guide.html and provide commented, 
 minimal,
 self-contained, reproducible code.
 
 
 
 
 

-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.11 (GNU/Linux)
Comment: Using GnuPG with Thunderbird - http://www.enigmail.net/

iQEcBAEBAgAGBQJRHlgpAAoJENvXNx4PUvmCel0H/jJ0KYNfEryzjnSIt9U9xrk/
cEgnFLRAzoWuBaMvTuPVthiOHIjIPWyiauoZtn9Vju6HALBQAY9h7SWpwMyVIm6y
VwHsRZ3Szbn1kRPkIbHpZFlMHtfFg9I9KlSywBPHjD1FVEUI8TXJSz4padjbmSWc
VHRK88o2YHNAA6ss7AyGjpXuFxwIvCNjxTVitPdaMBqnmoJtO3jIOgcPAp18pLlK
uYv1pvcVKY6yizNwwneMhIO0dqmZMIjL/TT6Z8IT2K7EUOqEtWXQ2DhSFOFBtYGi
olae3hK9I8FI+wS5oD7i86PUNPB8KDya3KnGk1USEKF69nesRbWXU3fuC6BaI+A=
=t0Hg
-END PGP SIGNATURE-

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

2013-02-15 Thread C Lin

Can anyone help explain to me why the two codes below have different result? I 
thought I can use log(time)~. to replace log(time)~dist+climb+timef.I am using 
CVlm from DAAG package. I think nihills is preloaded with the package. Thanks 
in advance.
 CVlm(df=nihills, form.lm=formula(log(time)~.),plotit=Observed,m=2)Analysis 
 of Variance Table
Response: log(time)  Df Sum Sq Mean Sq F value  Pr(F)dist   1  
 6.346.34  384.31 4.6e-14 ***climb  1   0.120.127.24  0.0145 *  
timef  1   0.190.19   11.29  0.0033 ** Residuals 19   0.310.02  
  ---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 


fold 1 Observations in test set: 11 Slieve Gullion McVeigh Classic 
Tollymore Mountain Moughanmore Hen  Cock Annalong Horseshoe  Rocky Meelbeg 
Meelmore Slieve Donard Seven Sevens Slieve GallionPredicted   -0.680
  -0.538 -0.610  -0.786 -0.849   7.87e-01 
-0.682-6.20e-01 -2.98e-01 1.41e+00  -5.87e-01cvpred 
 -0.762  -0.614 -0.727  -0.769 -0.799   
6.67e-01 -0.648-7.89e-01 -5.30e-02 1.36e+00  
-4.52e-01log(time)   -0.762  -0.614 -0.727  
-0.769 -0.799   6.67e-01 -0.648-7.89e-01 -5.30e-02 
1.36e+00  -4.52e-01CV residual  0.000   0.000  
0.000   0.000  0.000   1.11e-16  0.000 1.11e-16 
-4.86e-17 4.44e-16  -5.55e-17
Sum of squares = 0Mean square = 0n = 11 
fold 2 Observations in test set: 12 Binevenagh Glenariff Mountain 
Donard  Commedagh Slieve Martin Monument Race Loughshannagh Horseshoe Donard 
Forest Flagstaff to Carling Slieve Bearnagh Lurig Challenge Scrabo Hill Race 
BARF Turkey TrotPredicted-1.78e-01  -4.89e-01   2.78e-02
-0.577 -7.27e-01  -0.623-0.571 
3.03e-01  -0.401   -7.20e-01   -0.889
-4.88e-01cvpred   -1.53e-01  -3.52e-01   3.79e-02
-0.597 -7.51e-01  -0.435-0.657 3.76e-01 
 -0.374   -8.33e-01   -1.125-3.38e-01log(time)
-1.53e-01  -3.52e-01   3.79e-02-0.597 -7.51e-01 
 -0.435-0.657 3.76e-01  -0.374   
-8.33e-01   -1.125-3.38e-01CV residual  -2.78e-17  
-5.55e-17  -6.94e-18 0.000  1.11e-16   
0.000 0.000-1.11e-16   0.0001.11e-16
0.000-5.55e-17
Sum of squares = 0Mean square = 0n = 12 
Overall (Sum over all 12 folds)   ms 1.18e-32 Warning messages:1: In 
predict.lm(subs.lm, newdata = df[rows.out, ]) :  prediction from a 
rank-deficient fit may be misleading2: In predict.lm(subs.lm, newdata = 
df[rows.out, ]) :  prediction from a rank-deficient fit may be misleading3: In 
CVlm(df = nihills, form.lm = formula(log(time) ~ .), plotit = Observed,  :  
 As there is 1 explanatory variable, cross-validation predicted values for a 
fold are not a linear function of corresponding overall predicted values.  
Lines that are shown for the different folds are approximate
 CVlm(df=nihills, 
 form.lm=formula(log(time)~dist+climb+timef),plotit=Observed,m=2)Analysis of 
 Variance Table
Response: log(time)  Df Sum Sq Mean Sq F value  Pr(F)dist   1  
 6.346.34  384.31 4.6e-14 ***climb  1   0.120.127.24  0.0145 *  
timef  1   0.190.19   11.29  0.0033 ** Residuals 19   0.310.02  
  ---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 


fold 1 Observations in test set: 11 Slieve Gullion McVeigh Classic 
Tollymore Mountain Moughanmore Hen  Cock Annalong Horseshoe   Rocky Meelbeg 
Meelmore Slieve Donard Seven Sevens Slieve GallionPredicted  -0.6801
-0.53759 -0.610 -0.7857   -0.84929  0.787 
-0.6824   -0.620-0.298 1.41 -0.587cvpred
 -0.7068-0.60517 -0.680 -0.7501   -0.80507  
1.053 -0.6856   -0.642-0.185 2.81 
-0.588log(time)  -0.7621-0.61413 -0.727 -0.7687 
  -0.79913  0.667 -0.6481   -0.789-0.053 
1.36 -0.452CV residual-0.0554-0.00896 
-0.047 -0.01860.00595 -0.386  0.0376   -0.147   
  0.132-1.45  0.136
Sum of squares = 2.31Mean square = 0.21n = 11 
fold 2 Observations in test set: 12 Binevenagh Glenariff Mountain 
Donard  Commedagh Slieve Martin Monument Race Loughshannagh Horseshoe Donard 
Forest Flagstaff to Carling Slieve Bearnagh Lurig Challenge Scrabo Hill Race 
BARF 

[R] sprintf in system command

2013-02-15 Thread Esam Tolba
hi all
I am using r (2.15.2) in windows 7 32bit
I want to execute an external program from r console. the program is a
command line program which needs the following format to start
  C:/Users/.../dssp-2.0.4-win32.exe -i
data_1.txt -o data_1.dssp
I used the system command as:
  system
('C:/Users/.../dssp-2.0.4-win32.exe -i data.txt -o data.dssp')
it worked.
Now I want to use the program on a list of files, so for that I used a
for loop and sprintf
   for (i in 1:10) {
   system
('C:/Users/.../dssp-2.0.4-win32.exe -i sprintf(data_%s.txt,i) -o
sprintf(data_%s.dssp,i)')
but I received the following error
No such file
Warning message:
running command 'C:/Users/.../dssp-2.0.4-win32.exe
-i sprintf(data_%s.txt,i) -o sprintf(data_%s.dssp,i)' had status 1


SO, what is my mistake? and is there a way to update the input file
name based on  the for loop counter

Best Regards,

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


Re: [R] sprintf in system command

2013-02-15 Thread jim holtman
I think what you want is something like this

 system(sprintf('C:/Users/.../dssp-2.0.4-win32.exe -i data%d.txt -o
data%d.dssp', i, i))

On Fri, Feb 15, 2013 at 12:33 PM, Esam Tolba eato...@gmail.com wrote:
 hi all
 I am using r (2.15.2) in windows 7 32bit
 I want to execute an external program from r console. the program is a
 command line program which needs the following format to start
   C:/Users/.../dssp-2.0.4-win32.exe -i
 data_1.txt -o data_1.dssp
 I used the system command as:
   system
 ('C:/Users/.../dssp-2.0.4-win32.exe -i data.txt -o data.dssp')
 it worked.
 Now I want to use the program on a list of files, so for that I used a
 for loop and sprintf
for (i in 1:10) {
system
 ('C:/Users/.../dssp-2.0.4-win32.exe -i sprintf(data_%s.txt,i) -o
 sprintf(data_%s.dssp,i)')
 but I received the following error
 No such file
 Warning message:
 running command 'C:/Users/.../dssp-2.0.4-win32.exe
 -i sprintf(data_%s.txt,i) -o sprintf(data_%s.dssp,i)' had status 1


 SO, what is my mistake? and is there a way to update the input file
 name based on  the for loop counter

 Best Regards,

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



-- 
Jim Holtman
Data Munger Guru

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

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


Re: [R] sprintf in system command

2013-02-15 Thread Rui Barradas

Hello,

In what follows I've used print(), not system().
The trick is to use paste0() to form the command.


for (i in 1:10) {
	cmd - paste0('C:/Users/.../dssp-2.0.4-win32.exe -i ', 
sprintf(data_%s.txt,i), ' -o ', sprintf(data_%s.dssp,i))

print(cmd)
}


Hope this helps,

Rui Barradas

Em 15-02-2013 17:33, Esam Tolba escreveu:

hi all
I am using r (2.15.2) in windows 7 32bit
I want to execute an external program from r console. the program is a
command line program which needs the following format to start
   C:/Users/.../dssp-2.0.4-win32.exe -i
data_1.txt -o data_1.dssp
I used the system command as:
   system
('C:/Users/.../dssp-2.0.4-win32.exe -i data.txt -o data.dssp')
it worked.
Now I want to use the program on a list of files, so for that I used a
for loop and sprintf
for (i in 1:10) {
system
('C:/Users/.../dssp-2.0.4-win32.exe -i sprintf(data_%s.txt,i) -o
sprintf(data_%s.dssp,i)')
but I received the following error
 No such file
 Warning message:
 running command 'C:/Users/.../dssp-2.0.4-win32.exe
-i sprintf(data_%s.txt,i) -o sprintf(data_%s.dssp,i)' had status 1


SO, what is my mistake? and is there a way to update the input file
name based on  the for loop counter

Best Regards,

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



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


Re: [R] 2 setGeneric's, same name, different method signatures

2013-02-15 Thread William Dunlap
Earlier you wrote
 changing me to b1 in A's bB(), allows A to use B's setGeneric()
and here mention intra-class name clashes (I'm not sure what you
mean by this).

In R, classes do not own (or contain) generic functions.  They can supply
methods for generic functions but a generic function is an object that
exists outside of any class.  C++ and Java do things differently.

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
 Behalf
 Of Greg Minshall
 Sent: Friday, February 15, 2013 3:35 AM
 To: Martin Morgan
 Cc: r-help@r-project.org
 Subject: Re: [R] 2 setGeneric's, same name, different method signatures
 
 Martin,
 
 fantastic.  thank you *very* much!  that clears lots of things up for
 me.
 
 (for the record: i think that setGeneric overwriting a previous is more
 surprising -- thus violating the principle of least surprise -- than one
 function overwriting a previous, in that we think of (or, one way to
 think of) OOP is that, after finagling to have no clashes on *class*
 names, we should not worry about other than intra-class name clashes.
 as i said, for the record.)
 
 thanks again.
 
 cheers, Greg
 
 Martin Morgan mtmor...@fhcrc.org wrote:
 
  Hi Greg --
 
  this setGeneric is over-writing the first, as would
 
  f = function() first
  f = function() second
  f() # second
 
  If you'd like to dispatch on a single argument, then
 
  setGeneric(one, function(x, ...) standardGeneric(one))
  setMethod(one, A, function(x, ...) A-method)
  setMetohd(one, B, function(x, y, ...) B-method)
 
  The '...' in the generic allow you to add arguments that are 'picked
  off' by methods. The user could provide any value for y, not only an
  object of class B.
 
  If you'd like to dispatch sometimes on two arguments then
 
  setGeneric(two, function(x, y, ...) standardGeneric(two))
  setMethod(two, c(A, ANY), function(x, y, ...) A,ANY-method)
  setMethod(two, c(B, B), function(x, y, ...) B,B-method)
 
  then two(new(A)), two(new(A), new(A)) and two(new(A),
  new(B)) end up in A,ANY,two-method while two(new(B), new(B))
  ends up in B,B,two-method. Other combinations are errors. One might
  instead not define A,ANY but instead
 
  setMethod(two, c(A, missing), function(x, y, ...) A,missing-method)
 
  and then two(new(A), new(A)) would be an error.
 
  Multiple dispatch is complicated, and perhaps best to avoid if possible.
 
  It's possible to write a generic and methods that dispatch on '...',
  with the requirement that all classes are the same; this in the spirit
  of comparing two B's and returning the smaller; see ?dotsMethods
  though again this is not a trivial use case.
 
  Hope that helps enough.
 
  Martin
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] reading data

2013-02-15 Thread arun
Hi,
#working directory data1 #changed name data to data1.  Added some files in each 
of sub directories a1, a2, etc.
 indx1- indx[indx!=]
lapply(indx1,function(x) list.files(x))
#[[1]]
#[1] a1.txt    m11kk.txt

#[[2]]
#[1] a2.txt    m11kk.txt

#[[3]]
#[1] a3.txt    m11kk.txt

#[[4]]
#[1] b1.txt    m11kk.txt

#[[5]]
#[1] b2.txt    b3.txt    m11kk.txt

[[6]]
[1] c1.txt    c2.txt    c3.txt    c4.txt   
[5] m11kk.txt

res-do.call(c,lapply(list.files(recursive=T)[grep(m11kk,list.files(recursive=T))],function(x)
 {names(x)-gsub(^(.*)\\/.*,\\1,x); lapply(x,function(y) 
read.table(y,header=TRUE,stringsAsFactors=FALSE,fill=TRUE))}))  #it seems like 
one of the rows of your file doesn't have 6 elements, so added fill=TRUE
head(res,2)
#$a1
 #    Id  M mm    x b  u  k  j    y    p    v
#1   aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#2 a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#3    aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#4   aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#5  aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#6    AA na  2 1972 0.0007000 11  3 AR   25  509  734

#$a2
 #    Id  M mm    x b  u  k  j    y    p    v
#1   aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#2 a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#3    aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#4   aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#5  aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#6    AA na  2 1972 0.0007000 11  3 AR   25  509  734


If you want the names to be group_a, group_b etc.
 names(res)-paste(group_,gsub(\\d+,,names(res)),sep=)
res[grep(group_b,names(res))]
$group_b
# Id  M mm    x b  u  k  j    y    p    v
#1   aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#2 a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#3    aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#4   aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#5  aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#6    AA na  2 1972 0.0007000 11  3 AR   25  509  734

#$group_b
 #    Id  M mm    x b  u  k  j    y    p    v
#1   aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#2 a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#3    aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#4   aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#5  aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#6    AA na  2 1972 0.0007000 11  3 AR   25  509  734

A.K.







- Original Message -
From: veracosta...@gmail.com veracosta...@gmail.com
To: smartpink...@yahoo.com
Cc: 
Sent: Friday, February 15, 2013 9:15 AM
Subject: reading data

Hi,
I post yesterday and you helped me. I have little problem.

At first, I never worked with regular expressions...

The code that you gave me it's ok, but my files are inside the folders 
a1,a2,a3. I try to explain better.

I have one folder named data. Inside this folder I have some other folders 
named a1,a2,b1,b2,...and inside of each one of that I have some files. I 
want only the file mm.txt (in all folders I have One file with this name).
The name of the folder give me the name of the group,but I need to read the 
file inside. And after, have group_a, group_b...because I need to work with 
this data grouped (and know the name of the group).

Thank you.


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


[R] Ho w Do I Get Cox Model Convergence After Multiple Imputation

2013-02-15 Thread Retsilisitsoe Moholisa
Due to missing data with some of my predictor variables I first do multiple
imputation as follows:
library(foreign)
library(Amelia)
library(norm)
set.seed(666)
M=10
impdat -
 
NVP[,c(X_t0,X_t,nvp,adstatus,t0rwfa,ageatran,whostage,t0rhfa,vlsupp,t0rwfh,t0rvl,t0rcd4pc,postrantb,resistance,id)]
round(apply(apply(impdat,c(1,2),is.na),2,mean),digits=3)   # percentage of
missing values
myimp - amelia(impdat, m=M, p2s=0,
noms=c(postrantb,vlsupp,whostage,resistance),cs=c(id),ts=c(X_t0),bounds=matrix(c(3,0,70,
4,0,100, 11,0,400, 12,0,50
),ncol=3,nrow=4,byrow=T),logs=nvp,polytime=3,splinetime=3,empri=1000,incheck=TRUE,tolerance=0.001)
The I do cox model as follows
# Categorization of relevant variables
for(i in 1:10){
myimp$imputations[[i]] -
cbind(myimp$imputations[[i]],cut(myimp$imputations[[i]]$t0rvl,breaks=c(-1,50,400,10)))
colnames(myimp$imputations[[i]])[15]-c(vlcat)
myimp$imputations[[i]] -
cbind(myimp$imputations[[i]],cut(myimp$imputations[[i]]$t0rcd4pc,breaks=c(-1,25,10)))
colnames(myimp$imputations[[i]])[16]-c(cd4pccat)
myimp$imputations[[i]] -
cbind(myimp$imputations[[i]],cut(myimp$imputations[[i]]$ageatran,breaks=c(-1,36,10)))
colnames(myimp$imputations[[i]])[17]-c(agecat)
myimp$imputations[[i]] -
cbind(myimp$imputations[[i]],cut(myimp$imputations[[i]]$t0rwfa,breaks=c(-1000,-3,-2,10)))
colnames(myimp$imputations[[i]])[18]-c(wfacat)
myimp$imputations[[i]]
-cbind(myimp$imputations[[i]],cut(myimp$imputations[[i]]$t0rwfh,breaks=c(-1000,-3,-2,10))
colnames(myimp$imputations[[i]])[19]-c(wfhcat)
myimp$imputations[[i]] -
cbind(myimp$imputations[[i]],cut(myimp$imputations[[i]]$t0rhfa,breaks=c(-1000,-3,-2,10)))
colnames(myimp$imputations[[i]])[20]-c(hfacat)
myimp$imputations[[i]] -
cbind(myimp$imputations[[i]],cut(myimp$imputations[[i]]$nvp,breaks=c(-1,1,10)))
colnames(myimp$imputations[[i]])[21]-c(nvpcat1)
myimp$imputations[[i]] -
cbind(myimp$imputations[[i]],cut(myimp$imputations[[i]]$nvp,breaks=c(-1,4,10)))
colnames(myimp$imputations[[i]])[22]-c(nvpcat10)}
head(myimp$imputations[[1]])
##
Then I perform cox regression as follows
m2_1-coxph(Surv(X_t0,X_t, vlsupp) ~ nvp  + as.factor(cd4pccat) +
as.factor(vlcat) + as.factor(agecat) + as.factor(whostage)  +
as.factor(hfacat) + as.factor(wfacat) + as.factor(wfhcat) +
as.factor(resistance) + as.factor(postrantb) +
cluster(id),data=myimp$imputations[[1]],method=breslow,robust=TRUE)
summary(m2_1)
The I get the following eWarning message:
In fitter(X, Y, strats, offset, init, control, weights = weights,  :
  Ran out of iterations and did not converge

I have tried to increase the maximum number of iterations,however I still
get the same warning message.
I have also looked at change my initial estimates but I get the following
error
Error in fitter(X, Y, strats, offset, init, control, weights = weights,  :
  Wrong length for inital values
Does anybody know have any suggestions as to I can get the model to
converge
Kind Regards
-- 
Retsilisitsoe Moholisa, MSc Medical Biochemistry
PhD Student in Pharmacometrics
Pharmacometrics Lab
Department of Clinical Pharmacology
Faculty of Health Sciences
University of Cape Town
Tel:021-404 7719
Cell:0712801254

[[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] mixed effects regression with non-independent data

2013-02-15 Thread Bone, Jonathan
Hi,
I will be performing mixed effects regression on subjects total scores from 2 
player games (prisoners dilemma) that played. I am aware that including both 
players score from a game will cause problems due to non-independence. Is there 
a way that to deal with this apart from randomly picking one subject from each 
game for the analysis (and so losing half the data). Is there a way to 
introduce this into the model instead, perhaps as a random effect??

Thanks,

Jonathan
[X]
[X]
[X]
[X]
[X]
[X]
[X]
[X]

[[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] Mixed models with missing data

2013-02-15 Thread Bone, Jonathan
Hi,

I am creating a mixed model based on a experiment where each subject has 2 
repeats. In some instances though there is only data for one of a given 
subjects repeats for most there is data for both. Can I still justify having 
subject as a random effect?

Thanks,

Jonathan

[X]
[X]
[X]
[X]
[X]
[X]
[X]
[X]

[[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] reading data

2013-02-15 Thread arun
HI,

Just to add:

res-do.call(c,lapply(list.files(recursive=T)[grep(m11kk,list.files(recursive=T))],function(x)
 {names(x)-gsub(^(.*)\\/.*,\\1,x); lapply(x,function(y) 
read.table(y,header=TRUE,stringsAsFactors=FALSE,fill=TRUE))}))  #it seems like 
one of the rows of your file doesn't have 6 elements, so added fill=TRUE
 names(res)-paste(group_,gsub(\\d+,,names(res)),sep=)
res[grep(group_b,names(res))]

I am not sure how you want the grouped data to look like.  If you want 
something like this:
res1-do.call(rbind,res)
res2-lapply(split(res1,gsub([.0-9],,row.names(res1))),function(x) 
{row.names(x)-1:nrow(x);x})
res2
#$group_a
 #     Id  M mm    x b  u  k  j    y    p    v
#1    aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#2  a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#3 aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#4    aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#5   aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#6 AA na  2 1972 0.0007000 11  3 AR   25  509  734
#7    aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#8  a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#9 aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#10   aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#11  aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#12    AA na  2 1972 0.0007000 11  3 AR   25  509  734
#13   aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#14 a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#15    aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#16   aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#17  aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#18    AA na  2 1972 0.0007000 11  3 AR   25  509  734

#$group_b
 #     Id  M mm    x b  u  k  j    y    p    v
#1    aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#2  a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#3 aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#4    aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#5   aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#6 AA na  2 1972 0.0007000 11  3 AR   25  509  734
#7    aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#8  a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#9 aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#10   aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#11  aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#12    AA na  2 1972 0.0007000 11  3 AR   25  509  734

#$group_c
 #    Id  M mm    x b  u  k  j    y    p    v
#1   aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#2 a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#3    aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#4   aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#5  aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#6    AA na  2 1972 0.0007000 11  3 AR   25  509  734


#or if you want it like this:
res2-split(res,names(res))

res2[[group_b]]
#$group_b
# Id  M mm    x b  u  k  j    y    p    v
#1   aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#2 a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#3    aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#4   aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#5  aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#6    AA na  2 1972 0.0007000 11  3 AR   25  509  734

#$group_b
 #    Id  M mm    x b  u  k  j    y    p    v
#1   aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#2 a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#3    aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#4   aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#5  aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#6    AA na  2 1972 0.0007000 11  3 AR   25  509  734

Hope this helps.
A.K.



- Original Message -
From: veracosta...@gmail.com veracosta...@gmail.com
To: smartpink...@yahoo.com
Cc: 
Sent: Friday, February 15, 2013 9:15 AM
Subject: reading data

Hi,
I post yesterday and you helped me. I have little problem.

At first, I never worked with regular expressions...

The code that you gave me it's ok, but my files are inside the folders 
a1,a2,a3. I try to explain better.

I have one folder named data. Inside this folder I have some other folders 
named a1,a2,b1,b2,...and inside of each one of that I have some files. I 
want only the file mm.txt (in all folders I have One file with this name).
The name of the folder give me the name of the group,but I need to read the 
file inside. And after, have group_a, group_b...because I need to work with 
this data grouped (and know the name of the group).

Thank you.

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

Re: [R] reading data

2013-02-15 Thread arun
HI,
No problem.
?c() for concatenate to vector or list().
If I use do.call(cbind,..) or do.call(rbind,...)

do.call(cbind,lapply(list.files(recursive=T)[grep(m11kk,list.files(recursive=T))],function(x)
 {names(x)-gsub(^(.*)\\/.*,\\1,x); lapply(x,function(y) 
read.table(y,header=TRUE,stringsAsFactors=FALSE,fill=TRUE))}))  
#   [,1]    [,2]    [,3]    [,4]    [,5]    [,6]   
#a1 List,11 List,11 List,11 List,11 List,11 List,11


 
do.call(rbind,lapply(list.files(recursive=T)[grep(m11kk,list.files(recursive=T))],function(x)
 {names(x)-gsub(^(.*)\\/.*,\\1,x); lapply(x,function(y) 
read.table(y,header=TRUE,stringsAsFactors=FALSE,fill=TRUE))}))  
# a1 
#[1,] List,11
#[2,] List,11
#[3,] List,11
#[4,] List,11
#[5,] List,11
#[6,] List,11
ie.
list within in a list

 
restrial-lapply(list.files(recursive=T)[grep(m11kk,list.files(recursive=T))],function(x)
 {names(x)-gsub(^(.*)\\/.*,\\1,x); lapply(x,function(y) 
read.table(y,header=TRUE,stringsAsFactors=FALSE,fill=TRUE))})
 str(restrial)
#List of 6
# $ :List of 1
  #..$ a1:'data.frame':    6 obs. of  11 variables:
  .#. ..$ Id: chr [1:6] aAA a aA aAA ...
  #.. ..$ M : chr [1:6] 1 1 2 1 ...
  #. ..$ mm: int [1:6] 2 2 1 2 3 2
  #. ..$ x : int [1:6] 739 2263 1 1965 3660 1972
  -
str(res)
#List of 6
# $ a1:'data.frame':    6 obs. of  11 variables:
 # ..$ Id: chr [1:6] aAA a aA aAA ...
  #..$ M : chr [1:6] 1 1 2 1 ...
 # ..$ mm: int [1:6] 2 2 1 2 3 2
 # ..$ x : int [1:6] 739 2263 1 1965 3660 1972
-

You mentioned about naming this to group_a,group_b. etc..
 names(res)-paste(group_,gsub(\\d+,,names(res)),sep=)
res2-split(res,names(res))

res3- lapply(res2,function(x) 
{names(x)-paste(gsub(.*_,,names(x)),1:length(x),sep=);x})
 res3$group_a
$a1
# Id  M mm    x b  u  k  j    y    p    v
#1   aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#2 a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#3    aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#4   aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#5  aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#6    AA na  2 1972 0.0007000 11  3 AR   25  509  734

#$a2
# Id  M mm    x b  u  k  j    y    p    v
#1   aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#2 a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#3    aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#4   aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#5  aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#6    AA na  2 1972 0.0007000 11  3 AR   25  509  734

#$a3
 #    Id  M mm    x b  u  k  j    y    p    v
#1   aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#2 a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#3    aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#4   aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#5  aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#6    AA na  2 1972 0.0007000 11  3 AR   25  509  734
A.K.

From: Vera Costa veracosta...@gmail.com
To: arun smartpink...@yahoo.com 
Sent: Friday, February 15, 2013 12:39 PM
Subject: Re: reading data


Thank you very much and sorry my questions.

But this code isn't grouping for letters sure? I mean, a1,a2,a3 is the same 
group, (the first letter give me the name of the group)

Another question, in do.call, you did do.call (c,.) .What is c?

Sorry



2013/2/15 arun smartpink...@yahoo.com

HI,

Just to add:


res-do.call(c,lapply(list.files(recursive=T)[grep(m11kk,list.files(recursive=T))],function(x)
 {names(x)-gsub(^(.*)\\/.*,\\1,x); lapply(x,function(y) 
read.table(y,header=TRUE,stringsAsFactors=FALSE,fill=TRUE))}))  #it seems like 
one of the rows of your file doesn't have 6 elements, so added fill=TRUE

 names(res)-paste(group_,gsub(\\d+,,names(res)),sep=)
res[grep(group_b,names(res))]

I am not sure how you want the grouped data to look like.  If you want 
something like this:
res1-do.call(rbind,res)
res2-lapply(split(res1,gsub([.0-9],,row.names(res1))),function(x) 
{row.names(x)-1:nrow(x);x})
res2
#$group_a

 #     Id  M mm    x b  u  k  j    y    p    v
#1    aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#2  a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#3 aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#4    aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#5   aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#6 AA na  2 1972 0.0007000 11  3 AR   25  509  734
#7    aAA  1  2  739 0.1257000  2  2 AA    2 8867 8926
#8  a  1  2 2263 0.0004000  2  2 AR    4 7640 8926
#9 aA  2  1    1 0.0845435  2 AA  2 6790 734,1092   NA
#10   aAA  1  2 1965 0.0007000  4  3 AR    2    11616 8926
#11  aAAA  1  3 3660 0.0008600 18  3 AA    2    20392  496
#12    AA na  2 1972 0.0007000 11  3 AR   25  509  734
#13   aAA  1  2  739 

[R] unbalanced design

2013-02-15 Thread Bond, Stephen
Please, help with a formula for dealing with unbalanced design:

To see the counts:
aggregate(dfa$CertId,by=list(type=dfa$ComType,stat=dfa$StatusCodeId),length)

  type stat x
1C1  6571
2C3 28957
3C8 12390
4C   11 12415
5E   13 9
6R   1351
7E   15  2079
8R   15  6692

I would like to have a slope for statuses 1,3,8,11,13 and two slopes for status 
15 one for type E and one for type R.
I tried nesting, but it assumes that all levels exist for each factor and 
complains about singular model matrix. Is there a theoretically proper way to 
deal with this or I should just relabel status 15 and make it 16 for type R and 
regress on status alone??
Thanks everybody


Stephen B

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


[R] convert list into a time series

2013-02-15 Thread twynne
I am trying to use the SeasonalMannKendall function in the Kendall package.
My dataset (alb_data) is in the same format as the example dataset 
(manaus) in the package.

  class(manaus)
[1] ts
  is.ts(manaus)
[1] TRUE
  typeof(manaus)
[1] double
  alb_data=read.table(R:/albemarle_manken.txt, header=T)
  head(alb_data)
   yearJanFebMarAprMayJunJulAug Sep
OctNovDec
1 1997NaNNaNNaNNaNNaNNaNNaNNaN 23.406 
16.166 16.057 16.803
2 1998 14.157 14.605 14.112 17.709 13.800 14.338 14.157 17.404 17.725 
15.429 16.090 18.061
3 1999 14.888 13.837 15.929 13.637 16.020 14.699 15.987 15.212 14.752 
15.935 13.397 21.725
4 2000 16.562 18.125 19.600 17.971 16.454 15.129 13.901 21.664 17.675 
13.793 13.464 16.452
5 2001 15.706 16.417 13.324 14.117 13.550 15.825 14.687 14.844 15.006 
14.793 13.489 15.726
6 2002 11.777 11.775 11.564 13.141 14.462 15.517 18.801 15.652 17.127 
14.070 14.626 17.964
  class(alb_data)
[1] data.frame
  is.ts(alb_data)
[1] FALSE
  typeof(alb_data)
[1] list


Can anyone please help me get my dataset into the correct format so that 
I can run the test?
Thank you for your time and help!




--
View this message in context: 
http://r.789695.n4.nabble.com/convert-list-into-a-time-series-tp4658716.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.


Re: [R] data formatting

2013-02-15 Thread arun


Dear Eliza,

Try this:

Lines1-readLines(textConnection(1911.01.01   7.87
1911.01.02   9.26 
1911.01.03   8.06 
1911.01.04   8.13 
1911.01.05  12.90 
1911.02.06   5.45 
1911.02.07   3.26 
1911.03.08   5.70 
1911.03.09   9.24 
1911.04.10   7.60 
1911.05.11  14.82 
1911.05.12  14.10 
1911.06.13   7.87 
1911.06.14   9.26 
1911.07.15   8.06 
1911.07.16   8.13 
1911.08.17  12.90 
1911.08.18   5.45 
1911.09.19   3.26 
1911.09.20   5.70 
1911.10.21   9.24 
1911.10.22   7.60 
1911.11.23  14.82 
1911.12.24  14.10)) 

Lines2-Lines1[Lines1!=]
library(stringr)
 str_count(Lines2,  )
# [1] 7 7 7 7 6 7 7 7 7 7 6 6 7 7 7 7 6 7 7 7 7 7 6 6


Lines2[str_count(Lines2, )==7]- str_replace(Lines2[str_count(Lines2, 
)==7],\\s+, ) #reduced 2 spaces

 Lines2[str_count(Lines2, )==6]- str_replace(Lines2[str_count(Lines2, 
)==6],\\s+,    ) #reduced 2 spaces
 str_count(Lines2, )
# [1] 5 5 5 5 4 5 5 5 5 5 4 4 5 5 5 5 4 5 5 5 5 5 4 4
substr(Lines2[substr(Lines2,6,6)==0|substr(Lines2,9,9)==0],6,6)- 
substr(Lines2[substr(Lines2,6,6)==0|substr(Lines2,9,9)==0],9,9)- 
str_count(Lines2, ) #see the difference in space.  This counts all the space. 
 Here 2 white space are added to replace 0
# [1] 7 7 7 7 6 7 7 7 7 6 5 5 6 6 6 6 5 6 6 6 5 5 4 4
Lines2
# [1] 1911. 1. 1 7.87 1911. 1. 2 9.26 1911. 1. 3 8.06
# [4] 1911. 1. 4 8.13 1911. 1. 5    12.90 1911. 2. 6 5.45
# [7] 1911. 2. 7 3.26 1911. 3. 8 5.70 1911. 3. 9 9.24
#[10] 1911. 4.10 7.60 1911. 5.11    14.82 1911. 5.12    14.10
#[13] 1911. 6.13 7.87 1911. 6.14 9.26 1911. 7.15 8.06
#[16] 1911. 7.16 8.13 1911. 8.17    12.90 1911. 8.18 5.45
#[19] 1911. 9.19 3.26 1911. 9.20 5.70 1911.10.21 9.24
#[22] 1911.10.22 7.60 1911.11.23    14.82 1911.12.24    14.10

A.K.

From: eliza botto eliza_bo...@hotmail.com
To: smartpink...@yahoo.com smartpink...@yahoo.com 
Sent: Friday, February 15, 2013 12:38 PM
Subject: data formatting



Dear Arun,
[text file is also attached if format is changed]
i need to data managing genius expertise on the following issue.
i have data like the following table

1911.01.01   7.87 ##(7 spaces between the columns)
1911.01.02   9.26 ##(7 spaces between the columns)
1911.01.03   8.06 ##(7 spaces between the columns)
1911.01.04   8.13 ##(7 spaces between the columns)
1911.01.05  12.90 ##(6 spaces between the columns)
1911.02.06   5.45 ##(7 spaces between the columns)
1911.02.07   3.26 ##(7 spaces between the columns)
1911.03.08   5.70 ##(7 spaces between the columns)
1911.03.09   9.24 ##(7 spaces between the columns)
1911.04.10   7.60 ##(7 spaces between the columns)
1911.05.11  14.82 ##(6 spaces between the columns)
1911.05.12  14.10 ##(6 spaces between the columns)
1911.06.13   7.87 ##(7 spaces between the columns)
1911.06.14   9.26 ##(7 spaces between the columns) 
1911.07.15   8.06 ##(7 spaces between the columns) 
1911.07.16   8.13 ##(7 spaces between the columns) 
1911.08.17  12.90 ##(6 spaces between the columns) 
1911.08.18   5.45 ##(7 spaces between the columns) 
1911.09.19   3.26 ##(7 spaces between the columns) 
1911.09.20   5.70 ##(7 spaces between the columns)
1911.10.21   9.24 ##(7 spaces between the columns)
1911.10.22   7.60 ##(7 spaces between the columns)
1911.11.23  14.82 ##(6 spaces between the columns)
1911.12.24  14.10 ##(6 spaces between the columns)
and i want it to be in the following manner and afterwards i want to save that 
file in .txt format.
 1911. 1. 1 7.87 ##(5 spaces between the columns)
 1911. 1. 2 9.26 ##(5 spaces between the columns)
 1911. 1. 3 8.06 ##(5 spaces between the columns)
 1911. 1. 4 8.13 ##(5 spaces between the columns)
 1911. 1. 5    12.90 ##(4 spaces between the columns)
 1911. 2. 6 5.45 ##(5 spaces between the columns)
 1911. 2. 7 3.26 ##(5 spaces between the columns)
 1911. 3. 8 5.70 ##(5 spaces between the columns)
 1911. 3. 9 9.24 ##(5 spaces between the columns)
 1911. 4.10 7.60 ##(5 spaces between the columns)
 1911. 5.11    14.82 ##(4 spaces between the columns)
 1911. 5.12    14.10 ##(4 spaces between the columns)
 1911. 6.13 7.87 ##(5 spaces between the columns)
 1911. 6.14 9.26 ##(5 spaces between the columns)
 1911. 7.15 8.06 ##(5 spaces between the columns)
 1911. 7.16 8.13 ##(5 spaces between the columns)
 1911. 8.17    12.90 ##(4 spaces between the columns)
 1911. 8.18 5.45 ##(5 spaces between the columns)
 1911. 9.19 3.26 ##(5 spaces between the columns)
 1911. 9.20 5.70 ##(5 spaces between the columns)
 1911.10.21 9.24 ##(5 spaces between the columns)
 1911.10.22 7.60 ##(5 spaces between the columns)
 1911.11.23    14.82 ##(4 spaces between the columns)
 1911.12.24    14.10 ##(4 spaces between the columns)
you could see that spaces between the columns 

[R] minimizing a numerical integration

2013-02-15 Thread Aya Anas
Dear all,


I am a new user to R and I am using pracma and nloptr libraries to minimize
a numerical integration subject to a single constraint . The integrand
itself is somehow a complicated function of x and y that is computed
through several steps.  i formulated the integrand in a separate function
called f which is a function of x y. I want to find the optimal value of x
such that the integration over y is minimum. Here is my code, it is not
working. I got an error: 'y' is missing


library('nloptr')
 library('pracma')



f - function(x,y) {#here i should put the commands representing my function

return(   )
}


 #constraint function


eval_g0 - function(x) {

return( )
}


# objective function


eval_f0 - function(x) {
romberg(f, 0.5, 0.5001)}


ARL1 - nloptr( x0=c(0.653),
eval_f=eval_f0,
lb = c(0),
ub = c(6),
eval_g_ineq = eval_g0,
opts = list(algorithm=NLOPT_LN_COBYLA, maxeval=1000),
 )
print( )



 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] 3D-plots of 2D-grids

2013-02-15 Thread David Winsemius

On Feb 15, 2013, at 1:14 AM, jas wrote:

 Hello David,
 thanks again for your reply.
 
 Two things remain unclear. That the data is disjointed is ok, as there are
 only values in hectares, where there are actually buildings and stuff,
 forest/nature is NA. 
 
 The presp-scan that I get from persp(grd) has no similarity to the image in
 2d. My guess is that the order and regularity in the data somehow gets lost
 in the process of making the matrix?! 

I think the disjointedness is the problem. There are many values where there is 
no adjacent value in one direction of another and so no tessellation can be 
formed. If you look at the image result and note the places where there are 
solid values in both x and y directions I think the overall patterns match 
up. The image result is more faithful to the data.

 
 The grid-extraction that I uploaded contains 40x20 cells. if its a bigger
 grid with 7million cells, does it still work the same way or can the vectors
 only be a maximum of 100 cells? 

Might be a performance problem although I don't think it is theoretically 
impossible. An 8000 x 8000 matrix consumed my full CPU resources and 
essentially locked up my session. I'm in the process of deciding when to halt 
it.

Are you aware that there are far more knowledgeable persons than I that hang 
out at the R-SIG-Geo list?

-- 
David

 
 Thank you for your help
 jas
 
 
 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/3D-plots-of-2D-grids-tp4658517p4658639.html
 Sent from the R help mailing list archive at Nabble.com.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

David Winsemius
Alameda, CA, USA

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


Re: [R] 3D-plots of 2D-grids

2013-02-15 Thread David Winsemius
Compare:

persp(matrix(c(1,1,NA,NA,
  1,1,NA,NA,
  1,NA,NA,NA,
   1,NA,NA,1), 4,4), zlim=c(0,2))

persp(matrix(c(1,1,0,0,
 1,1,0,0,
 1,0,0,0,
  1,0,0,1), 4,4), zlim=c(0,2))

-- 
David.


On Feb 15, 2013, at 12:01 PM, David Winsemius wrote:

 
 On Feb 15, 2013, at 1:14 AM, jas wrote:
 
 Hello David,
 thanks again for your reply.
 
 Two things remain unclear. That the data is disjointed is ok, as there are
 only values in hectares, where there are actually buildings and stuff,
 forest/nature is NA. 
 
 The presp-scan that I get from persp(grd) has no similarity to the image in
 2d. My guess is that the order and regularity in the data somehow gets lost
 in the process of making the matrix?! 
 
 I think the disjointedness is the problem. There are many values where there 
 is no adjacent value in one direction of another and so no tessellation can 
 be formed. If you look at the image result and note the places where there 
 are solid values in both x and y directions I think the overall patterns 
 match up. The image result is more faithful to the data.
 
 
 The grid-extraction that I uploaded contains 40x20 cells. if its a bigger
 grid with 7million cells, does it still work the same way or can the vectors
 only be a maximum of 100 cells? 
 
 Might be a performance problem although I don't think it is theoretically 
 impossible. An 8000 x 8000 matrix consumed my full CPU resources and 
 essentially locked up my session. I'm in the process of deciding when to halt 
 it.
 
 Are you aware that there are far more knowledgeable persons than I that hang 
 out at the R-SIG-Geo list?
 
 -- 
 David
 
 
 Thank you for your help
 jas
 
 
 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/3D-plots-of-2D-grids-tp4658517p4658639.html
 Sent from the R help mailing list archive at Nabble.com.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 David Winsemius
 Alameda, CA, USA
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

David Winsemius
Alameda, CA, USA

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


Re: [R] Mixed models with missing data

2013-02-15 Thread Bert Gunter
This has basically nothing to do with R, so please don't post here.

You may wish to try the r-sig-mixed-models list, however. They are
more sympathetic to such questions -- and what are likely to be the
torrent from you that follows.

-- Bert

On Fri, Feb 15, 2013 at 9:17 AM, Bone, Jonathan
jonathan.bone...@ucl.ac.uk wrote:
 Hi,

 I am creating a mixed model based on a experiment where each subject has 2 
 repeats. In some instances though there is only data for one of a given 
 subjects repeats for most there is data for both. Can I still justify having 
 subject as a random effect?

 Thanks,

 Jonathan

 [X]
 [X]
 [X]
 [X]
 [X]
 [X]
 [X]
 [X]

 [[alternative HTML version deleted]]

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



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

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

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


Re: [R] 2 setGeneric's, same name, different method signatures

2013-02-15 Thread William Dunlap
I thought you were thinking of the R class system (the S3 and S4 ones
anyway) as if it were C++'s.  It is quite different.

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com


 -Original Message-
 From: Greg Minshall [mailto:minsh...@umich.edu]
 Sent: Friday, February 15, 2013 1:37 PM
 To: William Dunlap
 Cc: r-help@r-project.org
 Subject: Re: 2 setGeneric's, same name, different method signatures
 
 William,
 
  and here mention intra-class name clashes (I'm not sure what you
  mean by this).
 
 sorry, i meant, in something like C++, if i have a class Foo and you
 have a class Bar, then i can invent whatever method names/signatures i
 want, independent of whatever method names/signatures *you* want.  *i*
 just need to make sure i don't introduce incompatible methods with the
 same name/signature *within* Foo.
 
 
 cheers, Greg

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] lattice 3x3 plot: force common y-limits accross rows and align x-axes

2013-02-15 Thread Boris.Vasiliev
Good afternoon,

I would like to ask for help in controlling y-axis limits and labels in
lattice doplots.  Unfortunately, the problem is somewhat convoluted,
please bear with the long explanation.

I would like to create a 3x3 lattice of dotplots, say subject ~ count.
The plot is conditioned on variables treatment and risk: subject ~ count
| treatment + risk.  In the experiment,  not all subjects were exposed
to all combinations of treatment and risk.  For each risk, I would like
to show subject ~ count | treatment and order the subjects by the total
count.  At the same time, I would like the x-axes to be the same in all
panels and aligned by columns.

Here is a sample data set:

# raw data
df - data.frame(subject=c('A','A','A','BB','BB','CCC','CCC','CCC',
   'DD','DD','A','A','A','','',
   'A','A','B','B'),
 
risk=c('high','high','high','high','high','high','high','high',
'med','med','med','med','med','med','med',
'low','low','low','low'),
 
treatment=c('none','optX','optZ','none','optZ','none','optX','optZ',
 
'none','optZ','none','optX','optZ','none','optZ',
 'none','optX','none','optZ'),
 count=c(5,10,2,3,5,8,1,2,
 3,7,10,2,5,15,2,
 7,7,10,8))
# re-level factors
df$risk - factor(df$risk,levels=c('low','med','high'))
df$treatment - factor(df$treatment,levels=c('none','optX','optZ'))


##  df
##subject risk treatment count
## 1A high  none 5
## 2A high  optX10
## 3A high  optZ 2
## 4   BB high  none 3
## 5   BB high  optZ 5
## 6  CCC high  none 8
## 7  CCC high  optX 1
## 8  CCC high  optZ 2
## 9   DD  med  none 3
## 10  DD  med  optZ 7
## 11   A  med  none10
## 12   A  med  optX 2
## 13   A  med  optZ 5
## 14  med  none15
## 15  med  optZ 2
## 16   A  low  none 7
## 17   A  low  optX 7
## 18   B  low  none10
## 19   B  low  optZ 8

One way to plot the data is to break-up the data into sub-frames, one
frame for each risk, order subjects by total counts, create dotplots,
and merge with trellis.c().  This almost works but in the merged plot I
cannot decrease column spacing to be small enough.  Also, the output of
trellis.c() would not work with useOuterStrips() which I really like.
My code is in TRY ONE below.

Another way to create the plot is specify y-limits for each panel with
custom prepanel and panel functions.  For each panel, the data-frame for
the panel row is isolated, subjects in the data-frame for the current
row are ordered by counts, panel y-limits are set to the re-ordered
levels, y-data for each panel is releveled, and data plotted with
standard panel.dotplot().  This somewhat works but lattice does not
honour the user-defined y-limits and labels are not correct.  I suspect
that it is not correct to use  y-relation=same in this case but free
and sliced do not give correct results too. My code in in TRY TWO
below.

If anybody can offer any assistance with this problem, it would be much
appreciated,
Sincerely,
Boris.


 BEGIN TRY ONE - MERGE LATTICE PLOTS 
library(lattice)
library(latticeExtra)
library(grid)

for (irisk in levels(df$risk)) {
  # subset data frame
  df.irisk - subset(df,risk==irisk)

  # order subjects by total count; store levels of subjectx variables
  # for later re-use in panel labels
  df.irisk$subjectx - df.irisk$subject[,drop=TRUE]
  df.irisk$subjectx - reorder(df.irisk$subjectx,df.irisk$count,sum)
  assign(paste('sbjx.',irisk,sep=''),levels(df.irisk$subjectx))

  # create dotplot and store it in oltc.{irisk} variable
  oltc.irisk - dotplot(subjectx~count|treatment,data=df.irisk,
layout=c(3,1),type=c('p','h'),
xlim=c(-1,16),origin=0,
xlab=,ylab=)
  assign(paste('oltc.',irisk,sep=''),oltc.irisk)
}

# combine everthing in one plot
oltc - c(low=oltc.low,med=oltc.med,high=oltc.high)
print(oltc)

# get rid of variable labels in middle and right column; decrease
# distance between columns.  But can't make inter-column spaces
# small enought and get rid of the panels in all but top rows.
laywid - trellis.par.get('layout.widths')
laywid$between - -5
laywid$axis.panel - 0.7
yscales - list(labels=list(sbjx.low,NULL,NULL,
sbjx.med,NULL,NULL,
sbjx.high,NULL,NULL))
oltd - update(oltc,scales=list(y=yscales),
par.settings=list(layout.widths=laywid))
print(oltd)

 END TRY ONE - MERGE LATTICE PLOTS 

 BEGIN TRY TWO - CUSTOM PREPANEL AND PANEL FUNCTIONS 

prepanel.dotplot.x - function(x,y,type,subscripts,...,data=NULL) {
  # find data-frame that corresponds to the entire row of 

Re: [R] file copy to password protected network drive

2013-02-15 Thread Kumar Mainali
TACC = Texas Advanced Computing Center, OS is Linux. I have used shell
commands for simple copy procedure. This time, I have to copy tens of
thousand files, and I can write an R code to find all of those files. I am
sure shell command can be used to get that done but I have almost no
knowledge of it.

On Fri, Feb 15, 2013 at 2:10 AM, Jeff Newmiller jdnew...@dcn.davis.ca.uswrote:

 I have no idea what TACC is, or what your OS is, or what file networking
 scheme your system is using, and these issues are all outside the topic
 area for this list. You should go talk to your network admin or local help
 desk about how to accomplish this task at the command line, and then if it
 involves anything more than ordinary file system access then you should
 Google for how to invoke shell commands on your OS (shell/system/system2).
 ---
 Jeff NewmillerThe .   .  Go Live...
 DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live
 Go...
   Live:   OO#.. Dead: OO#..  Playing
 Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
 /Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
 ---
 Sent from my phone. Please excuse my brevity.

 Kumar Mainali kpmain...@gmail.com wrote:

 I am trying to copy files to a password protected drive which is
 Ranch at
 TACC from another network drive. I am logged in to the source drive and
 can
 run R there. The following code does not even find the destination
 folder.
 
 file.copy(sourcedrive/file.tar, 
 usern...@ranch.tacc.utexas.edu/uniqueID/destinationfolder/file.tar,
 overwrite = FALSE)
 
 Thanks in advance,
 Kumar
 
 On Thu, Feb 14, 2013 at 4:41 AM, e-letter inp...@gmail.com wrote:
 
  Readers,
 
  For this data set:
 
  testvalues-c(10,20,30,40)
 
  How to amend the plot instruction:
 
  plot(testvalues,ann=FALSE,type='l',yaxt='n',xaxt='n')
 
  so that x axis ticks labels can be added to existing graph with
  arbitrary value such as 0,100,200,300)?
 
  Thanks in advance.
 
  --
  r2151
 
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide
  http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
 




-- 
Section of Integrative Biology
University of Texas at Austin
Austin, Texas 78712, USA

[[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] 2 setGeneric's, same name, different method signatures

2013-02-15 Thread Hadley Wickham
This is unfortunately reinforced by the (Not So) Short Introduction
to S4 Object Oriented Programming in R - I wouldn't recommend that
document to learn about S4.

The most important thing to get about OO in R is that methods belong
to generic functions, not like classes, as in most other programming
languages.  If you don't get that, you are really going to struggle
with S3 and S4 and you will wonder what on earth the designers of R
were thinking when they created them.

Hadley

On Fri, Feb 15, 2013 at 3:42 PM, William Dunlap wdun...@tibco.com wrote:
 I thought you were thinking of the R class system (the S3 and S4 ones
 anyway) as if it were C++'s.  It is quite different.

 Bill Dunlap
 Spotfire, TIBCO Software
 wdunlap tibco.com


 -Original Message-
 From: Greg Minshall [mailto:minsh...@umich.edu]
 Sent: Friday, February 15, 2013 1:37 PM
 To: William Dunlap
 Cc: r-help@r-project.org
 Subject: Re: 2 setGeneric's, same name, different method signatures

 William,

  and here mention intra-class name clashes (I'm not sure what you
  mean by this).

 sorry, i meant, in something like C++, if i have a class Foo and you
 have a class Bar, then i can invent whatever method names/signatures i
 want, independent of whatever method names/signatures *you* want.  *i*
 just need to make sure i don't introduce incompatible methods with the
 same name/signature *within* Foo.


 cheers, Greg

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



-- 
Chief Scientist, RStudio
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.


Re: [R] lattice 3x3 plot: force common y-limits accross rows and align x-axes

2013-02-15 Thread Duncan Mackay

Hi Boris

Not sure what you mean exactly
try

library(latticeExtra)
useOuterStrips(dotplot(count ~ subject|risk*treatment,df))

if you want to change the order of the subjects in each panel and an 
index column and plot the index column instead of subject and change 
the scales to suit.


HTH

Duncan

Duncan Mackay
Department of Agronomy and Soil Science
University of New England
Armidale NSW 2351
Email: home: mac...@northnet.com.au



At 07:54 16/02/2013, you wrote:

Good afternoon,

I would like to ask for help in controlling y-axis limits and labels in
lattice doplots.  Unfortunately, the problem is somewhat convoluted,
please bear with the long explanation.

I would like to create a 3x3 lattice of dotplots, say subject ~ count.
The plot is conditioned on variables treatment and risk: subject ~ count
| treatment + risk.  In the experiment,  not all subjects were exposed
to all combinations of treatment and risk.  For each risk, I would like
to show subject ~ count | treatment and order the subjects by the total
count.  At the same time, I would like the x-axes to be the same in all
panels and aligned by columns.

Here is a sample data set:

# raw data
df - data.frame(subject=c('A','A','A','BB','BB','CCC','CCC','CCC',
   'DD','DD','A','A','A','','',
   'A','A','B','B'),

risk=c('high','high','high','high','high','high','high','high',
'med','med','med','med','med','med','med',
'low','low','low','low'),

treatment=c('none','optX','optZ','none','optZ','none','optX','optZ',

'none','optZ','none','optX','optZ','none','optZ',
 'none','optX','none','optZ'),
 count=c(5,10,2,3,5,8,1,2,
 3,7,10,2,5,15,2,
 7,7,10,8))
# re-level factors
df$risk - factor(df$risk,levels=c('low','med','high'))
df$treatment - factor(df$treatment,levels=c('none','optX','optZ'))


##  df
##subject risk treatment count
## 1A high  none 5
## 2A high  optX10
## 3A high  optZ 2
## 4   BB high  none 3
## 5   BB high  optZ 5
## 6  CCC high  none 8
## 7  CCC high  optX 1
## 8  CCC high  optZ 2
## 9   DD  med  none 3
## 10  DD  med  optZ 7
## 11   A  med  none10
## 12   A  med  optX 2
## 13   A  med  optZ 5
## 14  med  none15
## 15  med  optZ 2
## 16   A  low  none 7
## 17   A  low  optX 7
## 18   B  low  none10
## 19   B  low  optZ 8

One way to plot the data is to break-up the data into sub-frames, one
frame for each risk, order subjects by total counts, create dotplots,
and merge with trellis.c().  This almost works but in the merged plot I
cannot decrease column spacing to be small enough.  Also, the output of
trellis.c() would not work with useOuterStrips() which I really like.
My code is in TRY ONE below.

Another way to create the plot is specify y-limits for each panel with
custom prepanel and panel functions.  For each panel, the data-frame for
the panel row is isolated, subjects in the data-frame for the current
row are ordered by counts, panel y-limits are set to the re-ordered
levels, y-data for each panel is releveled, and data plotted with
standard panel.dotplot().  This somewhat works but lattice does not
honour the user-defined y-limits and labels are not correct.  I suspect
that it is not correct to use  y-relation=same in this case but free
and sliced do not give correct results too. My code in in TRY TWO
below.

If anybody can offer any assistance with this problem, it would be much
appreciated,
Sincerely,
Boris.


 BEGIN TRY ONE - MERGE LATTICE PLOTS 
library(lattice)
library(latticeExtra)
library(grid)

for (irisk in levels(df$risk)) {
  # subset data frame
  df.irisk - subset(df,risk==irisk)

  # order subjects by total count; store levels of subjectx variables
  # for later re-use in panel labels
  df.irisk$subjectx - df.irisk$subject[,drop=TRUE]
  df.irisk$subjectx - reorder(df.irisk$subjectx,df.irisk$count,sum)
  assign(paste('sbjx.',irisk,sep=''),levels(df.irisk$subjectx))

  # create dotplot and store it in oltc.{irisk} variable
  oltc.irisk - dotplot(subjectx~count|treatment,data=df.irisk,
layout=c(3,1),type=c('p','h'),
xlim=c(-1,16),origin=0,
xlab=,ylab=)
  assign(paste('oltc.',irisk,sep=''),oltc.irisk)
}

# combine everthing in one plot
oltc - c(low=oltc.low,med=oltc.med,high=oltc.high)
print(oltc)

# get rid of variable labels in middle and right column; decrease
# distance between columns.  But can't make inter-column spaces
# small enought and get rid of the panels in all but top rows.
laywid - trellis.par.get('layout.widths')
laywid$between - -5
laywid$axis.panel - 0.7
yscales - 

Re: [R] Plot a Matrix as an Image with ggplot

2013-02-15 Thread Alaios
You were right,
the following two work


testData-matrix(data=round(runif(25)),nrow=5,ncol=5,dimnames=list(Var1=c(1:5),Var2=c(1:5)))
 ggplot(melt(testData), aes(Var2,Var1, fill=value))+xlab(MHz) + 
ylab(Threshold) +    geom_raster() + scale_fill_gradient(low=#FF, 
high=#00)



still there are two minor unresolved issues.

a. How to reduce the extra space in the plot between legend and the tiles?
b. How to make color bar depict only the two values of the game, 0 and 1 , 
instead of 0, 0.25, 0.50, 0.75, 1?

I would like to thank you in advanec for your help
Regards
Alex


- Original Message -
From: Dennis Murphy djmu...@gmail.com
To: Alaios ala...@yahoo.com
Cc: 
Sent: Friday, February 15, 2013 10:19 AM
Subject: Re: [R] Plot a Matrix as an Image with ggplot

Your aesthetic is fill, not color. Change scale_color_gradient to
scale_fill_gradient and you'll get what you expect.

D.

On Thu, Feb 14, 2013 at 11:31 PM, Alaios ala...@yahoo.com wrote:


 Hi,
 thanks
 I changed slightly the code to be reproducible from everyone . I have  tried 
 ggplot
 but I need a bit of help to tweak it a bit

 So you can run the following in your computer

   testData-matrix(data=round(runif(25)),nrow=5,ncol=5,dimnames=list(1:5,1:5))

   p-ggplot(melt(testData), aes(Var2,Var1, fill=value))+xlab(MHz) + 
ylab(Threshold) +    geom_raster()


   p


 What I want to improve is
 a: make the colorbar so only two specific colors lets say black and white and 
 only two values 0 and 1 and somewhere the string as title of the color bar 
 text to appear.

 I have tried something like
 p+   scale_color_gradient(low=red,high=blue)
 Fehler in unit(tic_pos.c, mm) : 'x' and 'units' must have length  0


 ggplot(melt(testData), aes(Var2,Var1, fill=value,colorbin=2))+xlab(MHz) + 
 ylab(Threshold) +    geom_raster()
 but this did not affect colorbar entries.


 b. reduce/remove the grayish border that appears between the legend and the 
 image plot

 Could you please help me with these two?

 Regards
 Alex



 
 From: John Kane jrkrid...@inbox.com
 To: Alaios ala...@yahoo.com; R help r-help@r-project.org
 Sent: Thursday, February 14, 2013 4:35 PM
 Subject: RE: [R] Plot a Matrix as an Image with ggplot

 The R-help list is rather picky about what attached. None of your attachments 
 arrived.

 The str() info is useful but please supply some sample data

 The easiest way to supply data  is to use the dput() function.  Example with 
 your file named testfile:
 dput(testfile)

 Then copy the output and paste into your email.  For large data sets, you can 
 just supply a representative sample.  Usually,
 dput(head(testfile, 100)) will be sufficient.

 How are you writing the code/or what format is the original email. Your code 
 in the body of the text is badly messed up -- you probably need to post only 
 in text. HTML etc is automatically dropped.




 John Kane
 Kingston ON Canada


 -Original Message-
 From: ala...@yahoo.com
 Sent: Thu, 14 Feb 2013 07:15:05 -0800 (PST)
 To: r-help@r-project.org
 Subject: [R] Plot a Matrix as an Image with ggplot

 Dear all,
 I am trying to plot a matrix I have  as an image

 str(matrixToPlot)
  num [1:21, 1:66] 0 0 0 0 0 0 0 0 0 0 .


  that contains only 0s and 1s,


 where the xlabel will be Labeled as
 str(xLabel)
  num [1:66] 1e+09 1e+09 1e+09 1e+09 1e+09 ...

 and the yLabels will be labeled as
 str(yLabel)
  num [1:21] -88 -87 -86 -85 -84 -83 -82 -81 -80 -79 ...


 I have found on the internet that I can do something like that with
 ggplot2.
 For example
 you can run the following



 library(reshape2)library(ggplot2)m
 =matrix(rnorm(20),5)ggplot(melt(m),aes(Var1,Var2,fill=value))+geom_raster()

 What I see missing here is to get my matrix and transform it to a data
 frame with labels for x and y axis so as ggplot can print the values
 nice. If you get the idea this matrix will be printed as a two tile
 pattern let's say with black tiles zeros and white tiles black where a
 color bar will depicting that black means zero and white one.

 To help you with my code you will find attached the three items I have
 discussed here, the matrix, the x and the y labels.

 Could you please help me with that?

 Alex

 
 FREE 3D MARINE AQUARIUM SCREENSAVER - Watch dolphins, sharks  orcas on your 
 desktop!
 Check it out at http://www.inbox.com/marineaquarium

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

2013-02-15 Thread arun
Hi, 

I couldn't find the data(manaus) in library(Kendall). 
I guess data(GuelphP) is similar to your dataset. 

class(GuelphP) 
[1] ts 
 typeof(GuelphP) 
#[1] double 
 head(GuelphP) 
#[1] 0.47 0.51 0.35 0.19 0.33   NA 
 GuelphP 
#       Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec 
#1972 0.470 0.510 0.350 0.190 0.330    NA 0.365 0.650 0.825 1.000 0.385 0.900 
#1973 0.295 0.140 0.220 0.200 0.140 0.400    NA 0.495 1.100 0.590 0.270 0.300 
#1974    NA 0.065 0.240 0.058 0.079 0.065 0.120 0.091 0.058 0.120 0.120 0.110 
#1975 0.460 0.150 0.086 0.028    NA 0.110 0.360 0.180 0.065 0.130 0.120 0.190 
#1976 0.150 0.107 0.047 0.055 0.080 0.071 0.121 0.108 0.169 0.066 0.079 0.104 
#1977 0.157 0.140 0.070 0.056 0.042 0.116 0.106 0.094 0.097 0.050 0.079 0.114 

 
 
alb_data-read.table(text=
year    Jan    Feb    Mar    Apr    May    Jun    Jul    Aug Sep    Oct    Nov  
  Dec
 1997    NaN    NaN    NaN    NaN    NaN    NaN    NaN    NaN 23.406 16.166 
16.057 16.803
 1998 14.157 14.605 14.112 17.709 13.800 14.338 14.157 17.404 17.725 15.429 
16.090 18.061
 1999 14.888 13.837 15.929 13.637 16.020 14.699 15.987 15.212 14.752 15.935 
13.397 21.725
 2000 16.562 18.125 19.600 17.971 16.454 15.129 13.901 21.664 17.675 13.793 
13.464 16.452
2001 15.706 16.417 13.324 14.117 13.550 15.825 14.687 14.844 15.006 14.793 
13.489 15.726
 2002 11.777 11.775 11.564 13.141 14.462 15.517 18.801 15.652 17.127 14.070 
14.626 17.964 
,sep=,header=TRUE)
typeof(alb_data) 
#[1] list 
 class(alb_data) 
#[1] data.frame 

alb_data1-ts(matrix(unlist(t(alb_data[,-1])),ncol=1,byrow=F),frequency=12,start=1997)
 
#or
alb_data1- ts(as.vector(unlist(t(alb_data[,-1]))),frequency=12,start=1997)
 alb_data1
#    Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
#1997    NaN    NaN    NaN    NaN    NaN    NaN    NaN    NaN 23.406 16.166
#1998 14.157 14.605 14.112 17.709 13.800 14.338 14.157 17.404 17.725 15.429
#1999 14.888 13.837 15.929 13.637 16.020 14.699 15.987 15.212 14.752 15.935
#2000 16.562 18.125 19.600 17.971 16.454 15.129 13.901 21.664 17.675 13.793
#2001 15.706 16.417 13.324 14.117 13.550 15.825 14.687 14.844 15.006 14.793
#2002 11.777 11.775 11.564 13.141 14.462 15.517 18.801 15.652 17.127 14.070
 #   Nov    Dec
#1997 16.057 16.803
#1998 16.090 18.061
#1999 13.397 21.725
#2000 13.464 16.452
#2001 13.489 15.726
#2002 14.626 17.964

class(alb_data1) 
#[1] ts 
 typeof(alb_data1) 
#[1] double 
is.ts(alb_data1) 
#[1] TRUE 

library(Kendall)
SeasonalMannKendall(na.omit(alb_data1))
#tau = -0.143, 2-sided pvalue =0.20287

A.K. 



- Original Message -
From: twynne timothy.wy...@noaa.gov
To: r-help@r-project.org
Cc: 
Sent: Friday, February 15, 2013 2:14 PM
Subject: [R] convert list into a time series

I am trying to use the SeasonalMannKendall function in the Kendall package.
My dataset (alb_data) is in the same format as the example dataset 
(manaus) in the package.

 class(manaus)
[1] ts
 is.ts(manaus)
[1] TRUE
 typeof(manaus)
[1] double
 alb_data=read.table(R:/albemarle_manken.txt, header=T)
 head(alb_data)
   year    Jan    Feb    Mar    Apr    May    Jun    Jul    Aug Sep    
Oct    Nov    Dec
1 1997    NaN    NaN    NaN    NaN    NaN    NaN    NaN    NaN 23.406 
16.166 16.057 16.803
2 1998 14.157 14.605 14.112 17.709 13.800 14.338 14.157 17.404 17.725 
15.429 16.090 18.061
3 1999 14.888 13.837 15.929 13.637 16.020 14.699 15.987 15.212 14.752 
15.935 13.397 21.725
4 2000 16.562 18.125 19.600 17.971 16.454 15.129 13.901 21.664 17.675 
13.793 13.464 16.452
5 2001 15.706 16.417 13.324 14.117 13.550 15.825 14.687 14.844 15.006 
14.793 13.489 15.726
6 2002 11.777 11.775 11.564 13.141 14.462 15.517 18.801 15.652 17.127 
14.070 14.626 17.964
 class(alb_data)
[1] data.frame
 is.ts(alb_data)
[1] FALSE
 typeof(alb_data)
[1] list


Can anyone please help me get my dataset into the correct format so that 
I can run the test?
Thank you for your time and help!




--
View this message in context: 
http://r.789695.n4.nabble.com/convert-list-into-a-time-series-tp4658716.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-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] 2 setGeneric's, same name, different method signatures

2013-02-15 Thread Greg Minshall
William,

 and here mention intra-class name clashes (I'm not sure what you
 mean by this).

sorry, i meant, in something like C++, if i have a class Foo and you
have a class Bar, then i can invent whatever method names/signatures i
want, independent of whatever method names/signatures *you* want.  *i*
just need to make sure i don't introduce incompatible methods with the
same name/signature *within* Foo.


cheers, Greg

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

2013-02-15 Thread Tanushree Tunstall
Hello All,
I am new to the list. I have been learning to use R recently and its been
great to see so much help available.

However I must admit, I have stumbled upon a problem with correlation
matrices and was hoping if someone could  help.

I have an excel workbook with a sheet per drug.
Each sheet has 40 patients' drug response measured over different time
points (several days).

I have already been able to create a cor() matrix of all the drugs for all
the patients on each day. However I have been asked to prepare a matrix of
all days of all people for all drugs. So since there are 26 drugs, I am
expected to form a  26x26 matrix for all this data.

I am not even sure if this is possible. Can anyone point me in the right
direction?


PS: If I make individual matrix and combine them using cbind(), it is not a
26x26 sqaure matrix.

Any help would be greatly appreciated.

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] Iterating through slots of an S4 object

2013-02-15 Thread Scott Robinson
I want to loop through slots of an S4 object and am unsure how to do so
since the only way I can find to access them is individually in the form
'object@slotName'. I have guessed at a few possibilities which did not work
and I have read some S4 object tutorials and things but still unsuccessful.

I presume it is possible though?

Any help would be much appreciated,

Scott

[[alternative HTML version deleted]]

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


Re: [R] data formatting

2013-02-15 Thread arun


HI Eliza,

Suppose you have 147 data files in the same working directory.   Here, I am 
using Eliza1.txt and a modified Eliza2.txt (attached).
list.files()
#[1] Eliza1.txt Eliza2.txt

lapply(list.files(),function(i) str_count(gsub( $,,readLines(i)), )) 
#count the spaces.  Used gsub as there were spaces at the end (possibly due to 
formatting error) #which was removed.  If there are no spaces at the end, you 
don't need ?gsub()
#[[1]]
 #[1] 7 7 7 7 6 7 7 7 7 7 6 6 7 7 7 7 6 7 7 7 7 7 6 6
#
#[[2]]
# [1] 7 7 7 7 6 7 7 7 7 7 6 6 7 7 7 7 6 7 7 7 7 7 6 6


res- lapply(list.files(),function(i) {Lines2-gsub( 
$,,readLines(i));Lines2[str_count(Lines2, )==7]- 
str_replace(Lines2[str_count(Lines2, )==7],\\s+, 
);Lines2[str_count(Lines2, )==6]- str_replace(Lines2[str_count(Lines2, 
)==6],\\s+,    
);substr(Lines2[substr(Lines2,6,6)==0|substr(Lines2,9,9)==0],6,6)- 
;substr(Lines2[substr(Lines2,6,6)==0|substr(Lines2,9,9)==0],9,9)- ;Lines2})

names(res)-gsub(\\..*,,list.files())
res
#$Eliza1
# [1] 1911. 1. 1 7.87 1911. 1. 2 9.26 1911. 1. 3 8.06
# [4] 1911. 1. 4 8.13 1911. 1. 5    12.90 1911. 2. 6 5.45
# [7] 1911. 2. 7 3.26 1911. 3. 8 5.70 1911. 3. 9 9.24
#[10] 1911. 4.10 7.60 1911. 5.11    14.82 1911. 5.12    14.10
#[13] 1911. 6.13 7.87 1911. 6.14 9.26 1911. 7.15 8.06
#[16] 1911. 7.16 8.13 1911. 8.17    12.90 1911. 8.18 5.45
#[19] 1911. 9.19 3.26 1911. 9.20 5.70 1911.10.21 9.24
#[22] 1911.10.22 7.60 1911.11.23    14.82 1911.12.24    14.10

#$Eliza2
# [1] 1911. 1. 1 4.87  1911. 1. 2 11.26 1911. 1. 3 6.06 
# [4] 1911. 1. 4 8.13  1911. 1. 5    11.90  1911. 2. 6 5.55 
# [7] 1911. 2. 7 3.16  1911. 3. 8 5.10  1911. 3. 9 9.34 
#[10] 1911. 4.10 7.10  1911. 5.11    14.92  1911. 5.12    14.20 
#[13] 1911. 6.13 7.77  1911. 6.14 9.36  1911. 7.15 8.66 
#[16] 1911. 7.16 8.23  1911. 8.17    11.90  1911. 8.18 15.45
#[19] 1911. 9.19 13.26 1911. 9.20 15.77 1911.10.21 19.34
#[22] 1911.10.22 7.66  1911.11.23    14.84  1911.12.24    14.11 
 lapply(res,function(x) str_count(x, ))
#$Eliza1
# [1] 7 7 7 7 6 7 7 7 7 6 5 5 6 6 6 6 5 6 6 6 5 5 4 4

#$Eliza2
# [1] 7 7 7 7 6 7 7 7 7 6 5 5 6 6 6 6 5 6 6 6 5 5 4 4
Hope this helps.
A.K.








From: eliza botto eliza_bo...@hotmail.com
To: smartpink...@yahoo.com smartpink...@yahoo.com 
Sent: Friday, February 15, 2013 4:47 PM
Subject: RE: data formatting



Thankyou very much for replying arun. i just need to know, what change will i 
have  to make if i am importing 147 data files into a list. what difference 
will it make on  the first command which is,
 Lines1-readLines(textConnection(1911.01.01   7.87
 1911.01.02   9.26 
 1911.01.03   8.06 
 1911.01.04   8.13 
 1911.01.05  12.90 
 1911.02.06   5.45 
 1911.02.07   3.26 
 1911.03.08   5.70 
 1911.03.09   9.24 
 1911.04.10   7.60 
 1911.05.11  14.82 
 1911.05.12  14.10 
 1911.06.13   7.87 
 1911.06.14   9.26 
 1911.07.15   8.06 
 1911.07.16   8.13 
 1911.08.17  12.90 
 1911.08.18   5.45 
 1911.09.19   3.26 
 1911.09.20   5.70 
 1911.10.21   9.24 
 1911.10.22   7.60 
 1911.11.23  14.82 
 1911.12.24  14.10)) 

thankyou so very much...

elisa


 Date: Fri, 15 Feb 2013 11:11:36 -0800
 From: smartpink...@yahoo.com
 Subject: Re: data formatting
 To: eliza_bo...@hotmail.com
 CC: r-help@r-project.org
 
 
 
 Dear Eliza,
 
 Try this:
 
 Lines1-readLines(textConnection(1911.01.01   7.87
 1911.01.02   9.26 
 1911.01.03   8.06 
 1911.01.04   8.13 
 1911.01.05  12.90 
 1911.02.06   5.45 
 1911.02.07   3.26 
 1911.03.08   5.70 
 1911.03.09   9.24 
 1911.04.10   7.60 
 1911.05.11  14.82 
 1911.05.12  14.10 
 1911.06.13   7.87 
 1911.06.14   9.26 
 1911.07.15   8.06 
 1911.07.16   8.13 
 1911.08.17  12.90 
 1911.08.18   5.45 
 1911.09.19   3.26 
 1911.09.20   5.70 
 1911.10.21   9.24 
 1911.10.22   7.60 
 1911.11.23  14.82 
 1911.12.24  14.10)) 
 
 Lines2-Lines1[Lines1!=]
 library(stringr)
  str_count(Lines2,  )
 # [1] 7 7 7 7 6 7 7 7 7 7 6 6 7 7 7 7 6 7 7 7 7 7 6 6
 
 
 Lines2[str_count(Lines2, )==7]- str_replace(Lines2[str_count(Lines2, 
 )==7],\\s+, ) #reduced 2 spaces
 
  Lines2[str_count(Lines2, )==6]- str_replace(Lines2[str_count(Lines2, 
 )==6],\\s+,    ) #reduced 2 spaces
  str_count(Lines2, )
 # [1] 5 5 5 5 4 5 5 5 5 5 4 4 5 5 5 5 4 5 5 5 5 5 4 4
 substr(Lines2[substr(Lines2,6,6)==0|substr(Lines2,9,9)==0],6,6)- 
 substr(Lines2[substr(Lines2,6,6)==0|substr(Lines2,9,9)==0],9,9)- 
 str_count(Lines2, ) #see the difference in space.  This counts all the 
 space.  Here 2 white space are added to replace 0
 # [1] 7 7 7 7 6 7 7 7 7 6 5 5 6 6 6 6 5 6 6 6 5 5 4 4
 Lines2
 # [1] 1911. 1. 1 7.87 1911. 1. 2 9.26 1911. 1. 3 8.06
 # [4] 1911. 1. 4 8.13 1911. 1. 5    12.90 1911. 2. 6 

Re: [R] How can I plot graphs together?

2013-02-15 Thread soon yi
look at dev.new() to specify plot window size

and then ?layout to specify number and size of each plot in the window




Jiaqi.Zhang wrote
 Hi, all,
 
 I am working on the following code to learn how to plot graphs together. I
 used the par(mfrow=c(1,3)) function to try to put all three plot() graphs
 together. But it always fail without any error message? Can anybody help
 me out?
 
 
 ## Synthetic Data Generation
 
 
 n = 500
 x1 = rnorm(n, 1, 100)
 x2 = rnorm(n, 10, 100)
 x3 = rnorm(n, 5, 1)
 x4 = rnorm(n, 10, 10)
 x5 = rbinom(n, 1, .4)
 x6 = rnorm(n, 30, 5)
 treatment = rbinom(n, 1, .15)
 data = cbind(x1,x2,x3,x4,x5,x6, treatment)
 dim(data)
 
 
 ## Propensity score matching
 ## nearest neighbor matching (1:1)
 
 
 require(MatchIt)
 # data1 is the subset of data with only the selected variables mentioned
 below
 data1 = data[,c(x1,x2,x3,x4,x5,x6, treatment)]
 # getting rid of missing values (below)
 data1 = as.data.frame(na.omit(data1)) 
 # matching is performed below using propensity scores given the covariates
 mentioned below
 m.out = matchit(treatment~x1+x2+x3+x4+x5+x6,method=nearest, data=data1,
 ratio = 1)
 # check the sample sizes (below)
 m.out 
 # Final matched data saved as final_data
 final_data = match.data(m.out) 
 # (here distance = propensity score)
 # check balance (below)
 par(mfrow=c(1,3))
 plot(m.out) # covariate balance
 plot(m.out, type = jitter) # propensity score locations
 plot(m.out, type = hist) #check matched treated vs matched control
 
 I hope to put 
 plot(m.out), plot(m.out, type = jitter), and plot(m.out, type = hist)
 in one graph.





--
View this message in context: 
http://r.789695.n4.nabble.com/How-can-I-plot-graphs-together-tp4658612p4658742.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] Iterating through slots of an S4 object

2013-02-15 Thread Duncan Murdoch

On 13-02-15 6:14 PM, Scott Robinson wrote:

I want to loop through slots of an S4 object and am unsure how to do so
since the only way I can find to access them is individually in the form
'object@slotName'. I have guessed at a few possibilities which did not work
and I have read some S4 object tutorials and things but still unsuccessful.


See ?slotNames.  For example,

for (slot in slotNames(x)) {
  cat(slot, :\n)
  print(slot(x, slot))
}

Duncan Murdoch



I presume it is possible though?

Any help would be much appreciated,

Scott

[[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] Plot a Matrix as an Image with ggplot

2013-02-15 Thread Dennis Murphy
Hi:

See if the following works for you:

library(reshape2)
library(ggplot2)
tdm - melt(testData)

ggplot(tdm, aes(x = Var2, y = Var1, fill = factor(value))) +
labs(x = MHz, y = Threshold, fill = Value) +
geom_raster() +
scale_fill_manual(breaks = levels(factor(tdm$value)),
  values = c(white, black)) +
theme(plot.background = element_rect(fill = grey90),
  legend.background = element_rect(fill = grey90)) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0))

It was necessary to create a different plot background color because
you wanted to remove the padding, which blended much of the plot
boundary with the original white background, and had to do something
similar to the legend background so that you could see the white
legend box.

The expand = c(0, 0) argument to the two scale functions answers your
second question, and converting value from numeric to vector answers
the first.

Dennis

On Fri, Feb 15, 2013 at 3:24 PM, Alaios ala...@yahoo.com wrote:
 You were right,
 the following two work


 testData-matrix(data=round(runif(25)),nrow=5,ncol=5,dimnames=list(Var1=c(1:5),Var2=c(1:5)))
  ggplot(melt(testData), aes(Var2,Var1, fill=value))+xlab(MHz) + 
 ylab(Threshold) +geom_raster() + scale_fill_gradient(low=#FF, 
 high=#00)



 still there are two minor unresolved issues.

 a. How to reduce the extra space in the plot between legend and the tiles?
 b. How to make color bar depict only the two values of the game, 0 and 1 , 
 instead of 0, 0.25, 0.50, 0.75, 1?

 I would like to thank you in advanec for your help
 Regards
 Alex


 - Original Message -
 From: Dennis Murphy djmu...@gmail.com
 To: Alaios ala...@yahoo.com
 Cc:
 Sent: Friday, February 15, 2013 10:19 AM
 Subject: Re: [R] Plot a Matrix as an Image with ggplot

 Your aesthetic is fill, not color. Change scale_color_gradient to
 scale_fill_gradient and you'll get what you expect.

 D.

 On Thu, Feb 14, 2013 at 11:31 PM, Alaios ala...@yahoo.com wrote:


 Hi,
 thanks
 I changed slightly the code to be reproducible from everyone . I have  tried 
 ggplot
 but I need a bit of help to tweak it a bit

 So you can run the following in your computer

   
 testData-matrix(data=round(runif(25)),nrow=5,ncol=5,dimnames=list(1:5,1:5))

   p-ggplot(melt(testData), aes(Var2,Var1, fill=value))+xlab(MHz) + 
 ylab(Threshold) +geom_raster()


   p


 What I want to improve is
 a: make the colorbar so only two specific colors lets say black and white 
 and only two values 0 and 1 and somewhere the string as title of the color 
 bar text to appear.

 I have tried something like
 p+   scale_color_gradient(low=red,high=blue)
 Fehler in unit(tic_pos.c, mm) : 'x' and 'units' must have length  0


 ggplot(melt(testData), aes(Var2,Var1, fill=value,colorbin=2))+xlab(MHz) + 
 ylab(Threshold) +geom_raster()
 but this did not affect colorbar entries.


 b. reduce/remove the grayish border that appears between the legend and the 
 image plot

 Could you please help me with these two?

 Regards
 Alex



 
 From: John Kane jrkrid...@inbox.com
 To: Alaios ala...@yahoo.com; R help r-help@r-project.org
 Sent: Thursday, February 14, 2013 4:35 PM
 Subject: RE: [R] Plot a Matrix as an Image with ggplot

 The R-help list is rather picky about what attached. None of your 
 attachments arrived.

 The str() info is useful but please supply some sample data

 The easiest way to supply data  is to use the dput() function.  Example with 
 your file named testfile:
 dput(testfile)

 Then copy the output and paste into your email.  For large data sets, you 
 can just supply a representative sample.  Usually,
 dput(head(testfile, 100)) will be sufficient.

 How are you writing the code/or what format is the original email. Your code 
 in the body of the text is badly messed up -- you probably need to post only 
 in text. HTML etc is automatically dropped.




 John Kane
 Kingston ON Canada


 -Original Message-
 From: ala...@yahoo.com
 Sent: Thu, 14 Feb 2013 07:15:05 -0800 (PST)
 To: r-help@r-project.org
 Subject: [R] Plot a Matrix as an Image with ggplot

 Dear all,
 I am trying to plot a matrix I have  as an image

 str(matrixToPlot)
  num [1:21, 1:66] 0 0 0 0 0 0 0 0 0 0 .


  that contains only 0s and 1s,


 where the xlabel will be Labeled as
 str(xLabel)
  num [1:66] 1e+09 1e+09 1e+09 1e+09 1e+09 ...

 and the yLabels will be labeled as
 str(yLabel)
  num [1:21] -88 -87 -86 -85 -84 -83 -82 -81 -80 -79 ...


 I have found on the internet that I can do something like that with
 ggplot2.
 For example
 you can run the following



 library(reshape2)library(ggplot2)m
 =matrix(rnorm(20),5)ggplot(melt(m),aes(Var1,Var2,fill=value))+geom_raster()

 What I see missing here is to get my matrix and transform it to a data
 frame with labels for x and y axis so as ggplot can print the values
 nice. If you get the idea this matrix will be printed as a two tile
 

[R] odd behavior within R2HTML

2013-02-15 Thread Erin Hodgess
Dear R People:

I'm using R2HTML but having a strange result.

Here is the original data:

resp trt block
90.3 A I
89.2 A II
98.2 A III
93.9 A IV
87.4 A V
97.9 A VI
92.5 B I
89.5 B II
90.6 B III
94.7 B IV
87.0 B V
95.8 B VI
85.5 C I
90.8 C II
89.6 C III
86.2 C IV
88.0 C V
93.4 C VI
82.5 D I
89.5 D II
85.6 D III
87.4 D IV
78.9 D V
90.7 D VI

And here are the commands:
 resin1.df - read.table(resin1.txt,header=TRUE)
 #set up R-B anova
 resin1.aov - aov(resp ~ trt + block, data=resin1.df)
 library(R2HTML
+ )
 HTMLStart(outdir=getwd(),file=resin2,echo=TRUE)

 *** Output redirected to directory:  c:/R64/R-2.15.2/bin/x64
 *** Use HTMLStop() to end redirection.[1] TRUE
HTML library(effects)
Loading required package: lattice
Loading required package: grid
Loading required package: MASS
Loading required package: nnet
Loading required package: colorspace

Attaching package: ‘effects’

The following object(s) are masked from ‘package:datasets’:

Titanic

HTML summary(allEffects(resin1.aov)
+ )
 model: resp ~ trt + block

 trt effect
trt
   ABCD
92.81667 91.68333 88.91667 85.76667

 Lower 95 Percent Confidence Limits
trt
   ABCD
90.46148 89.32815 86.56148 83.41148

 Upper 95 Percent Confidence Limits
trt
   ABCD
95.17185 94.03852 91.27185 88.12185

 block effect
block
 I IIIII IV  V VI
87.700 89.750 91.000 90.550 85.325 94.450

 Lower 95 Percent Confidence Limits
block
  I  II III  IV   V  VI
84.8155 86.8655 88.1155 87.6655 82.4405 91.5655

 Upper 95 Percent Confidence Limits
block
  I  II III  IV   V  VI
90.5845 92.6345 93.8845 93.4345 88.2095 97.3345
HTML HTMLStop()
[1] c:/R64/R-2.15.2/bin/x64/resin2_main.html
 file.show(resin2_main.html)


When I look in the resin2_main.html file, the results from the summary
are not there.

Has anyone run into this before, please?

Thanks for any help!

Sincerely,
Erin
R-2.15.2, Windows, 64-bit


-- 
Erin Hodgess
Associate Professor
Department of Computer and Mathematical Sciences
University of Houston - Downtown
mailto: erinm.hodg...@gmail.com

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


[R] Interpret R-squared and cor in R

2013-02-15 Thread Janesh Devkota
Hi I am trying to find the relationship between two variables.

First I fitted a linear model between two variables and I found the
following results:
Residual standard error: 0.03253 on 2498 degrees of freedom
Multiple R-squared: 0.5551, Adjusted R-squared: 0.5549
F-statistic:  3116 on 1 and 2498 DF,  p-value:  2.2e-16

Then I used the cor function to see the correlation between two variable
I get the following result
-0.7450344

How can we interpret the result based on R-squared and correlation ? From
the p-value we can see that there is very strong relationship between
variables as it is  way less that 0.001

Can anyone kindly explain the difference between Multiple R squared,
adjusted R-squared and correlation and how to report these values while
writing a report ?

Thank you so much.

[[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] Interpret R-squared and cor in R

2013-02-15 Thread Jeremy Miles
On 15 February 2013 21:26, Janesh Devkota janesh.devk...@gmail.com wrote:

 Hi I am trying to find the relationship between two variables.

 First I fitted a linear model between two variables and I found the
 following results:
 Residual standard error: 0.03253 on 2498 degrees of freedom
 Multiple R-squared: 0.5551, Adjusted R-squared: 0.5549
 F-statistic:  3116 on 1 and 2498 DF,  p-value:  2.2e-16

 Then I used the cor function to see the correlation between two variable
 I get the following result
 -0.7450344


r is a correlation (it actually stands for regression).

R (upper case) is a multiple correlation. But you only have one predictor,
so it's a correlation.

R squared is R (or r), squared.  So -0.7450433^2 = 0.555.



 How can we interpret the result based on R-squared and correlation ? From
 the p-value we can see that there is very strong relationship between
 variables as it is  way less that 0.001



The p-value doesn't tell you about the strength of the relationship.


 Can anyone kindly explain the difference between Multiple R squared,
 adjusted R-squared and correlation and how to report these values while
 writing a report ?


I can suggest a number of books that do this much better than I could in an
email. But you probably have a favorite of your own.

Jeremy

[[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] subplot (Hmisc) and radial.plot (plotrix) problem

2013-02-15 Thread Andrew Roberts

Folks,

I am having problems with a plot I want to create to give an impression 
of changes in an ordinal scale measure (1-5)  at three time points  (0, 
14 and 21 days). I can produce a radial plot of bare vectors but getting 
this to appear on the base plot is not possible as it always seems to 
end up below the plot area and even outside the plot window.


It seems I have not understood something about the coordinate system, Is 
anyone able to shed any light on this?


many thanks,

Andrew Roberts

# Example ###
library(Hmisc)
library(plotrix)
library(plyr)
# Create a data frame of ordinal results #
X -c(1,3,4,5,3,2,4,4,5,2,1,3,3,2,4)
Y -c(1,2,2,4,1,1,2,3,3,1,2,2,4,1,2)
Z -c(1,2,1,2,1,1,3,2,2,2,1,2,2,1,1)
data -as.data.frame(cbind(X,Y,Z))
data$xydiff -data$X - data$Y
data$yzdiff -data$Y - data$Z

# Create the background frame for the subplots #
plot(1, type=n, xlim=c(0,21),
ylim = c(1,5),
main=Score Changes,
xlab=Days , ylab=Ashworth Scale)
# Work through each level of the scale #
for(i in 1:5){
  t  -subset(data, X==i)
  t2 -count(t,vars=xydiff)
  t2$ang -tan(t2$xydiff/14)
  ang -as.matrix(t2)[,3]
  freq -as.matrix(t2)[,2]
  subplot(radial.plot(freq,ang, labels= , clockwise=TRUE, 
show.grid=FALSE,radial.lim=c(0,8),

point.symbols=12,point.col=green, show.centroid=TRUE,
show.grid.labels=TRUE, show.radial.grid=FALSE, 
line.col=red,

lwd=3,add=TRUE),
  0,i)
  rm(ang,freq,t,t2)
}
# But this exmple works! ##
subplot( hist(rnorm(100)), 1, 4)

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


Re: [R] Interpret R-squared and cor in R

2013-02-15 Thread David Winsemius


On Feb 15, 2013, at 9:26 PM, Janesh Devkota wrote:


Hi I am trying to find the relationship between two variables.

First I fitted a linear model between two variables and I found the
following results:
Residual standard error: 0.03253 on 2498 degrees of freedom
Multiple R-squared: 0.5551, Adjusted R-squared: 0.5549
F-statistic:  3116 on 1 and 2498 DF,  p-value:  2.2e-16

Then I used the cor function to see the correlation between two  
variable

I get the following result
-0.7450344

How can we interpret the result based on R-squared and correlation ?  
From

the p-value we can see that there is very strong relationship between
variables as it is  way less that 0.001

Can anyone kindly explain the difference between Multiple R squared,
adjusted R-squared and correlation and how to report these values  
while

writing a report ?


This is not an on-topic question for this mailing list. The  
CrossValidated website is more likely to respond as you might have  
wished. (Please read the Posting Guide.)


--

David Winsemius, MD
Alameda, CA, USA

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