[R] timestamp shifted by hour(s) while mering zoo objects

[Keith]

Dear R-users,

I have two regular hourly time series data which were recorded in time 
zone GMT+1, and now I would like to merge them together for further 
analyses. Here I used zoo and merge.zoo for my purposes and everything 
worked fine except the timestamp shifted 2 hours after merging which 
bugs me a little bit. Here is the example:


data01
00:00:00 01.01.2007, 8.0250
01:00:00 01.01.2007, 8.0167
02:00:00 01.01.2007, 10.0917
03:00:00 01.01.2007, 8.6750
04:00:00 01.01.2007, 6.3250

data02
00:00:00 01.01.2007, 257.58
01:00:00 01.01.2007, 239.92
02:00:00 01.01.2007, 234.00
03:00:00 01.01.2007, 220.00
04:00:00 01.01.2007, 206.92

which are both read into zoo object, data01 and data02, separately by 
setting tz = "GMT+1". However, while merging function is operated, the 
result is


merge.zoo(data01, data02)
 data01 data02
2007-01-01 02:00:00  8.0250 257.58
2007-01-01 03:00:00  8.0167 239.92
2007-01-01 04:00:00 10.0917 234.00
2007-01-01 05:00:00  8.6750 220.00
2007-01-01 06:00:00  6.3250 206.92

which is 2 hours shifted comparing to the original data. I am wondering 
if it's the problem of tz parameter. Hence, I re-read the data by 
setting tz = "GMT", and the merging result is


merge.zoo(data01, data02)
 data01 data02
2007-01-01 01:00:00  8.0250 257.58
2007-01-01 02:00:00  8.0167 239.92
2007-01-01 03:00:00 10.0917 234.00
2007-01-01 04:00:00  8.6750 220.00
2007-01-01 05:00:00  6.3250 206.92

which is 1 hour shifted. I only noticed this but don't know why and how 
to fix it. Does anyone have idea about this issue?


Best regards,
Keith

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

Re: [R] timestamp shifted by hour(s) while mering zoo objects

[Keith]

Thanks Gabor,

Just one comment on the error you received on the Vista machine. I also 
noticed this problem on Windows as well, and found there is no 
specification of "GMT+1" on Windows XP. The closest setting for "GMT+1" 
on windows is "Africa/Lagos", however there is no this problem on *unix 
machine. This is system dependent problem, and I have to admit it's 
somehow annoying.



Gabor Grothendieck wrote:

On my Vista system I get an error message saying that there is no such
timezone as GMT+1.

You may be better off not using time zones and just adjusting the
times yourself.  You could just use chron and avoid the entire time
zone problem in the first place.

In the first of the two approaches below we simply use GMT everywhere
in the first approach and adjust the result by one hour manually.

In the second approach we use chron which has no time zones and again
just adjust the times as we need.

We have added one hour in both cases here but you may need to add or
subtract a different number.

Sys.setenv(TZ = "GMT")

library(zoo)

data01.Lines <- "junk
00:00:00 01.01.2007, 8.0250
01:00:00 01.01.2007, 8.0167
02:00:00 01.01.2007, 10.0917
03:00:00 01.01.2007, 8.6750
04:00:00 01.01.2007, 6.3250"

data02.Lines <- "junk
00:00:00 01.01.2007, 257.58
01:00:00 01.01.2007, 239.92
02:00:00 01.01.2007, 234.00
03:00:00 01.01.2007, 220.00
04:00:00 01.01.2007, 206.92"

# POSIXct

data01 <- read.zoo(textConnection(data01.Lines), sep=",",
format="%H:%M:%S %d.%m.%Y", tz = "GMT", strip.white=TRUE, skip=1)

data02 <- read.zoo(textConnection(data02.Lines), sep=",",
format="%H:%M:%S %d.%m.%Y", tz = "GMT", strip.white=TRUE, skip=1)

m <- merge(data01, data02)

# add 3600 seconds to each time
time(m) <- time(m) + 3600
m

###

# chron

library(chron)

data01c <- read.zoo(textConnection(data01.Lines), sep=",",
format="%H:%M:%S %d.%m.%Y", FUN = as.chron, strip.white=TRUE, skip=1)

data02c <- read.zoo(textConnection(data02.Lines), sep=",",
format="%H:%M:%S %d.%m.%Y", FUN = as.chron, strip.white=TRUE, skip=1)

mc <- merge(data01c, data02c)

time(mc) <- time(mc) + as.numeric(times("01:00:00"))
mc


On Sun, Mar 7, 2010 at 3:36 PM, Keith  wrote:

Thanks Gabor,

You're right. The problem comes from the environment variable TZ. I just
tried the Sys.getenv("TZ") and it's nothing there. After I have set the
environment variable TZ as the same as the data, let's say
Sys.setenv(TZ="GMT+1"), the problem is gone.

In order to complete the problem I've mentioned, here are the data and the
code:
data01.txt
% time[GMT+1:00] Temperature[°C]
00:00:00 01.01.2007, 8.0250
01:00:00 01.01.2007, 8.0167
02:00:00 01.01.2007, 10.0917
03:00:00 01.01.2007, 8.6750
04:00:00 01.01.2007, 6.3250

data02.txt
% time[GMT+1:00] Conductance[µS]
00:00:00 01.01.2007, 257.58
01:00:00 01.01.2007, 239.92
02:00:00 01.01.2007, 234.00
03:00:00 01.01.2007, 220.00
04:00:00 01.01.2007, 206.92

data01 <- read.zoo("data01.txt", sep=",", format="%H:%M:%S %d.%m.%Y",
tz="GMT+1", strip.white=TRUE, skip=1)

data02 <- read.zoo("data02.txt", sep=",", format="%H:%M:%S %d.%m.%Y",
tz="GMT+1", strip.white=TRUE, skip=1)

merge.zoo(data01, data02)

Besides, thanks for your recommendation, and I'll have a check the R News
4/1.

Regards,
Keith

Gabor Grothendieck wrote:

Without reproducible code (that means we can copy your code from your
post, paste it into our session and see the same problem that you see)
there is not much that can be said that addresses your specific
situation but in terms of general advice:

- the inappropriate use of time zones is a frequent source of errors
in R for the unwary and you should read R News 4/1 to find out more
about this.

- if after reading that you still want to use POSIXct your best bet is
to set the time zone of your session to GMT and work entirely in GMT:
Sys.setenv(TZ = "GMT")


On Sun, Mar 7, 2010 at 12:20 PM, Keith  wrote:

Dear R-users,

I have two regular hourly time series data which were recorded in time
zone
GMT+1, and now I would like to merge them together for further analyses.
Here I used zoo and merge.zoo for my purposes and everything worked fine
except the timestamp shifted 2 hours after merging which bugs me a little
bit. Here is the example:

data01
00:00:00 01.01.2007, 8.0250
01:00:00 01.01.2007, 8.0167
02:00:00 01.01.2007, 10.0917
03:00:00 01.01.2007, 8.6750
04:00:00 01.01.2007, 6.3250

data02
00:00:00 01.01.2007, 257.58
01:00:00 01.01.2007, 239.92
02:00:00 01.01.2007, 234.00
03:00:00 01.01.2007, 220.00
04:00:00 01.01.2007, 206.92

which are both read into zoo object, data01 and data02, separately by
setting tz = "GMT+1". However, while merging

[R] Error in pruning hclust object using maptree package

[Keith]

Dear R users,

Due to too many children in my clustering result, I would like to prune 
the tree-like hclust object into a certain groups. However, an error:


> clip.clust: no data provided for hclust object

always shows up.

Firstly, I tried the example in the maptree package, and it worked well. 
Then, I tried to used the example in stats package to generate a hclust 
object and tried the clip.clust method again, and the error still 
occurred. Here is the example:


> hc <- hclust(dist(USArrests), "ave")
> pr <- clip.clust(hc, k=4)
Error in clip.clust(hc, k = 4) :
  clip.clust: no data provided for hclust object

Does anyone have the idea? My environment is
R version 2.10.1 (2009-12-14)
x86_64-pc-linux-gnu
[1] maptree_1.4-6

Best regards,
Keith

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

Re: [R] Error in pruning hclust object using maptree package

[Keith]

My bad.

I did read the help. However, I didn't pay enough attention to the data 
argument, and it seems it doesn't need to use data argument in the 
earlier version.


I also noticed this argument afterwards. Now it works for me as well. 
Thanks a lot.


Uwe Ligges wrote:



On 15.03.2010 16:37, Keith wrote:

Dear R users,

Due to too many children in my clustering result, I would like to prune
the tree-like hclust object into a certain groups. However, an error:

 > clip.clust: no data provided for hclust object

always shows up.

Firstly, I tried the example in the maptree package, and it worked well.
Then, I tried to used the example in stats package to generate a hclust
object and tried the clip.clust method again, and the error still
occurred. Here is the example:

 > hc <- hclust(dist(USArrests), "ave")
 > pr <- clip.clust(hc, k=4)
Error in clip.clust(hc, k = 4) :
clip.clust: no data provided for hclust object


Have you read ?clip.clust ??
It says for argument "data":
 data clustered dataset for hclust application.

whic means you need to provide "data" as the error message suggests as 
well.


 pr <- clip.clust(hc, k=4, data=USArrests)
 pr$membership

works for me.

Uwe Ligges








Does anyone have the idea? My environment is
R version 2.10.1 (2009-12-14)
x86_64-pc-linux-gnu
[1] maptree_1.4-6

Best regards,
Keith

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Plotting hourly time-series data loaded from file using plot.ts

[Keith]

Hello everyone,

I am just a tyro in R and would like your kindly help for some
problems which I've been struggling for a while but still in vain.

I have a time-series file (with some missing value ) which looks like

time[sec] , Factor1 , Factor2
00:00:00 01.01.2007 , 0. , 0.176083
01:00:00 01.01.2007 , 0. , 0.176417
[ ... ]
11:00:00 10.06.2007 , 0. , 0.148250
12:00:00 10.06.2007 , NA , 0.147000
13:00:00 10.06.2007 , NA , 0.144417
[ ... ]

and I would like to do some basic time-series analyses using R. The
first idea is to plot these time-series events and the main problem
was the handling of the date/time format in the 1st column. I was
using the script below to deal with:

data <- 
read.table("file",header=TRUE,sep=",",colClasses=c("character","numeric","numeric"))
data$time.sec. <- as.POSIXct(data$time.sec.,format="%H:%M:%S %d.%m.%Y")
dataTs <- as.ts(data)
plot.ts(dataTs)

Then, the plot showed up with 3 subplots in one plot. The 1st is the
linear line with the x-axis being just the sequence of orders and
y-axis being wrong numbers which is completely wrong. The 2nd and the
3rd are correct but the x-axis is still wrong. Does anyone know how to
plot correct Factor1 and Factor2 with respect to the 1st column time
format? Or, probably should I use some other packages? Besides, how
can I plot these two time-series data (Factor1 and Factor2) in two
separate plots?

Best regards,
Keith

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

Re: [R] Plotting hourly time-series data loaded from file using plot.ts

[Keith]

Thanks Gabor,

I tried a little bit, and your example works. However, it seems that
the read.zoo has a limitation of records up to around 300 !? I took
your suggestion and modified a little bit in order to read from the
file which contains about 9000 records:

dataTs <- read.zoo("filename.csv", header=TRUE, sep=",", format =
"%H:%M:%S %m.%d.%Y", tz = "", strip.white = TRUE)

and the R always shows up the message:

Error in read.zoo("filename.csv", header = TRUE, sep = ",", format =
"%H:%M:%S %m.%d.%Y",  :
  index contains NAs

At the beginning, I thought it is the problem of NA, and tried to
removed the records with NA. Still, the message appeared until I
reduce the number of records to around 280 and it works well with or
without NAs.

Does anyone has any idea to solve the problem?

Regards,
Keith

On Wed, Jul 15, 2009 at 5:08 PM, Gabor
Grothendieck wrote:
> Try the zoo package:
>
> Lines <- "time[sec] , Factor1 , Factor2
> 00:00:00 01.01.2007 , 0. , 0.176083
> 01:00:00 01.01.2007 , 0. , 0.176417
> 11:00:00 10.06.2007 , 0. , 0.148250
> 12:00:00 10.06.2007 , NA , 0.147000
> 13:00:00 10.06.2007 , NA , 0.144417"
>
> library(zoo)
> library(chron)
> z <- read.zoo(textConnection(Lines), sep = ",", header = TRUE,
>        format = "%H:%M:%S %m.%d.%Y", tz = "", strip.white = TRUE)
> plot(z)
>
> and read the three vignetttes (pdf documents) that come with it.
>
>
> On Wed, Jul 15, 2009 at 10:07 AM, Keith wrote:
>> Hello everyone,
>>
>> I am just a tyro in R and would like your kindly help for some
>> problems which I've been struggling for a while but still in vain.
>>
>> I have a time-series file (with some missing value ) which looks like
>>
>> time[sec] , Factor1 , Factor2
>> 00:00:00 01.01.2007 , 0. , 0.176083
>> 01:00:00 01.01.2007 , 0. , 0.176417
>> [ ... ]
>> 11:00:00 10.06.2007 , 0. , 0.148250
>> 12:00:00 10.06.2007 , NA , 0.147000
>> 13:00:00 10.06.2007 , NA , 0.144417
>> [ ... ]
>>
>> and I would like to do some basic time-series analyses using R. The
>> first idea is to plot these time-series events and the main problem
>> was the handling of the date/time format in the 1st column. I was
>> using the script below to deal with:
>>
>> data <- 
>> read.table("file",header=TRUE,sep=",",colClasses=c("character","numeric","numeric"))
>> data$time.sec. <- as.POSIXct(data$time.sec.,format="%H:%M:%S %d.%m.%Y")
>> dataTs <- as.ts(data)
>> plot.ts(dataTs)
>>
>> Then, the plot showed up with 3 subplots in one plot. The 1st is the
>> linear line with the x-axis being just the sequence of orders and
>> y-axis being wrong numbers which is completely wrong. The 2nd and the
>> 3rd are correct but the x-axis is still wrong. Does anyone know how to
>> plot correct Factor1 and Factor2 with respect to the 1st column time
>> format? Or, probably should I use some other packages? Besides, how
>> can I plot these two time-series data (Factor1 and Factor2) in two
>> separate plots?
>>
>> Best regards,
>> Keith
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>

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

Re: [R] Plotting hourly time-series data loaded from file using plot.ts

[Keith]

Yep, it's my mistake while manipulating the data. Now it works fine
except there is a warning message coming out:

Warning message:
In zoo(rval, ix) :
  some methods for “zoo” objects do not work if the index entries in
‘order.by’ are not unique

I checked the faq and some other documents, it seems this warning
implies that I have duplicated data. However, I checked the data and
it's hourly recorded throughout the year 2007, no duplication found.

Besides, while I plot the dataset, it looked what as I expected but
also showed the warning messages:

1: In zoo(rval, x.index[i]) :
  some methods for “zoo” objects do not work if the index entries in
‘order.by’ are not unique
2: In zoo(rval, x.index[i]) :
  some methods for “zoo” objects do not work if the index entries in
‘order.by’ are not unique
3: In zoo(rval, x.index[i]) :
  some methods for “zoo” objects do not work if the index entries in
‘order.by’ are not unique
4: In zoo(rval, x.index[i]) :
  some methods for “zoo” objects do not work if the index entries in
‘order.by’ are not unique

These looks the same problem again. Does anyone has any idea solving this?

Regards,
Keith

On Thu, Jul 16, 2009 at 4:51 PM, Gabor
Grothendieck wrote:
> By the way, note that read.zoo passes the ... arguments to read.table
> and so can use the same skip= and nrows= arguments that read.table
> uses.  These can be used to read in a subset of the rows.
>
> On Thu, Jul 16, 2009 at 10:35 AM, Gabor
> Grothendieck wrote:
>> There is no such limitation.   There is likely a data problem with
>> one or more records past the 280th one.
>>
>> Suggest you remove the first 280 and then divide the remaining
>> in half and try each half and keep dividing that way until you have
>> located the offending record or records.
>>
>> On Thu, Jul 16, 2009 at 8:51 AM, Keith wrote:
>>> Thanks Gabor,
>>>
>>> I tried a little bit, and your example works. However, it seems that
>>> the read.zoo has a limitation of records up to around 300 !? I took
>>> your suggestion and modified a little bit in order to read from the
>>> file which contains about 9000 records:
>>>
>>> dataTs <- read.zoo("filename.csv", header=TRUE, sep=",", format =
>>> "%H:%M:%S %m.%d.%Y", tz = "", strip.white = TRUE)
>>>
>>> and the R always shows up the message:
>>>
>>> Error in read.zoo("filename.csv", header = TRUE, sep = ",", format =
>>> "%H:%M:%S %m.%d.%Y",  :
>>>  index contains NAs
>>>
>>> At the beginning, I thought it is the problem of NA, and tried to
>>> removed the records with NA. Still, the message appeared until I
>>> reduce the number of records to around 280 and it works well with or
>>> without NAs.
>>>
>>> Does anyone has any idea to solve the problem?
>>>
>>> Regards,
>>> Keith
>>>
>>> On Wed, Jul 15, 2009 at 5:08 PM, Gabor
>>> Grothendieck wrote:
>>>> Try the zoo package:
>>>>
>>>> Lines <- "time[sec] , Factor1 , Factor2
>>>> 00:00:00 01.01.2007 , 0. , 0.176083
>>>> 01:00:00 01.01.2007 , 0. , 0.176417
>>>> 11:00:00 10.06.2007 , 0. , 0.148250
>>>> 12:00:00 10.06.2007 , NA , 0.147000
>>>> 13:00:00 10.06.2007 , NA , 0.144417"
>>>>
>>>> library(zoo)
>>>> library(chron)
>>>> z <- read.zoo(textConnection(Lines), sep = ",", header = TRUE,
>>>>        format = "%H:%M:%S %m.%d.%Y", tz = "", strip.white = TRUE)
>>>> plot(z)
>>>>
>>>> and read the three vignetttes (pdf documents) that come with it.
>>>>
>>>>
>>>> On Wed, Jul 15, 2009 at 10:07 AM, Keith wrote:
>>>>> Hello everyone,
>>>>>
>>>>> I am just a tyro in R and would like your kindly help for some
>>>>> problems which I've been struggling for a while but still in vain.
>>>>>
>>>>> I have a time-series file (with some missing value ) which looks like
>>>>>
>>>>> time[sec] , Factor1 , Factor2
>>>>> 00:00:00 01.01.2007 , 0. , 0.176083
>>>>> 01:00:00 01.01.2007 , 0. , 0.176417
>>>>> [ ... ]
>>>>> 11:00:00 10.06.2007 , 0. , 0.148250
>>>>> 12:00:00 10.06.2007 , NA , 0.147000
>>>>> 13:00:00 10.06.2007 , NA , 0.144417
>>>>> [ ... ]
>>>>>
>>>>> and I would like to do some basic time-series analyses using R. The
>>>>&

[R] Lag representation in ccf() while zoo object is used?

[Keith]

Dear All,

I have 2 time-series data sets and would like to check the cross
correlation. These data sets were set as a zoo object, called data,
and in general look like:

V1 V2
2007-01-01 00:00:00  0.0   0.176083
2007-01-01 01:00:00  0.0   0.176417
2007-01-01 02:00:00  0.0   0.175833
2007-01-01 03:00:00  0.0   0.175833
2007-01-01 04:00:00  0.3   0.176000
2007-01-01 05:00:00  1.8   0.176250
2007-01-01 06:00:00  2.0   0.177583
2007-01-01 07:00:00  0.2   0.178333
2007-01-01 08:00:00  0.3   0.178167
2007-01-01 09:00:00  3.2   0.178417

When I applied the ccf method, ccf(data$V1, data$V2), I noticed the
lag is every 3600 which is a little surprising to me. I was thinking
the lag should be 1, but it seems the lag unit becomes 3600. I guess
the number 3600 representing 3600 "seconds" because of the zoo object.
I am not sure if I'm right and would like someone here could certify
this (or not). Besides, does anyone know any default argument to
adjust the 3600 into 1 while plotting? The only idea I have is to
divide the lag manually by 3600 and then plot it later.

Regards,
Keith

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

[Keith]

Dear R-users,

I tried to read recorded data into R, but confronted the problem about
time zone. The data was recorded in GMT+1 (without summer time
changing), and at first I tried to do in the way:

data <- read.zoo("data.txt", header=FALSE, sep=",", format="%H:%M:%S
%d.%m.%Y", strip.white=TRUE, tz="GMT+1", skip=3)

However, it failed. Then, I found the description of "Time Zones" in R
help and noticed the Olson database is stored in the specific folder
under Windows OS which is the OS I'm working on. I noticed there is no
file related to GMT+1 time zone if I understood correctly. So far my
workaround is to set the tz to Africa/Lagos where is the city in
Africa sharing the same time zone. However, it's a little bit weird
because the data was not collected in that city nor the place near by.

My question is that if there is other way to set the tz to GMT+1
instead using my workaround to have a neutral meaning in my script?

Regards,
Keith

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

Re: [R] Question

[Keith Jewell]

On 15/04/2018 17:26, Marc Girondot via R-help wrote:

Le 15/04/2018 à 17:56, alireza daneshvar a écrit :

break-down point


Can you explain more what you plan to do and give an example of what you 
have tried to do until now to do a "break down point" in R. Perhaps a 
"break down point" is common in your field, but I have no idea about 
what it is !


https://en.wikipedia.org/wiki/Breakdown_(1997_film)

Sincerely

Marc

Perhaps the OP means "breakpoint" in which case the 'strucchange' 
package might be relevant 



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

[R] Problem with choose.files(default=..., multi=FALSE)

[Keith Jewell]

I'm very hesitant to suggest that there's a bug in such a venerable R 
function, but I can't see what I'm doing wrong. Any comments are welcome


When using choose.files() where:
default = something
multi = FALSE
selected file path is shorter than the default
... then the returned value is at least as long as the default, 
characters from default appearing (wrongly) at the end of the returned 
value.


Example, in which all but the first choose.files() select 
"M:\\test\\target.dat". Note the last result.


> pathlong <- choose.files(caption = "long")
> pathlong # long file name to use as default for short selection
[1] 
"M:\\test\\Averyveryveryveryverylongfoldername\\Averyveryveryveryverylongfoldername\\Averyveryveryveryverylongfoldername\\target.dat"

> choose.files(caption = "short")  # no default without multi works
[1] "M:\\test\\target.dat"
> choose.files(default=pathlong, caption = "short") # default without 
multi= works

[1] "M:\\test\\target.dat"
> choose.files(caption = "short", multi = FALSE) # multi = FALSE 
without default works

[1] "M:\\test\\target.dat"
> choose.files(default=pathlong, caption = "short", multi = TRUE) # 
multi = TRUE with default works

[1] "M:\\test\\target.dat"
> choose.files(default=pathlong, caption = "short", multi = FALSE) # 
multi = FALSE with default fails
[1] 
"M:\\test\\target.dat\\ryveryverylongfoldername\\Averyveryveryveryverylongfoldername\\Averyveryveryveryverylongfoldername\\target.dat"


> # in case it's relevant
> sessionInfo()
R version 3.4.0 (2017-04-21)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows Server 2008 R2 x64 (build 7601) Service Pack 1

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United 
Kingdom.1252
[3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C 


[5] LC_TIME=English_United Kingdom.1252

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

[9] base

other attached packages:
 [1] CBRIutils_1.0   stringr_1.2.0   svSocket_0.9-57 TinnR_1.0-5 
R2HTML_2.3.2
 [6] Hmisc_4.0-3 ggplot2_2.2.1   Formula_1.2-1   survival_2.41-3 
lattice_0.20-35


loaded via a namespace (and not attached):
 [1] RColorBrewer_1.1-2  htmlTable_1.9   digest_0.6.12 
htmltools_0.3.6
 [5] splines_3.4.0   scales_0.4.1grid_3.4.0 
checkmate_1.8.2
 [9] devtools_1.12.0 knitr_1.15.1munsell_0.4.3 
compiler_3.4.0
[13] tibble_1.3.0nnet_7.3-12 acepack_1.4.1 
Matrix_1.2-10
[17] svMisc_0.9-70   plyr_1.8.4  base64enc_0.1-3 
data.table_1.10.4
[21] stringi_1.1.5   magrittr_1.5gtable_0.2.0 
colorspace_1.3-2
[25] foreign_0.8-68  cluster_2.0.6   gridExtra_2.2.1 
htmlwidgets_0.8
[29] withr_1.0.2 lazyeval_0.2.0  backports_1.0.5 
memoise_1.1.0

[33] rpart_4.1-11Rcpp_0.12.10latticeExtra_0.6-28
>

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

Re: [R] Problem with choose.files(default=..., multi=FALSE)

[Keith Jewell]

Thanks for confirming that I wasn't being stupid :-}

When using default=pathlong I get the _correct_ starting directory...
(M:\test\Averyveryveryveryverylongfoldername\Averyveryveryveryverylongfoldername\Averyveryveryveryverylongfoldername) 

... both in the environment I indicated originally (Windows Server 2008 
R2 x64) and also in Windows 10 x64


Keith Jewell

On 09/05/2017 17:49, Duncan Murdoch wrote:

On 09/05/2017 12:06 PM, Keith Jewell wrote:

I'm very hesitant to suggest that there's a bug in such a venerable R
function, but I can't see what I'm doing wrong. Any comments are welcome


Yes, it looks like a bug.  One other thing I find a little strange: the
starting directory seems wrong when I have the pathlong default.  Did
you see that?  (I'm in Windows 10, not the same version as you.)

Duncan Murdoch



When using choose.files() where:
 default = something
 multi = FALSE
 selected file path is shorter than the default
... then the returned value is at least as long as the default,
characters from default appearing (wrongly) at the end of the returned
value.

Example, in which all but the first choose.files() select
"M:\\test\\target.dat". Note the last result.

 > pathlong <- choose.files(caption = "long")
 > pathlong # long file name to use as default for short selection
[1]
"M:\\test\\Averyveryveryveryverylongfoldername\\Averyveryveryveryverylongfoldername\\Averyveryveryveryverylongfoldername\\target.dat"

 > choose.files(caption = "short")  # no default without multi works
[1] "M:\\test\\target.dat"
 > choose.files(default=pathlong, caption = "short") # default without
multi= works
[1] "M:\\test\\target.dat"
 > choose.files(caption = "short", multi = FALSE) # multi = FALSE
without default works
[1] "M:\\test\\target.dat"
 > choose.files(default=pathlong, caption = "short", multi = TRUE) #
multi = TRUE with default works
[1] "M:\\test\\target.dat"
 > choose.files(default=pathlong, caption = "short", multi = FALSE) #
multi = FALSE with default fails
[1]
"M:\\test\\target.dat\\ryveryverylongfoldername\\Averyveryveryveryverylongfoldername\\Averyveryveryveryverylongfoldername\\target.dat"


 > # in case it's relevant
 > sessionInfo()
R version 3.4.0 (2017-04-21)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows Server 2008 R2 x64 (build 7601) Service Pack 1

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United
Kingdom.1252
[3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C

[5] LC_TIME=English_United Kingdom.1252

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

other attached packages:
  [1] CBRIutils_1.0   stringr_1.2.0   svSocket_0.9-57 TinnR_1.0-5
R2HTML_2.3.2
  [6] Hmisc_4.0-3 ggplot2_2.2.1   Formula_1.2-1   survival_2.41-3
lattice_0.20-35

loaded via a namespace (and not attached):
  [1] RColorBrewer_1.1-2  htmlTable_1.9   digest_0.6.12
htmltools_0.3.6
  [5] splines_3.4.0   scales_0.4.1grid_3.4.0
checkmate_1.8.2
  [9] devtools_1.12.0 knitr_1.15.1munsell_0.4.3
compiler_3.4.0
[13] tibble_1.3.0nnet_7.3-12 acepack_1.4.1
Matrix_1.2-10
[17] svMisc_0.9-70   plyr_1.8.4  base64enc_0.1-3
data.table_1.10.4
[21] stringi_1.1.5   magrittr_1.5gtable_0.2.0
colorspace_1.3-2
[25] foreign_0.8-68  cluster_2.0.6   gridExtra_2.2.1
htmlwidgets_0.8
[29] withr_1.0.2 lazyeval_0.2.0  backports_1.0.5
memoise_1.1.0
[33] rpart_4.1-11Rcpp_0.12.10latticeExtra_0.6-28
 >

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] expand.model.frame issue

[Keith Jewell]

Following on from a stackoverflow question "Why does this simple 
function calling `lm(..., subset)` fail?"


--
myfun <- function(form., data., subs.) lm(form., data., subs.)
myfun(mpg ~ cyl + hp, mtcars, TRUE)
## Error in eval(expr, envir, enclos) : object 'subs.' not found
-

The answer to the stated question was in ?lm "If not found in data, the 
variables are taken from environment(formula), typically the environment 
from which lm is called"; the environment of the formula (mpg ~ cyl + 
hp) does not contain 'subs.'. A fix is quite straightforward, set the 
environment of the formula to that of the function, which does contain 
'subs.'. There are multiple ways of doing that, this works but to me 
seems a bit "clunky":

---
myfun <- function(form., data., subs.) lm(as.formula(deparse(form.)), 
data., subs.)

myfun(mpg ~ cyl + hp, mtcars, TRUE)
--
To me this seems more elegant, but then I have no taste :-}
--
myfun <- function(form., data., subs.){
  environment(form.) <- environment()
  lm(form., data., subs.)}
myfun(mpg ~ cyl + hp, mtcars, TRUE)
--

But the OP went on to consider `expand.model.frame` e.g.
-
myfun <- function(form., data., subs.){
  environment(form.) <- environment()
  model <- lm(form., data., subs.)
  print(ls(envir = environment(formula(model
  expand.model.frame(model, ~drat)}
myfun(mpg ~ cyl + hp, mtcars, TRUE)
## [1] "data." "form." "model" "subs."
## Error in eval(expr, envir, enclos) : object 'subs.' not found
-

myfun can be fixed by (e.g.) avoiding the subset argument of lm

myfun <- function(form., data., subs.){
  environment(form.) <- environment()
  model <- lm(form., data.[subs.,])
  expand.model.frame(model, ~drat)}
myfun(mpg ~ cyl + hp, mtcars, TRUE)

... but this message is about the apparent inconsistency between the 
behaviour of expand.model.frame and the help text which says:

?expand.model.frame:
---
Usage

expand.model.frame(model, extras,
   envir = environment(formula(model)),
   na.expand = FALSE)

envir   an environment to evaluate things in
-

In the example of the `expand.model.frame` issue above the result of the 
'ls()' clearly shows that 'subs.' is in that environment, but 
expand.model.frame fails to find it.


Am I misunderstanding?
Or is there an error in the help text?
Or is there a bug in expand.model.frame?

=
I don't think this is relevant, but for completeness
> sessionInfo()
R version 3.3.0 (2016-05-03)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows Server 2008 R2 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United 
Kingdom.1252
[3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C 


[5] LC_TIME=English_United Kingdom.1252

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


other attached packages:
 [1] CBRIutils_1.0   stringr_1.0.0   svSocket_0.9-57 TinnR_1.0-5 
R2HTML_2.3.1Hmisc_3.17-4ggplot2_2.1.0

 [8] Formula_1.2-1   survival_2.39-4 lattice_0.20-33

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.5 magrittr_1.5cluster_2.0.4 
splines_3.3.0   devtools_1.11.1
 [6] munsell_0.4.3   colorspace_1.2-6plyr_1.8.3 
nnet_7.3-12 grid_3.3.0
[11] data.table_1.9.6gtable_0.2.0latticeExtra_0.6-28 
withr_1.0.1 svMisc_0.9-70
[16] digest_0.6.9Matrix_1.2-6gridExtra_2.2.1 
RColorBrewer_1.1-2  acepack_1.3-3.3
[21] rpart_4.1-10memoise_1.0.0   stringi_1.1.1 
scales_0.4.0foreign_0.8-66

[26] chron_2.3-47

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

Re: [R] [FORGED] Re: Generate list if sequence form two vector element

[Keith Jewell]

or as a one-liner
mapply(pmin(a, b), pmax(a,b), FUN=seq, SIMPLIFY=FALSE)

On 22/06/2016 10:23, peter dalgaard wrote:

There's also

mapply(a, b, FUN=seq, SIMPLIFY=FALSE)

(turn off simplication so that you don't unexpectedly get a matrix whenever all 
elements of results have same length. This also affects apply()-based 
solutions.)

...except that according to original spec, one should ensure a < b. So

myseq <- function(a,b) if(a
On 22 Jun 2016, at 10:42 , Jim Lemon  wrote:

Now why didn't I think of that?

apply(matrix(c(a,b),ncol=2),1,function(x)x[1]:x[2])

Jim

On Wed, Jun 22, 2016 at 6:14 PM, Rolf Turner  wrote:

On 22/06/16 20:00, Jim Lemon wrote:


Hi Tanvir,
Not at all elegant, but:

make.seq<-function(x) return(seq(x[1],x[2]))
apply(matrix(c(a,b),ncol=2),1,make.seq)



Not sure that this is more "elegant" but it's a one-liner:

lapply(1:length(a),function(i,a,b){a[i]:b[i]},a=a,b=b)

cheers,

Rolf



On Wed, Jun 22, 2016 at 5:32 PM, Mohammad Tanvir Ahamed via R-help
 wrote:


Hi,
I want to do the follow thing

Input :
a <- c(1,3,6,9)


b<-c(10,7,20,2)


Expected outcome :

d<-list(1:10,3:7,6:20,2:9)


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

[R] [R-pkgs] New package: simstudy

[Goldfeld, Keith]

Greetings –

A new package “simstudy” is now available on CRAN. What started as a small 
number of functions that enabled me to quickly generate simple data sets for 
teaching and power/sample size calculations has grown into a more robust set of 
tools that allows users to simulate more complex data sets in order to explore 
modeling techniques or better understand data generating processes. The user 
specifies a set of relationships between covariates in table form (the table 
can be built interactively or created externally as a csv file), and generates 
data based on these specifications. The final data sets can represent data from 
randomized control trials, observed (non-randomized) studies, repeated measure 
(longitudinal) designs, and cluster randomized trials. Missingness can be 
generated using various mechanisms (MCAR, MAR, NMAR). Currently, data can be 
generated from normal/Gaussian, binary, Poisson, truncated Poisson, Gamma, and 
uniform distributions. Survival data can also be generated.

I will be adding functionality over time, and will be particularly interested 
in knowing what userRs would be interested in having me add. I look forward to 
hearing your comments.

- Keith


Keith Goldfeld
Department of Population Health
School of Medicine, New York University
227 East 30th Street, 6th Floor
New York, NY  10016



This email message, including any attachments, is for the sole use of the 
intended recipient(s) and may contain information that is proprietary, 
confidential, and exempt from disclosure under applicable law. Any unauthorized 
review, use, disclosure, or distribution is prohibited. If you have received 
this email in error please notify the sender by return email and delete the 
original message. Please note, the recipient should check this email and any 
attachments for the presence of viruses. The organization accepts no liability 
for any damage caused by any virus transmitted by this email.
=

[[alternative HTML version deleted]]

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

Re: [R] Strange message after reading multiple scripts from one folder

[Keith Jewell]

I can't immediately see it in the help text but it seems that source 
returns a list with two named elements; value and visible.


I surmise that it is returned using withVisible (qv).

KJ

On 29/07/2016 13:26, jim holtman wrote:

Hard to tell without seeing the scripts.  Do you have a matrix in your
scripts that have "value" and "visible" as row names?  You probably have
some statement that is causing output and so the problem is "your" as to
how to avoid the message.  So look at your scripts to see if anything
refers to either "value" or "visible", and then you might find the cause of
your problem.


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.

On Fri, Jul 29, 2016 at 6:52 AM, Frank S.  wrote:


Dear list,

I have one folder named "scripts_JMbayes", wich contains 10 R scripts.
I can read them properly by doing:


pathnames <- list.files(pattern="[.]R", path="Mydir/scripts_JMbayes",

full.names = TRUE)

sapply(pathnames, USE.NAMES = FALSE, FUN = source,)


However, R generates the following message:

 [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9]  [,10]
value   ? ? ? ? ? ? ? ? ? ?
visible FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

What does it mean and what should I change to avoid this message?
Any help would be appreciated!

Best,

Frank


 [[alternative HTML version deleted]]

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



[[alternative HTML version deleted]]



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

Re: [R] Problem with particular file in XML package?

[Keith Jewell]

The OP asked "Has anyone else had trouble with the XML package lately 
and if so, how did you resolve it?"


For what it's worth...

I failed to install XML using install() with the defaults; I can't 
remember the exact error message, something about access denied.


Downloading XML_3.98-1.1.zip with Internet Explorer 
 
gave a more informative error:


Threat Source: http://cran.r-project.or...trib/3.2/XML_3.98-1.1.zip

The file requested could not be scanned by Sophos Anti-Virus. This means 
it could be encrypted or may contain errors that prevent full scanning. 
As a result, the file was blocked from downloading.

=
... so my access was blocked by our systems anti-virus.

Our IT people bypassed the anti-virus to download the zip from which I 
successfully installed.


I note that in the library as installed there is a file
   ...\XML\exampleData\dtd.zip
Double-clicking in Windows Explorer gives an error:
=
Windows cannot open the folder. The Compressed (zipped) Folder ... is 
invalid.

=

I speculate that Sophos Anti-Virus could not scan dtd.zip because it 
tried to open it as a zipped folder and failed. I don't know if it 
really is a zipped folder or if that's just its name :-O


This 2013 thread seems relevant 
https://stat.ethz.ch/pipermail/r-sig-mac/2013-November/010494.html



On 28/05/2015 07:38, Prof Brian Ripley wrote:

This really should have been sent to the package maintainer.  But that
the zip file is corrupt has been reported several times, and does not
block installation for anyone else, so your (plural) diagnosis is wrong.

On 28/05/2015 03:56, Gen wrote:

I have been attempting to install the R devtools package at work.  The
version of R is 3.1.2 (Pumpkin Helmet).  However, the installation of
devtools fails because devtools depends on rversions which in turn
depends
upon the XML package (XML_3.98-1.1.tar.gz), and the XML package is not
importing correctly for us.

One of our system administrators tried scanning through the files in the
XML package, and he said that the particular file:
/src/contrib/XML_3.98-1.1.tar.gz/XML/inst/exampleData/dtd.zip looks
corrupted.  The actual error message he received was: "Archive parsing
failed!  (Data is corrupted)."  For the record, I tried downloading an
older version of the XML package (XML_3.95-0.1.tar.gz) but that was also
without success -- this time there was a separate error message about not
being able to locate xml2-config.  (Perhaps XML_3.95-0.1.tar.gz is
just not
compatible with R version 3.1.2?)

I tried browsing over to the "CRAN checks" link for the XML package and
noticed several red warning messages under the "Status" column -- not
sure
if that is typical?  Has anyone else had trouble with the XML package
lately and if so, how did you resolve it?  Would it be possible to remove
the potentially corrupted file and then re-upload the package source
XML_3.98-1.1.tar.gz to the CRAN webpage?  Thanks for your
help/suggestions!

[[alternative HTML version deleted]]

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


PLEASE do, including what it says about HTML mail, 'at a minimum'
information required and upgrading before posting: R 3.1.2 is already 2
versions obsolete.



__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] predict.poly for multivariate data

[Keith Jewell]

A recent stackoverflow post 
<http://stackoverflow.com/questions/31134985> "How do you make R poly() 
evaluate (or “predict”) multivariate new data (orthogonal or raw)?" made 
me look at poly and polym again.


predict.poly doesn't work with multivariate data because for such data 
poly calls polym which does not:

a) return an object inheriting from class poly;
b) return the coefficients needed to make predictions;
c) accept coefficients as an argument or include code to make predictions.

This does lead to some wrong answers without warnings. e.g.
## vanilla poly and polym ###
library(datasets)
alm <- lm(stack.loss ~ poly(Air.Flow, Water.Temp, degree=3), stackloss)
# "correct" prediction values [1:10]
alm$fitted.values[1:10]
# predict(alldata)[1:10] gives correct values
predict(alm, stackloss)[1:10]
# predict(partdata) gives wrong values
predict(alm, stackloss[1:10,])
#
I guess - but haven't confirmed - that with multivariate newdata 
predict.lm(alm, newdata) calculates new orthogonal polynomials based on 
newdata rather than applying the original coefficients.


Below I append versions of:
a) polym edited to address the three points above;
b) poly slightly edited to reflect the changes in polym;
c) predict.poly unaltered, just to get it in the same environment as 
polym and poly for testing.

After implementing these the three sets of predictions above all agree.

I'm very ready to believe that I've got the wrong end of the stick 
and/or my suggestions can be improved so I welcome correction.


Otherwise, how do I go about getting these changes implemented?
I see stats is maintained by R Core Team. Are they likely to pick it up 
from here, or do I need to take any other action?


Best regards

Keith Jewell

### polym ##
polym <- function (..., degree = 1, coefs = NULL, raw = FALSE)
# add coefs argument
{
  if(is.null(coefs)) {
dots <- list(...)
nd <- length(dots)
if (nd == 0)
  stop("must supply one or more vectors")
if (nd == 1)
  return(poly(dots[[1L]], degree, raw = raw))
n <- sapply(dots, length)
if (any(n != n[1L]))
  stop("arguments must have the same length")
z <- do.call("expand.grid", rep.int(list(0:degree), nd))
s <- rowSums(z)
ind <- (s > 0) & (s <= degree)
z <- z[ind, ]
s <- s[ind]
aPoly <- poly(dots[[1L]], degree, raw = raw) # avoid 2 calcs
res <- cbind(1, aPoly)[, 1 + z[, 1]]
# attribute "coefs" = list of coefs from individual variables
if (!raw) coefs <- list(attr(aPoly, "coefs"))
for (i in 2:nd) {
  aPoly <- poly(dots[[i]], degree, raw = raw)
  res <- res * cbind(1, aPoly)[, 1 + z[, i]]
  if (!raw) coefs <- c(coefs, list(attr(aPoly, "coefs")))
}
colnames(res) <- apply(z, 1L, function(x) paste(x, collapse = "."))
attr(res, "degree") <- as.vector(s)
if (!raw) attr(res, "coefs") <- coefs
class(res) <- c("poly", "matrix") # add poly class
res
  }
  else
  {
nd <- length(coefs)# number of variables
newdata <- as.data.frame(list(...)) # new data
if (nd != ncol(newdata)) stop("wrong number of columns in newdata")
z <- do.call("expand.grid", rep.int(list(0:degree), nd))
s <- rowSums(z)
ind <- (s > 0) & (s <= degree)
z <- z[ind, ]
res <- cbind(1, poly(newdata[[1]], degree=degree, 
coefs=coefs[[1]]))[, 1 + z[, 1]]
for (i in 2:nd) res <- res*cbind(1, poly(newdata[[i]], 
degree=degree, coefs=coefs[[i]]))[, 1 + z[, i]]

colnames(res) <- apply(z, 1L, function(x) paste(x, collapse = "."))
res
  }
}
##

 poly ##
poly <- function (x, ..., degree = 1, coefs = NULL, raw = FALSE)
{
  dots <- list(...)
  if (nd <- length(dots)) {
if (nd == 1 && length(dots[[1L]]) == 1L)
  degree <- dots[[1L]]
# pass coefs argument as well
else return(polym(x, ..., degree = degree, coefs=coefs, raw = raw))
  }
  if (is.matrix(x)) {
m <- unclass(as.data.frame(cbind(x, ...)))
# pass coefs argument as well
return(do.call("polym", c(m, degree = degree, raw = raw, 
list(coefs=coefs

  }
  if (degree < 1)
stop("'degree' must be at least 1")
  if (anyNA(x))
stop("missing values are not allowed in 'poly'")
  n <- degree + 1
  if (raw) {
Z <- outer(x, 1L:degree, "^")
colnames(Z) <- 1L:degree
attr(Z, "degree") <- 1L:degree
class(Z) <- c("poly", "matrix")
return(Z)
  }
  if (is.null(coefs)) {
if (degree >= length(unique(x)))
  stop("'degree' must be less than number of unique 

Re: [R] rgl::writeWebGL( , prefix = , )

[Keith Jewell]

Thanks Duncan, your suggestions led me to a solution.

Perhaps this could be reflected in the help, but I'll leave that 
decision to you.


It comes down to the template. As well as including a single line for 
each scene containing

  paste("%", prefix, "WebGL%"); e.g. %WebGL% or %AWebGL%
the  tag must contain an onload attribute with an element for 
each scene

  paste0(prefix, "webGLStart();")
e.g. 

While I'm suggesting additions to the help, it took me a little while to 
work out that when writing multiple scenes the result from one 
writeWebGL was the template for the next. E.g. in the above example with 
two scenes:


outfile <- writeWebGL(dir=getwd(), template = file.path(getwd(), 
"template.html"), prefix="")

# position to next scene
outfile <- writeWebGL(dir=getwd(), template = outfile, prefix="A")

Thanks for your help (and a really nice package!)

Keith J

On 03/02/2015 15:14, Duncan Murdoch wrote:

On 03/02/2015 9:43 AM, keith.jew...@campdenbri.co.uk wrote:

Dear all,

I am using writeWebGL to create an HTML page containing an interactive
3D plot. It works fine with the default prefix="" but fails when I
specify a prefix "for different scenes displayed on the same web page"
(quoting ?writeWebGL).  I'm sure I'm misreading the help, and would
appreciate guidance.

Briefly, it works fine with the default writeWebGL( ,prefix="", ) and
the template containing %WebGL%
I have not been able to make it work with any other value of prefix;
e.g. writeWebGL( ,prefix="A",) and the template containing %AWebGL%

Here is code illustrating the problem.

First create three templates:
a) Vanilla: copied system.file(file.path("WebGL", "template.html"),
package="rgl") to file.path(getwd(), "template.html")

b) First attempt: ?writeWebGL says # "[the template] should contain a
single line containing paste("%", prefix, "WebGL%"), e.g. %WebGL% with
the default empty prefix"
paste("%", "A", "WebGL%")
# [1] "% A WebGL%"
so file.path(getwd(), "templateA.html") is a copy of (a) replacing
%WebGL% with % A WebGL%

c) Second attempt: file.path(getwd(), "templateB.html") is  a copy of
(a) replacing %WebGL% with %AWebGL%

then, in R
#---
library(rgl)
plot3d(1:5, 1:5, 1:5) # generate rgl scene
#---
# a) vanilla
writeWebGL(dir=getwd(), template = file.path(getwd(),
"template.html"), prefix="")
# works OK; result opens and works in IE
#
# b) First attempt, my reading of ?writeWebGL
writeWebGL(dir=getwd(), template = file.path(getwd(),
"templateA.html"), prefix="A")
# Error in writeWebGL(dir = getwd(), template = file.path(getwd(),
"templateA.html"),  :
#   template ‘m://templateA.html’ does not contain
%AWebGL%
# so it looks as if the help is trivially wrong, it should be paste0
paste0("%", "A", "WebGL%")


Yes, that's right.  I'll fix it.


# [1] "%AWebGL%"
#
# c) second attempt using %AWebGL%
writeWebGL(dir=getwd(), template = file.path(getwd(),
"templateB.html"), prefix="A")
# runs without error in R but IE displays "You must enable Javascript
to view this page properly."
#--

I don't understand why (c) is different from (a).


There may be an error in the generated Javascript.  In Firefox, you
could ask to see the browser console log, and it would report if there
was an error on the page; sometimes those make the Javascript fail, and
it falls back to the error message you saw.  I don't know how/if you can
do that in IE.




Here are the system details:

R version 3.1.0 (2014-04-10)
Platform: i386-w64-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=English_United Kingdom.1252
[2] LC_CTYPE=English_United Kingdom.1252
[3] LC_MONETARY=English_United Kingdom.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United Kingdom.1252

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

other attached packages:
[1] knitr_1.8   animation_2.3   rgl_0.95.1158   CBRIutils_1.0


That's an old version of rgl; current on CRAN is 0.95.1201.
<http://cran.r-project.org/src/contrib/rgl_0.95.1201.tar.gz> (CRAN OSX
currently has an old binary; I don't recommend that you use it.  I don't
know why they haven't updated to the current one.) r-forge has even
newer versions, but I'm in the middle of some changes there, so I don't
recommend using that version right now.

Duncan Murdoch
<http://cran.r-project.org/src/contrib/rgl_0.95.1201.tar.gz>

  [5] stringr_0.6.2   svSocket_0.9-55 TinnR_1.0-5 R2HTML_2.3.1
  [9] Hmisc_3.12-2Formula_1.1-1   survival_2.37-7

loaded via 

Re: [R] Windows R doesn't recognize shortcuts ?

[Keith Jewell]

On 23/07/2014 14:21, Duncan Murdoch wrote:

On 23/07/2014 9:08 AM, ce wrote:

Hi All,

In Windows 7 , R installation:

R version 3.1.1 Patched (2014-07-14 r66149) -- "Sock it to Me"
Copyright (C) 2014 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

it doesn't recognize shortcuts in path :

>  list.files(path = "cygwin")
character(0)

cygwin is a shortcut,  in properties window  Target shows :
C:\Users\me\cygwin64\home\me
Real path works :

list.files(path = "C:/Users/me/cygwin64/home/me")
   [1] "1010week.sh""10week.sh"
"a.R""aa.sh"



I don't think R should recognize that.  Windows wouldn't recognize it
either, if you used "dir cygwin" in a shell, for example.

Duncan Murdoch


The "shortcut" which appears as cygwin is actually a file
called cygwin.lnk

readWindowsShellLink {R.utils} reads such files:
> readWindowsShellLink(con="cygwin.lnk")

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 check to see if a variable is within a range of another variable

[Keith Jewell]

On 01/10/2014 23:54, Peter Alspach wrote:

Tena koe Kate

If kateDF is a data.frame with your data, then

apply(kateDF, 1, function(x) isTRUE(all.equal(x[2], x[1], check.attributes = 
FALSE, tolerance=0.1)))

comes close to (what I think) you want (but not to what you have illustrated in 
your 'eventual outcome').  Anyhow, it may be enough to allow you to get there.

HTH 

Peter Alspach



I seem to need to specify all.equal(..., scale) to get correct values 
for some data/tolerances:

--
aDF <- data.frame(a = 1:10, b=10:1)

apply(aDF, 1, function(x)
  isTRUE(all.equal(x[2], x[1], check.attributes = FALSE, tolerance = 5, 
scale = 1))

  )
#  [1] FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE
apply(aDF, 1, function(x)
  isTRUE(all.equal(x[2], x[1], check.attributes = FALSE, tolerance = 5))
)
#  [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE
--
I'm probably being stupid, but from reading ?all.equal I would have 
expected scale = 1 and the default scale = NULL to give identical 
results for the length one numerics being passed to all.equal.


Can anyone explain?

KJ

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

Re: [R] Is xyz point inside 3d convex hull?

[Keith Jewell]

Back in 2009 I posted some code to this list, see:
<http://tolstoy.newcastle.edu.au/R/e8/help/09/12/8784.html>

I submitted the function 'inhull' to the geometry package maintainer, 
but I don't think it was ever included.


HTH

Keith J
On 12/10/2014 21:24, Don McKenzie wrote:

Check the R-news archive with approrpriate keywords.  There was a long exchange 
awhile back when I asked a similar question.

On Oct 12, 2014, at 1:20 PM, Camilo Mora  wrote:


Hi everyone,

I wonder if there is a code in r that can generate a 3d convex hull from a 
data-frame containing 3 columns and then use another database with the same 
three columns and for each row determine if the xyz point is inside or not the 
convex hull generated with the first database?

The package geometry allows to calculate a hull and it's volume. I was planning 
to calculate the volume of the convex hull after adding each point in the 
second database and if the hull gets bigger then the point is out and if not 
then the point is in. A problem with this method is that I have over 10 million 
points and the calculation for each point will take a lot of time.

Any guidance will be greatly appreciated,

Thanks,

Camilo

[[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] Is it possible to define another kind of NA

[Keith Jewell]

On 13/11/2014 11:08, Marc Girondot wrote:

Le 13/11/2014 01:26, MacQueen, Don a écrit :

Along the lines of what Bert Gunter said, the ideal way to represent 
I agree that LDL is a special case of what could be named ODL (Out of
detection limit).
To answer to Bert Gunter, indeed if LDL (or ODL) values are changed into
NA, the results will be biased. That's why I would like to introduce
another category. I don't plan to just transform them as NA.

But thinking again about this problem, a LDL must be always associated
with one value (or two in the case of ODL) that indicates the detection
limit. In a dataset, all values have not necessarily the same limit
depending on the experimental conditions.
The best solution that I find is to use attributes to indicate the
limits. A NA attribute for a NA value will be treated as a "true" NA.
For exemple:

 > values <- c(NA, 29, 30, NA, 3)
 > attributes(values) <- list(ODL=c(NA, "[10, 40]", "[0, 40]", "[0,
40]", "[0, 40]"))
 > values
[1] NA 29 30 NA  3
attr(,"ODL")
[1] NA "[10, 40]" "[0, 40]"  "[0, 40]"  "[0, 40]"
 > values[3]
[1] 30
 > attributes(values)$ODL[3]
[1] "[0, 40]"
 > values[1]
[1] NA
 > attributes(values)$ODL[1]
[1] NA

The attributes are retained in data.frame. So it seems to be a good
solution.

 > essai <- data.frame(c1=values)
 > essai
   c1
1 NA
2 29
3 30
4 NA
5  3
 > essai$c1
[1] NA 29 30 NA  3
attr(,"ODL")
[1] NA "[10, 40]" "[0, 40]"  "[0, 40]"  "[0, 40]"

Thanks to the list members,

Marc

I strongly recommend you re-read and take action on Bert Gunter's 
comment, quoted here for truth!

-
Ouch!

The values are **NOT** missing -- they are (left) censored, and need
to be handled by appropriate censored data methods. I suggest you
(all!) either read up on this or consult someone locally who has
knowledge of such methods.

-- Bert

You are re-inventing the wheel and yours will probably end up square!
R already has facilities for handling censored data, e.g. Surv in the 
survival package (which despite its name is applicable to applications 
other than survival analysis).


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

Re: [R] Issue with results from 'summary' function in R

[Keith Jewell]

On 18/09/2015 13:08, Praveen Surendran wrote:

Hi all,

Attached table (that contains summary for a genetic association study) was read 
using the command:

test <- read.table('testDat.txt',header=FALSE,stringsAsFactors=FALSE)

Results from the summary of the attached table is provided below:


summary(test$V5)

Min. 1st Qu.  MedianMean 3rd Qu.Max.
   22070   22070   22070   22070   22070   22070

As we can see column 5 of this table contains only one value - 22072
I am confused as to why I am getting a value 22070 in the summary of this 
column.

I tested this using versions of R including - R version 3.2.1 (2015-06-18) -- 
"World-Famous Astronaut"

Thank you for looking at this issue.
Kind Regards,

Praveen.


> summary(22072, digits=5)
   Min. 1st Qu.  MedianMean 3rd Qu.Max.
  22072   22072   22072   22072   22072   22072

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

Re: [R] formula argument evaluation

[Keith Jewell]

On 12/04/2016 11:24, Adrian Dușa wrote:

I have a simple function such as:

foo <- function(x) {
 call <- lapply(match.call(), deparse)
 testit <- capture.output(tryCatch(eval(x), error = function(e) e))
 if (grepl("Error", testit)) {
 return(call$x)
 }
}

and I would like to detect a formula when x is not an object:

# this works

foo(A + B)

[1] "A + B"

# but this doesn't

foo(A + B => C)

Error: unexpected '=' in "foo(A + B ="

Can I prevent it from evaluating the "=" sign?
The addition sign "+" hasn't been evaluated, and I was hoping the "=" would
not get evaluated either. The "=>" sign is important for other purposes,
not related to this example.

Thank you in advance,
Adrian

--
Adrian Dusa
University of Bucharest
Romanian Social Data Archive
Soseaua Panduri nr.90
050663 Bucharest sector 5
Romania

[[alternative HTML version deleted]]


Did you mean
> foo (A + B >= C)
??

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

Re: [R] Antwort: Re: RStudio: Place for Storing Options

[Keith Jewell]

RStudio seems to not pay due regard to the distinction between APPDATA 
and LOCALAPPDATA.


See 
https://support.rstudio.com/hc/en-us/community/posts/200650543-File-history-and-project-history-stored-in-LOCALAPPDATA


On 23/02/2017 08:00, g.maub...@weinwolf.de wrote:

Hi Martin,

the command

%localappdata%\RStudio-Desktop

gives on my machine

"The command is written wrong or could not be found.".

I found "RStudio-Desktop" under

C:\Users\\AppData\Local\RStudio-Desktop

There references on created notebooks and presentations are stored in the
folder "RStudio-Desktop". RStudio config is not documented yet.

Kind regards

Georg




Von:Martin Maechler 
An: Jeff Newmiller ,
Kopie:  Martin Maechler ,
, R-help mailing list 
Datum:  23.02.2017 08:37
Betreff:Re: [R] RStudio: Place for Storing Options




Jeff Newmiller 
 on Sat, 11 Feb 2017 08:09:36 -0800 writes:


 > For the record, then, Google listened to my incantation of
 > "rstudio configuration file" and the second result was:

 >
https://support.rstudio.com/hc/en-us/articles/200534577-Resetting-RStudio-Desktop-s-State


 > RStudio Desktop is also open source, so you can download
 > the source code and look at the operating-system-specific
 > bits (for "where") if the above link goes out of date or
 > disappears.

Thanks a lot, Jeff!

And for the archives:  On reasonable OS's,  the hidden
directory/folder containing all the info is
   ~/.rstudio-desktop/
and if "things are broken" the recommendation is to rename that
mv ~/.rstudio-desktop  ~/backup-rstudio-desktop
and (zip and) send along with your e-mail to the experts for diagnosis.


 > On Thu, 9 Feb 2017, Martin Maechler wrote:

 >>
 >>> Ulrik Stervbo  on Thu, 9
 >>> Feb 2017 14:37:57 + writes:
 >>
 >> > Hi Georg, > maybe someone here knows, but I think you
 >> are more likely to get answers to > Rstudio related
 >> questions with RStudio support: >
 >> https://support.rstudio.com/hc/en-us
 >>
 >> > Best, > Ulrik
 >>
 >> Indeed, thank you, Ulrik.
 >>
 >> In this special case, however, I'm quite sure many
 >> readers of R-help would be interested in the answer; so
 >> once you receive an answer, please post it (or a link to
 >> a public URL with it) here on R-help, thank you in
 >> advance.
 >>
 >> We would like to be able to *save*, or sometimes *set* /
 >> *reset* such options "in a scripted manner", e.g. for
 >> controlled exam sessions.
 >>
 >> Martin Maechler, ETH Zurich
 >>
 >> > On Thu, 9 Feb 2017 at 12:35 
 >> wrote:
 >>
 >> >> Hi All, >> I would like to make a backup of my RStudio
 >> IDE options I configure using >> "Tools/Global Options"
 >> from the menu bar. Searching the >> web did not reveal
 >> anything.
 >>
 >> >> Can you tell me where RStudio IDE does store its
 >> configuration?
 >>
 >> >> Kind regards >> Georg
 >>
 >> __
 >> R-help@r-project.org mailing list -- To UNSUBSCRIBE and
 >> more, see https://stat.ethz.ch/mailman/listinfo/r-help
 >> PLEASE do read the posting guide
 >> http://www.R-project.org/posting-guide.html and provide
 >> commented, minimal, self-contained, reproducible code.
 >>

 >
---
 > Jeff Newmiller The .  .  Go Live...
 > DCN: Basics: ##.#.  ##.#.  Live
 > Go...  Live: OO#.. Dead: OO#..  Playing Research Engineer
 > (Solar/Batteries O.O#.  #.O#.  with /Software/Embedded
 > Controllers) .OO#.  .OO#.  rocks...1k
 >



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

Re: [R] Location of HTML help files

[Keith Ponting]

On Wed, Jul 16, 2008 at 7:27 PM, Jan Smit  wrote:
> I am using R 2.7.1 under Windows XP, installed at C:/R/R-2.7.1.
>
> The location of the HTML SearchEngine is
> file:///C:/R/R-2.7.1/doc/html/search/SearchEngine.html. Now, when I
type 
a
> phrase, say "reshape", in the search text field, the Search Results
page
> suggest that the location of the reshape HTML help file is
> file:///C:/R/library/stats/html/reshape.html, while in reality it is
> file:///C:/R/R-2.7.1/library/stats/html/reshape.html.
>
> Is there an easy way in which I can fix this?

I too had this problem with Firefox 3.0 and 3.0.1 (but not with 3.0 RC1
by the way).
A work-around which works for me is to go directly to the search engine
page (enter URI
file:///C:/Program%20Files%20(x86)/R/R-2.7.1/doc/html/search/SearchEngin
e.html on my Windows Vista installation), rather than going there by
following the link on the R documentation page
(file:///C:/Program%20Files%20(x86)/R/R-2.7.1/doc/html/index.html) 

Keith Ponting
Aurix Ltd, Malvern WR14 3SZ  UK

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

[Keith Jewell]

Hi everyone,

I thought that for a selfStart function, these two should be exactly 
equivalent
> nls(Aform, DF)
> nls(Aform, DF, start=getInitial(Aform, DF))
but in this example that is not the case in R (although it is in S-plus 
V6.2)
--
SSbatch<-selfStart(
model=function(Batch, Coeffs)
{
Coeffs[Batch]

}
,initial=function(mCall, data, LHS)
{
 # Estimate coefficients as mean of each batch
xy <- sortedXyData(mCall[["Batch"]], LHS, data)
 Batch <- data[[as.character(mCall[["Batch"]])]]
 # check Batch is successive integers starting at 1
 if ((min(xy$x) !=1) | (any(diff(xy$x)!=1))) stop(
   "Batch is not a successive integers sequence")
 Lval <- list(xy$y)
 names(Lval) <- mCall["Coeffs"]
 Lval
}
)
DF <- data.frame(A=c(0.9, 1.1, 1.9, 2.0, 2.1, 2.9, 3.0), 
Batch=c(1,1,2,2,2,3,3))
Aform <- formula(A~SSbatch(Batch,cA))
nls(Aform, DF, start=getInitial(Aform, DF))
nls(Aform, DF)

Don't ask why I'd want such a silly selfStart, that's a long story.
I guess wherever I would have used nls(Aform, DF)
I could use nls(Aform, DF, start=getInitial(Aform, DF))
but that seems clumsy.

Can anyone point out my mistake? Or is this a limitation of nls in R (I 
hesitate to use the b*g word).

Thanks in advance,

Keith Jewell
--

I don't think it's relevant but, for completeness:

> sessionInfo()

 version 2.7.0 (2008-04-22)
i386-pc-mingw32

locale:
LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United 
Kingdom.1252;LC_MONETARY=English_United 
Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

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

other attached packages:
[1] xlsReadWrite_1.3.2 svSocket_0.9-5 svIO_0.9-5 R2HTML_1.58 
svMisc_0.9-5   svIDE_0.9-5

loaded via a namespace (and not attached):
[1] tools_2.7.0 VGAM_0.7-7

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

Re: [R] Unexpected nls behaviour: Solved

[Keith Jewell]

Hi Everyone,

I'd omitted the non-optional 'parameters' argument to selfStart. Making this 
change to SSbatch gives the same (successful) result from the two calls to 
nls.

SSbatch<-selfStart(
model=function(Batch, Coeffs)
{
Coeffs[Batch]
}
,initial=function(mCall, data, LHS)
{
 # Estimate coefficients as mean of each batch
xy <- sortedXyData(mCall[["Batch"]], LHS, data)
 Batch <- data[[as.character(mCall[["Batch"]])]]
 # check Batch is successive integers starting at 1
 if ((min(xy$x) !=1) | (any(diff(xy$x)!=1))) stop(
   "Batch is not a successive integers sequence")
 Lval <- list(xy$y)
 names(Lval) <- mCall["Coeffs"]
 Lval
}
, parameters = c("Coeffs")
)

Sorry for wasting anyones time.

Keith Jewell
-
"Keith Jewell" <[EMAIL PROTECTED]> wrote in message news:...
> Hi everyone,
>
> I thought that for a selfStart function, these two should be exactly 
> equivalent
>> nls(Aform, DF)
>> nls(Aform, DF, start=getInitial(Aform, DF))
> but in this example that is not the case in R (although it is in S-plus 
> V6.2)
> --
> SSbatch<-selfStart(
> model=function(Batch, Coeffs)
> {
> Coeffs[Batch]
>
> }
> ,initial=function(mCall, data, LHS)
> {
> # Estimate coefficients as mean of each batch
>xy <- sortedXyData(mCall[["Batch"]], LHS, data)
> Batch <- data[[as.character(mCall[["Batch"]])]]
> # check Batch is successive integers starting at 1
> if ((min(xy$x) !=1) | (any(diff(xy$x)!=1))) stop(
>   "Batch is not a successive integers sequence")
> Lval <- list(xy$y)
> names(Lval) <- mCall["Coeffs"]
> Lval
> }
> )
> DF <- data.frame(A=c(0.9, 1.1, 1.9, 2.0, 2.1, 2.9, 3.0), 
> Batch=c(1,1,2,2,2,3,3))
> Aform <- formula(A~SSbatch(Batch,cA))
> nls(Aform, DF, start=getInitial(Aform, DF))
> nls(Aform, DF)
> 
> Don't ask why I'd want such a silly selfStart, that's a long story.
> I guess wherever I would have used nls(Aform, DF)
> I could use nls(Aform, DF, start=getInitial(Aform, DF))
> but that seems clumsy.
>
> Can anyone point out my mistake? Or is this a limitation of nls in R (I 
> hesitate to use the b*g word).
>
> Thanks in advance,
>
> Keith Jewell
> --
>
> I don't think it's relevant but, for completeness:
>
>> sessionInfo()
>
> version 2.7.0 (2008-04-22)
> i386-pc-mingw32
>
> locale:
> LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United 
> Kingdom.1252;LC_MONETARY=English_United 
> Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252
>
> attached base packages:
> [1] stats graphics  grDevices datasets  tcltk utils methods 
> base
>
> other attached packages:
> [1] xlsReadWrite_1.3.2 svSocket_0.9-5 svIO_0.9-5 R2HTML_1.58 
> svMisc_0.9-5   svIDE_0.9-5
>
> loaded via a namespace (and not attached):
> [1] tools_2.7.0 VGAM_0.7-7
>  
>
>

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

[Keith Jewell]

Hi All,

I've hit a problem using nls. I think it may be a restriction in the 
applicability of nls and I may have found a fix, but I've been wrong before.

This example is simplified to the essentials. My real application is much 
more complicated.

Take a function of matrix 'x' with additional arguments:
matrix 'aMat' whose values are _not_ to be determined by nls
vector 'Coeffs' whose vales _are_ to be determined.
For simplicity, this isn't a selfStart function with an 'initial' attribute, 
but that doesn't change things.

Myfunc<-function(x, aMat, Coeffs)
{
#
# result = quadratic response in x with
# terms selected by aMat
#
aMat[aMat!=0] <- Coeffs
rowSums((x%*%aMat)%*%t(x))
}

If aMat is passed in by name (e.g. aMat = bMat) nls fails.
e.g.
#
# data frame with some noise
DF <- data.frame(x1 = runif(20, 1, 20), x2=runif(20, 1, 20))
DF$y <- 1 +DF$x1 +DF$x2 +DF$x1*DF$x2 +DF$x1^2 + DF$x2^2 + rnorm(20)
#
# matrix to pass in as aMat
bMat <- matrix(c(1,1,0,0), 2, 2)
#
# and nls fails
nls(y ~ Myfunc(cbind(x1, x2), bMat, aVec), DF, start=list(aVec=c(1,2)))
#
# pass in the same matrix other than by name and it works
nls(y ~ Myfunc(cbind(x1, x2), matrix(c(1,1,0,0), 2, 2), aVec), DF, 
start=list(aVec=c(1,2)))

I think the problem lies in this line in nls

  for (var in varNames[!varIndex]) mf[[var]] <- eval(as.name(var), data, 
env)

This adds values for some named arguments (bMat) as columns of the data 
frame. The problem is that generally they don't have the same number of 
rows. I've made it work for my example by replacing that line with this 
line, which adds values for those arguments to the data frame as parameters 
rather than as a column

  attributes(mf)[["parameters"]] <- 
c(attributes(mf)[["parameters"]],lapply(varNames[!varIndex], function(var) 
eval(as.name(var), data, env)) )

Problem is, I really don't know nls internals enough to be sure I haven't 
broken something.
And anyway, if this is really an improvement I ought to share it, but don't 
know how.

Or I could have totally the wrong end of the stick...

Comments, corrections and advice are welcome.

Thanks in advance,

Keith Jewell
---
> sessionInfo()
R version 2.7.0 (2008-04-22)
i386-pc-mingw32

locale:
LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United 
Kingdom.1252;LC_MONETARY=English_United 
Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

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

other attached packages:
[1] xlsReadWrite_1.3.2 svSocket_0.9-5 svIO_0.9-5 R2HTML_1.58
[5] svMisc_0.9-5   svIDE_0.9-5

loaded via a namespace (and not attached):
[1] tools_2.7.0 VGAM_0.7-7

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

Re: [R] Opening R from Tinn without setting directory each time

[Keith Jewell]

Hi,

See http://sourceforge.net/forum/forum.php?thread_id=1741502&forum_id=481901

I'm running under Windows Server 2003, Standard x64 Edition and (I think) I 
can't make any "system wide" registry changes. In my R installation there 
are (I think) no relevant registry entries. I think that may be the reason 
that:
a) The workaround mentioned in that SouirceForge thread (start R before 
Tinn-R) doesn't work for me.
b) The fix of running RSetReg.exe doesn't work for me.

As far as I know, the new version of Tinn-R (1.19.5.0) isn't released yet.
So I still have to redefine the R location each time I start Tinn-R (but 
it's still worth it!).

Hope that helps someone

Keith Jewell

"bartjoosen" <[EMAIL PROTECTED]> wrote in message 
news:[EMAIL PROTECTED]
>
> If you first install a newer version of R and aftwards remove a previous
> version, Tinn-R gives this behaviour. Doesn't matter if you activate the
> option to register.
>
> Best regards
>
> Bart
>
>
> Philippe Grosjean wrote:
>>
>>
>> Paul Chatfield wrote:
>>> Hi - someone has just e-mailed me direct with the answer which it'd be
>>> helpful to paste just so future users who have the same issue can see.
>>> Just
>>> follow the advice below and it works perfectly.
>>>
>>> Open a command window (Run;cmd) and cd to the bin directory of your R
>>> installation (cd C:/Program Files).  Run the program RSetReg.exe and
>>> that's it, Tinn-R should be able to start R.  When you update R repeat
>>> the
>>> process.
>>
>> This shouldn't be needed if you activate the option to register the
>> installed/upgraded version of R in the registry (somewhere at the end of
>>   the installer's questions).
>> Best,
>>
>> Philippe Grosjean
>>
>>>
>>> Paul Chatfield wrote:
>>>> Hi - I can access R from Tinn-R by going to 
>>>> Options->Main->Application/R
>>>> and setting the search path, but each time I exit Tinn-R I have to
>>>> redefine the search path.  Is there no way of fixing that directory as
>>>> default?  I have installed R under its default directory C:/Program
>>>> Files/R/R-2.7.1 and Tinn under a variety of different places to try to
>>>> rectify the problem though currently under C:/Program Files/Tinn-R. 
>>>> Any
>>>> ideas what I'm missing?
>>>>
>>>> Thanks
>>>>
>>>> Paul

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

Re: [R] Fix for nls bug???

[Keith Jewell]

Dear Prof. Ripley,

Thank you for your helpful reply. I will download and try R-patched ASAP.

I take your point, I should have tried the latest version (R-patched) before 
posting.

With respect to R-patched, would you recommend its use routinely, or only in 
investigation of "unexpected behaviour"?

Thanks again,

Keith Jewell
-
"Prof Brian Ripley" <[EMAIL PROTECTED]> wrote in message 
news:[EMAIL PROTECTED]
> Your example works in R-patched, as a consequence of investigations of a 
> different problem.  (See the comments in the posting guide about updating 
> your R and trying the very latest versions.)
>
> Windows binaries for R-patched are available on CRAN.
>
> On Tue, 5 Aug 2008, Keith Jewell wrote:
>
>> Hi All,
>>
>> I've hit a problem using nls. I think it may be a restriction in the
>> applicability of nls and I may have found a fix, but I've been wrong 
>> before.
>>
>> This example is simplified to the essentials. My real application is much
>> more complicated.
>>
>> Take a function of matrix 'x' with additional arguments:
>> matrix 'aMat' whose values are _not_ to be determined by nls
>> vector 'Coeffs' whose vales _are_ to be determined.
>> For simplicity, this isn't a selfStart function with an 'initial' 
>> attribute,
>> but that doesn't change things.
>>
>> Myfunc<-function(x, aMat, Coeffs)
>> {
>> #
>> # result = quadratic response in x with
>> # terms selected by aMat
>> #
>> aMat[aMat!=0] <- Coeffs
>> rowSums((x%*%aMat)%*%t(x))
>> }
>>
>> If aMat is passed in by name (e.g. aMat = bMat) nls fails.
>> e.g.
>> #
>> # data frame with some noise
>> DF <- data.frame(x1 = runif(20, 1, 20), x2=runif(20, 1, 20))
>> DF$y <- 1 +DF$x1 +DF$x2 +DF$x1*DF$x2 +DF$x1^2 + DF$x2^2 + rnorm(20)
>> #
>> # matrix to pass in as aMat
>> bMat <- matrix(c(1,1,0,0), 2, 2)
>> #
>> # and nls fails
>> nls(y ~ Myfunc(cbind(x1, x2), bMat, aVec), DF, start=list(aVec=c(1,2)))
>> #
>> # pass in the same matrix other than by name and it works
>> nls(y ~ Myfunc(cbind(x1, x2), matrix(c(1,1,0,0), 2, 2), aVec), DF,
>> start=list(aVec=c(1,2)))
>>
>> I think the problem lies in this line in nls
>>
>>  for (var in varNames[!varIndex]) mf[[var]] <- eval(as.name(var), data,
>> env)
>>
>> This adds values for some named arguments (bMat) as columns of the data
>> frame. The problem is that generally they don't have the same number of
>> rows. I've made it work for my example by replacing that line with this
>> line, which adds values for those arguments to the data frame as 
>> parameters
>> rather than as a column
>>
>>  attributes(mf)[["parameters"]] <-
>> c(attributes(mf)[["parameters"]],lapply(varNames[!varIndex], 
>> function(var)
>> eval(as.name(var), data, env)) )
>>
>> Problem is, I really don't know nls internals enough to be sure I haven't
>> broken something.
>> And anyway, if this is really an improvement I ought to share it, but 
>> don't
>> know how.
>>
>> Or I could have totally the wrong end of the stick...
>>
>> Comments, corrections and advice are welcome.
>>
>> Thanks in advance,
>>
>> Keith Jewell
>> ---
>>> sessionInfo()
>> R version 2.7.0 (2008-04-22)
>> i386-pc-mingw32
>>
>> locale:
>> LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United
>> Kingdom.1252;LC_MONETARY=English_United
>> Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252
>>
>> attached base packages:
>> [1] stats graphics  grDevices datasets  tcltk utils methods
>> base
>>
>> other attached packages:
>> [1] xlsReadWrite_1.3.2 svSocket_0.9-5 svIO_0.9-5 R2HTML_1.58
>> [5] svMisc_0.9-5   svIDE_0.9-5
>>
>> loaded via a namespace (and not attached):
>> [1] tools_2.7.0 VGAM_0.7-7
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide 
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
> -- 
> Brian D. Ripley,  [EMAIL PROTECTED]
> Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
> University of Oxford, Tel:  +44 1865 272861 (self)
> 1 South Parks Road, +44 1865 272866 (PA)
> Oxford OX1 3TG, UKFax:  +44 1865 272595
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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

Re: [R] Fix for nls bug??? [already fixed in R-patched]

[Keith Jewell]

Dear Prof. Ripley,

Thanks for that.

Just to wrap up the thread, I confirm that my problem is fully fixed in 
R-patched.

Best regards,

Keith Jewell
---
"Prof Brian Ripley" <[EMAIL PROTECTED]> wrote in message 
news:[EMAIL PROTECTED]
> On Fri, 8 Aug 2008, Keith Jewell wrote:
>
>> Dear Prof. Ripley,
>>
>> Thank you for your helpful reply. I will download and try R-patched ASAP.
>>
>> I take your point, I should have tried the latest version (R-patched) 
>> before
>> posting.
>>
>> With respect to R-patched, would you recommend its use routinely, or only 
>> in
>> investigation of "unexpected behaviour"?
>
> I use it routinely.
>
> In this case I was not expecting this to be fixed, but just tried your 
> example in my usual R-patched (when it worked) and then in 2.7.1 (to 
> confirm I could reproduce it) -- I then had to figure out what fix also 
> fixed this one -- as it was
>
> o   nls() was only finding its 'weights' argument in the case when
> all the variables in the formula were of the same length and
> hence that model.frame() could be used.
>
> it would have been obvious to no one.  Hence the advice to try R-patched 
> (rather than just read about the changes) can pay off.
>
>>
>> Thanks again,
>>
>> Keith Jewell
>> -
>> "Prof Brian Ripley" <[EMAIL PROTECTED]> wrote in message
>> news:[EMAIL PROTECTED]
>>> Your example works in R-patched, as a consequence of investigations of a
>>> different problem.  (See the comments in the posting guide about 
>>> updating
>>> your R and trying the very latest versions.)
>>>
>>> Windows binaries for R-patched are available on CRAN.
>>>
>>> On Tue, 5 Aug 2008, Keith Jewell wrote:
>>>
>>>> Hi All,
>>>>
>>>> I've hit a problem using nls. I think it may be a restriction in the
>>>> applicability of nls and I may have found a fix, but I've been wrong
>>>> before.
>>>>
>>>> This example is simplified to the essentials. My real application is 
>>>> much
>>>> more complicated.
>>>>
>>>> Take a function of matrix 'x' with additional arguments:
>>>> matrix 'aMat' whose values are _not_ to be determined by nls
>>>> vector 'Coeffs' whose vales _are_ to be determined.
>>>> For simplicity, this isn't a selfStart function with an 'initial'
>>>> attribute,
>>>> but that doesn't change things.
>>>>
>>>> Myfunc<-function(x, aMat, Coeffs)
>>>> {
>>>> #
>>>> # result = quadratic response in x with
>>>> # terms selected by aMat
>>>> #
>>>> aMat[aMat!=0] <- Coeffs
>>>> rowSums((x%*%aMat)%*%t(x))
>>>> }
>>>>
>>>> If aMat is passed in by name (e.g. aMat = bMat) nls fails.
>>>> e.g.
>>>> #
>>>> # data frame with some noise
>>>> DF <- data.frame(x1 = runif(20, 1, 20), x2=runif(20, 1, 20))
>>>> DF$y <- 1 +DF$x1 +DF$x2 +DF$x1*DF$x2 +DF$x1^2 + DF$x2^2 + rnorm(20)
>>>> #
>>>> # matrix to pass in as aMat
>>>> bMat <- matrix(c(1,1,0,0), 2, 2)
>>>> #
>>>> # and nls fails
>>>> nls(y ~ Myfunc(cbind(x1, x2), bMat, aVec), DF, start=list(aVec=c(1,2)))
>>>> #
>>>> # pass in the same matrix other than by name and it works
>>>> nls(y ~ Myfunc(cbind(x1, x2), matrix(c(1,1,0,0), 2, 2), aVec), DF,
>>>> start=list(aVec=c(1,2)))
>>>>
>>>> I think the problem lies in this line in nls
>>>>
>>>>  for (var in varNames[!varIndex]) mf[[var]] <- eval(as.name(var), data,
>>>> env)
>>>>
>>>> This adds values for some named arguments (bMat) as columns of the data
>>>> frame. The problem is that generally they don't have the same number of
>>>> rows. I've made it work for my example by replacing that line with this
>>>> line, which adds values for those arguments to the data frame as
>>>> parameters
>>>> rather than as a column
>>>>
>>>>  attributes(mf)[["parameters"]] <-
>>>> c(attributes(mf)[["parameters"]],lapply(varNames[!varIndex],
>>>> function(var)
>>>> eval(as.name(var), data, env)) )
>>>>
>>>> Problem is, I really don't know nls interna

[R] "nested" getInitial calls; variable scoping problems

[Keith Jewell]

Hi All,

Another nls related problem (for background, I'm migrating a complicated 
modelling package from S-plus to R).

Below I've reduced this to the minimum necessary to demonstrate my problem 
(I think); the real situation is more complicated.

Two similar selfStart functions, ssA and ssB.
The 'initial' function for ssB modifies its arguments a little and then 
calls getInital for ssA.
In addition to the "x" and the fitted coefficients ("Coeff"), ssA and ssB 
have arguments ("A") which are not the same length as the dataframe, so 
cannot be passed as columns of that dataframe.

The initial function for ssA uses eval(... parent.frame()) to find "A", 
which is fine when called from .GlobalEnv. But when called from the initial 
function of ssB it can't find (the modified version of) "A" .

If A is assigned as a parameter of the dataframe then getInitial.formula 
returns "A" as initial estimates of "Coeff", so that doesn't work.
I've explored the evaluation frame parent-child structure created by the 
nested getInitial calls, but it doesn't seem helpful, and I certainly 
couldn't trust my understanding of it.
I'm considering making up the matched call and calling the 'initial' 
function of ssA directly, in the same manner as getInitial.selfStart
attr(ssA, "initial"))(mCall = mCall, data = data, LHS = LHS)
or perhaps if I call getInitial thus...
getInitial(ssA, data, mCall, LHS = NULL, ...)
it will call getInitial.selfStart without calling getInitial.formula.
In that case perhaps eval(... parent.frame())  will find the argument?
Or, at least I could attach the argument to the dataframe as a parameter.

But this all seems a but clumsy, and I feel there must be a better way. Any 
comments and/or advice will be welcome.

Thanks in advance;

Keith Jewell
---
code showing the problem:
ssA <- selfStart(
  model = function(x, Coeff, A)
{
paste(x, Coeff, A)
},
  initial = function(mCall, data, LHS)
{
x <- eval(mCall[["x"]], data, parent.frame())
A <- eval(mCall[["A"]], data, parent.frame())
paste("CoeffA", x, A)
},
  parameters = c("Coeff")
   )
ssB <- selfStart(
  model = function(x, Coeff, A)
{
paste(x, Coeff, A)
},
  initial = function(mCall, data, LHS)
{
x <- eval(mCall[["x"]], data, parent.frame())
A <- eval(mCall[["A"]], data, parent.frame())
Amod <- paste(A, "mod in B")
getInitial(y ~ ssA(x, Coeff, Amod), data)
},
  parameters = c("Coeff")
 )
getInitial(y ~ ssA("this", "that", "other"), data= data.frame(x=c("test")))
getInitial(y ~ ssB("this", "that", "other"), data= data.frame(x=c("test")))
---
> sessionInfo()
R version 2.7.1 Patched (2008-08-15 r46352)
i386-pc-mingw32

locale:
LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United 
Kingdom.1252;LC_MONETARY=English_United 
Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

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

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

Re: [R] "nested" getInitial calls; variable scoping problems: Solved??

[Keith Jewell]

Hi All,

I've found a solution to my problems which I think is optimal, but I've been 
wrong before (!!) so comments are welcome.

To avoid scoping problems in the 'inner' getInitial I modify the formula so 
that all variable values are specified explitly in the formula, the inner 
getInitial need not search for any values so there are no scoping problems:
---
ssA <- selfStart(
  model = function(x, Coeff, A)
{
paste(x, Coeff, A)
},
  initial = function(mCall, data, LHS)
{
x <- eval(mCall[["x"]], data, parent.frame())
A <- eval(mCall[["A"]], data, parent.frame())
paste("CoeffA", x, A)
},
  parameters = c("Coeff")
   )
ssB <- selfStart(
  model = function(x, Coeff, A)
{
paste(x, Coeff, A)
},
  initial = function(mCall, data, LHS)
{
A <- eval(mCall[["A"]], data, parent.frame())
Amod <- paste(A, "mod in B")
#---
# Doing a getInitial for (for example)  y' ~ ssA(x, Coeff, Amod)
# to avoid scoping problems, specify values in the formula for LHS and
# all ssA arguments except Coeff; i.e. all except attr(ssA, "pnames")
#--
object <- y ~ ssA(x, Coeff, A)# define a formula object
object[[2]] <- eval(LHS, data, parent.frame())# left hand side 
values
# name RHS arguments by match.call, then assign values  to x and A
object[[3]] <- match.call(get(as.character(object[[3]][[1]])), 
object[[3]])
object[[3]]$x <-  as.character(eval(mCall[["x"]], data, parent.frame())) 
# x
object[[3]]$A <-  Amod  # A
getInitial(object, data) # do the getInitial
 },
  parameters = c("Coeff")
 )
getInitial(y ~ ssA("this", "that", "other"), data= data.frame(x=c("test"), 
y=1))
getInitial(y ~ ssB("this", "that", "other"), data= data.frame(x=c("test"), 
y=1))
getInitial(y ~ ssB(x, "that", "other"), data= data.frame(x=c("test1", 
"test2"), y=1:2))
--

Hope that helps someone,

Keith Jewell

"Keith Jewell" <[EMAIL PROTECTED]> wrote in message 
news:[EMAIL PROTECTED]
> Hi All,
>
> Another nls related problem (for background, I'm migrating a complicated 
> modelling package from S-plus to R).
>
> Below I've reduced this to the minimum necessary to demonstrate my problem 
> (I think); the real situation is more complicated.
>
> Two similar selfStart functions, ssA and ssB.
> The 'initial' function for ssB modifies its arguments a little and then 
> calls getInital for ssA.
> In addition to the "x" and the fitted coefficients ("Coeff"), ssA and ssB 
> have arguments ("A") which are not the same length as the dataframe, so 
> cannot be passed as columns of that dataframe.
>
> The initial function for ssA uses eval(... parent.frame()) to find "A", 
> which is fine when called from .GlobalEnv. But when called from the 
> initial function of ssB it can't find (the modified version of) "A" .
>
> If A is assigned as a parameter of the dataframe then getInitial.formula 
> returns "A" as initial estimates of "Coeff", so that doesn't work.
> I've explored the evaluation frame parent-child structure created by the 
> nested getInitial calls, but it doesn't seem helpful, and I certainly 
> couldn't trust my understanding of it.
> I'm considering making up the matched call and calling the 'initial' 
> function of ssA directly, in the same manner as getInitial.selfStart
>attr(ssA, "initial"))(mCall = mCall, data = data, LHS = LHS)
> or perhaps if I call getInitial thus...
>getInitial(ssA, data, mCall, LHS = NULL, ...)
> it will call getInitial.selfStart without calling getInitial.formula.
> In that case perhaps eval(... parent.frame())  will find the argument?
> Or, at least I could attach the argument to the dataframe as a parameter.
>
> But this all seems a but clumsy, and I feel there must be a better way. 
> Any comments and/or advice will be welcome.
>
> Thanks in advance;
>
> Keith Jewell
> ---
> code showing the problem:
> ssA <- selfStart(
>  model = function(x, Coeff, A)
>{
>paste(x, Coeff, A)
>},
>  initial = function(mCall, data, LHS)
>{
>x <- eval(mCall[["x"]], data, parent.frame())
>A <- eval(mCall[["A"]], data, parent.frame())
>paste("CoeffA", x, A)
>},
>  parameters = c("Coeff"

Re: [R] Upgrade 'R'

[Keith Jewell]

As possible help to others, and also as a request for comments on how I 
might do things better, I describe how I've recently altered my system to 
handle this.

I'm on a Windows Server 2003 network and the R installation is accessible to 
many others. Everyone has "read" access to all installation files, but only 
I have write access. I do _not_ have Administrator privileges on the server, 
so I cannot make/change registry entries.

R versions are in an R folder tree, which also holds the version independent 
library folder.

//Server02/stats/R
//Server02/stats/R/R-2.7.0
//Server02/stats/R/R-2.7.1
//Server02/stats/R/R-2.7.1pat
//Server02/stats/R/R-2.7.2
   :
//Server02/stats/R/library

In each version I have edited /etc/Rprofile.site to include the line
 .libPaths("//Server02/stats/R/library")

The "default" libraries (base, boot, class...) are installed into the 
relevant version specific .../R-n.n.n/library/ folder by the windows 
installer program (e.g. R-2.7.2-win32.exe).
Occasionaly, and after installing a new R-version,  I update all the 
downloaded libraries in the version independent //Server02/stats/R/library/ 
folder with a simple update.packages().

HTH. Comments welcome.

Keith Jewell
---
"Leon Yee" <[EMAIL PROTECTED]> wrote in message 
news:[EMAIL PROTECTED]
> Hello, Kevin
>
>You can get some hints by browsing in this mailist with the subject of 
> " Upgrading R means I lose my packages", which were posted several days 
> ago.
>
> HTH
>
> Leon
>
>
> [EMAIL PROTECTED] wrote:
>> More and more I am getting warnings from packages that I install that the 
>> package was built with 2.7.2 (I am running 2.7.1). I would like to 
>> upgrade but don't want to loose all of the packages that I have installed 
>> and the settings. Is there a way to just "upgrade" without uninstalling 
>> and reinstalling 'R'?
>>
>> Thank you.
>>
>> Kevin

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

Re: [R] Upgrade 'R'

[Keith Jewell]

"Prof Brian Ripley" <[EMAIL PROTECTED]> wrote in message 
news:[EMAIL PROTECTED]
> On Tue, 2 Sep 2008, Keith Jewell wrote:
>
>> As possible help to others, and also as a request for comments on how I
>> might do things better, I describe how I've recently altered my system to
>> handle this.
>>
>> I'm on a Windows Server 2003 network and the R installation is accessible 
>> to
>> many others. Everyone has "read" access to all installation files, but 
>> only
>> I have write access. I do _not_ have Administrator privileges on the 
>> server,
>> so I cannot make/change registry entries.
>>
>> R versions are in an R folder tree, which also holds the version 
>> independent
>> library folder.
>>
>> //Server02/stats/R
>> //Server02/stats/R/R-2.7.0
>> //Server02/stats/R/R-2.7.1
>> //Server02/stats/R/R-2.7.1pat
>> //Server02/stats/R/R-2.7.2
>>   :
>> //Server02/stats/R/library
>>
>> In each version I have edited /etc/Rprofile.site to include the line
>> .libPaths("//Server02/stats/R/library")
>>
>> The "default" libraries (base, boot, class...) are installed into the
>> relevant version specific .../R-n.n.n/library/ folder by the windows
>> installer program (e.g. R-2.7.2-win32.exe).
>> Occasionaly, and after installing a new R-version,  I update all the
>> downloaded libraries in the version independent 
>> //Server02/stats/R/library/
>> folder with a simple update.packages().
>>
>> HTH. Comments welcome.
>
> Please read the rw-FAQ for more details and insight, especially in using 
> Renviron.site (here you want to set R_LIBS_SITE: see ?libPaths) and using 
> update.packages(checkBuilt=TRUE, ask=FALSE)
>
> For example, my sysadmins have Windows R installed locally (via the MSI 
> installer), and then have on each machine in etc/Renviron.site have
>
> R_LIBS_SITE=N:/R/library/2.7
> R_LIBS_USER=P:/R/win-library/2.7
>
> This allows both common packages and user-specific packages to be stored 
> on SMB-mounted drives.  (N: is common to all machines, P: is 'personal' --  
> we use mapped drives rather than shares to allow quick changes of server)
>
> There are two reasons we install locally.  The first is performance: if 20 
> machines in a lab start R up simultaneously the network load is high.  The 
> second is that our security settings (and I believe the defaults these 
> days) disallow use of CHM help on network drives -- we remove the chtml 
> directories from packages installed on R_LIBS_SITE, so R defaults to text 
> help for those packages.
>

Thanks Professor Ripley. I'd looked at the R Windows FAQ but  not carefully
enough.Using Renviron.site does seem much more elegant than Rprofile.site.
I've now included a line in Renviron.site
R_LIBS_SITE=//Server02/stats/R/library

I also note your use of site and personal libraries specific to "first 
decimal"
versions of R. Would you recommend others to follow this practice?

Thanks again,

Keith Jewell.

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

[keith parker]

Is anyone aware of R code for calculating sample size in the context of
tolerance intervals? Thanks, Keith

[[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 should a SelfStart function handle illegal parameter values?

[Keith Jewell]

Hi Everyone,

I'm trying to write selfStart non-linear models for use with nls. In these 
models some combinations of parameter values are illegal; the function value 
is undefined.

That's OK when calling the function directly [e.g.  SSmodel(x, pars...)]; I 
return an appropriate non-value such as NA or Inf.

However, when called from nls [e.g. nls(y~SSmodel(x, pars...), ...)] those 
non-values lead to errors such as (but not limited to):
Error in numericDeriv(form[[3L]], names(ind), env) :
  Missing value or an infinity produced when evaluating the model
or (if I provide a gradient attribute)
Error in qr.default(.swts * attr(rhs, "gradient")) :
  NA/NaN/Inf in foreign function call (arg 1)

I can't see a way of making nls either stick to legal parameter values, or 
accept NA/NaN/Inf as indicating "bad" parameter values.

I really do want to use nls rather than a bounded optimisation tool (such as 
optim) because this fits into a much bigger picture predicated on nls.

I'd appreciate any suggestions.

Keith Jewell

> sessionInfo()
R version 2.9.1 (2009-06-26)
i386-pc-mingw32

locale:
LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United 
Kingdom.1252;LC_MONETARY=English_United 
Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

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

other attached packages:
[1] xlsReadWrite_1.3.3 svSocket_0.9-43svMisc_0.9-48  TinnR_1.0.3 
R2HTML_1.59-1  Hmisc_3.6-1

loaded via a namespace (and not attached):
[1] cluster_1.12.0  grid_2.9.1  lattice_0.17-25 stats4_2.9.1 
VGAM_0.7-9

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 should a SelfStart function handle illegal parameter values?

[Keith Jewell]

Hi Everyone,

I'm trying to write selfStart non-linear models for use with nls. In these 
models some combinations of parameter values are illegal; the function value 
is undefined.

That's OK when calling the function directly [e.g.  SSmodel(x, pars...)]; I 
return an appropriate non-value such as NA or Inf.

However, when called from nls [e.g. nls(y~SSmodel(x, pars...), ...)] those 
non-values lead to errors such as (but not limited to):
Error in numericDeriv(form[[3L]], names(ind), env) :
  Missing value or an infinity produced when evaluating the model
or (if I provide a gradient attribute)
Error in qr.default(.swts * attr(rhs, "gradient")) :
  NA/NaN/Inf in foreign function call (arg 1)

A toy example demonstrating my problem (legal values of param are >1):
#---
SSexample<-selfStart(
 model=function(x, param) x^log(param-1),
 initial = function(mCall, data, LHS){
   val<- 1.001
   names(val) <- mCall[c("param")]
   val
   },
 parameters=c("param")
)
#
nls(y~SSexample(x, par), data=data.frame(x=1:10,y=rnorm(10)))
#-

(repeat the last line a few times and you'll get the error).

I can't see a way of making nls either stick to legal parameter values, or 
accept NA/NaN/Inf as indicating "bad" parameter values.

I really do want to use nls rather than a bounded optimisation tool (such as 
optim) because this fits into a much bigger picture predicated on nls.

I'd appreciate any suggestions.

Keith Jewell

> sessionInfo()
R version 2.9.1 (2009-06-26)
i386-pc-mingw32

locale:
LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United
Kingdom.1252;LC_MONETARY=English_United
Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

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

other attached packages:
[1] xlsReadWrite_1.3.3 svSocket_0.9-43svMisc_0.9-48  TinnR_1.0.3
R2HTML_1.59-1  Hmisc_3.6-1

loaded via a namespace (and not attached):
[1] cluster_1.12.0  grid_2.9.1  lattice_0.17-25 stats4_2.9.1
VGAM_0.7-9

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 use nls when [selfStart] function returns NA or Inf??

[Keith Jewell]

Hi Everyone,

I posted this a couple of weeks ago with no responses. My interface (via 
gmane) seemed a bit flakey at the time, so I'm venturing to repost with some 
additional information.

I'm trying to write selfStart non-linear models for use with nls. In these 
models some combinations of parameter values are illegal; the function value 
is undefined.
That's OK when calling the function directly [e.g.  SSmodel(x, pars...)]; it 
is appropriate to return an appropriate non-value such as NA or Inf.
However, when called from nls [e.g. nls(y~SSmodel(x, pars...), ...)] those 
non-values lead to errors such as (but not limited to):
Error in numericDeriv(form[[3L]], names(ind), env) :
  Missing value or an infinity produced when evaluating the model
or (if I provide a gradient attribute)
Error in qr.default(.swts * attr(rhs, "gradient")) :
  NA/NaN/Inf in foreign function call (arg 1)

A toy example demonstrating my problem (legal values of param are >1):
#---
SSexample<-selfStart(
 model=function(x, param) x^log(param-1),
 initial = function(mCall, data, LHS){
   val<- 1.001
   names(val) <- mCall[c("param")]
   val
   },
 parameters=c("param")
)
#
nls(y~SSexample(x, par), data=data.frame(x=1:10,y=rnorm(10)))
#-

(repeat the last line a few times and you'll get the error).

I can't see a way of making nls either
a) stick to legal parameter values (which I'd have trouble specifying 
anyway), or
b) (ideally) accept NA/NaN/Inf as indicating "bad" parameter values, 
equivalent to very large errors in 'y' values

I really do want to use nls rather than a bounded optimisation tool (such as 
optim) because this fits into a much bigger picture predicated on nls.

I'd appreciate any suggestions.

Keith Jewell
==
sessionInfo() is given below.

I think the toy example above is enough to demonstrate my problem, but in 
case it is relevant (I don't think it is) here is some more info about the 
models I'm fitting:

I'm fitting sigmoidal models to microbial growth over time. The specific 
model giving me problems at the moment is only one of a whole class of such 
models which I need to work with. I specify it here to illustrate that it is 
not always obvious what the bounds on the parameters are.

SSbaranyiBR94<-selfStart(
model=function(Time, y0, ymax, mu, lambda, m = 1, v = mu)
{
 #+
 # From the paper Baranyi J. & Roberts T. A. (1994). A dynamic approach to 
predicting bacterial growth in food. Int J. of Fd. Micro. 23. 277-294
 # Papers equations (6), (5b), (4b)
 # eq 4b: y(Time) = y0 + mu' * A(Time) - ln(1+(exp(m' * mu' * A(Time)) - 
1)/exp(m' * (ymax-y0)))/m'
 # eq 5b: A(Time) = Time + ln((exp(-v * Time) + q0)/(1+ q0))/v
 # from eq 6: q0 = 1/(exp(mu'*lambda)-1)
 #
 #   m' = m*ln(10)
 #  mu' = mu/ln(10)
 #
 # Cited paper gives defaults m = 1 and v = mu
 #-
  m <- m * log(10)
  mu <- mu/log(10)
 .value <- ymax - log(1 + (exp(m * (ymax - y0)) - 1)/exp(m * mu * (Time +
log((exp( - v * Time) + (1/(exp(v * lambda) - 1)))/(1 + (1/(exp(v *
lambda) - 1/v)))/m
   .value
}
, initial = function(mCall, data, LHS)
{
# 
},
parameters=c("y0", "ymax", "mu", "lambda")
)

> sessionInfo()
R version 2.9.1 (2009-06-26)
i386-pc-mingw32

locale:
LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United
Kingdom.1252;LC_MONETARY=English_United
Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

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

other attached packages:
[1] xlsReadWrite_1.3.3 svSocket_0.9-43svMisc_0.9-48  TinnR_1.0.3
R2HTML_1.59-1  Hmisc_3.6-1

loaded via a namespace (and not attached):
[1] cluster_1.12.0  grid_2.9.1  lattice_0.17-25 stats4_2.9.1
VGAM_0.7-9

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 use nls when [selfStart] function returns NA or Inf??

[Keith Jewell]

Thanks Gabor, but my problem isn't finding reasonable starting parameter 
values, it's preventing nls "giving up" when it tries parameter values 
resulting in NA or Inf.

I know queries are often over-specific and the appropriate response is 
"don't start there", so I'm trying to balance between simplifying to 
essentials, and providing enough background. Here's some more background.

I've had this whole system working for several years in S-plus where it's 
routinely and successfully applied to a variety of models and data sets. I'm 
gradually porting it all into R

I'm using several sigmoidal models; the 4-optimisable-parameter Baranyi 
model I mentioned is just one of them. I have algorithms to find starting 
values for each of the sigmoidal models. However, even with reasonable 
starting values, nls sometimes tries "illegal" parameter values resulting in 
a "crash" of the nls process.

Later I fit a more complicated model in which each of the parameters in the 
sigmoidal is a function of several other variables; these models typically 
have dozens of parameters. I use the "simple" sigmoidal fits to generate 
starting values for the complicated model. Even when the simple sigmoidal 
fits have all converged, nls sometimes tries "illegal" parameter values when 
fitting the complicated model, so (even) better starting values for the 
original sigmoidals wouldn't solve my ultimate problem.

I think I've correctly identified my problem as preventing nls "giving up" 
when it tries parameter values resulting in NA or Inf - but I have been 
wrong before :-}

I've considered modifying the sigmoidals to return extreme numeric values 
instead of NA or Inf; but
a) it's non-trivial to choose an appropriate extreme
b) it might "break" numericDeriv
c) it offends me (!) to return a  (wrong!) value when the right answer is NA

Any suggestions/comments gratefully received.

Keith Jewell
===
"Gabor Grothendieck"  wrote in message 
news:971536df0909210923r3fd13fb0we72850bf73232...@mail.gmail.com...

With a small number of parameters just use brute force on grid
to calculate starting values.  See nls2 package.

On Mon, Sep 21, 2009 at 12:17 PM, Keith Jewell  
wrote:
> Hi Everyone,
>
> I posted this a couple of weeks ago with no responses. My interface (via
> gmane) seemed a bit flakey at the time, so I'm venturing to repost with 
> some
> additional information.
>
> I'm trying to write selfStart non-linear models for use with nls. In these
> models some combinations of parameter values are illegal; the function 
> value
> is undefined.
> That's OK when calling the function directly [e.g. SSmodel(x, pars...)]; 
> it
> is appropriate to return an appropriate non-value such as NA or Inf.
> However, when called from nls [e.g. nls(y~SSmodel(x, pars...), ...)] those
> non-values lead to errors such as (but not limited to):
> Error in numericDeriv(form[[3L]], names(ind), env) :
> Missing value or an infinity produced when evaluating the model
> or (if I provide a gradient attribute)
> Error in qr.default(.swts * attr(rhs, "gradient")) :
> NA/NaN/Inf in foreign function call (arg 1)
>
> A toy example demonstrating my problem (legal values of param are >1):
> #---
> SSexample<-selfStart(
> model=function(x, param) x^log(param-1),
> initial = function(mCall, data, LHS){
> val<- 1.001
> names(val) <- mCall[c("param")]
> val
> },
> parameters=c("param")
> )
> #
> nls(y~SSexample(x, par), data=data.frame(x=1:10,y=rnorm(10)))
> #-
>
> (repeat the last line a few times and you'll get the error).
>
> I can't see a way of making nls either
> a) stick to legal parameter values (which I'd have trouble specifying
> anyway), or
> b) (ideally) accept NA/NaN/Inf as indicating "bad" parameter values,
> equivalent to very large errors in 'y' values
>
> I really do want to use nls rather than a bounded optimisation tool (such 
> as
> optim) because this fits into a much bigger picture predicated on nls.
>
> I'd appreciate any suggestions.
>
> Keith Jewell
> ==
> sessionInfo() is given below.
>
> I think the toy example above is enough to demonstrate my problem, but in
> case it is relevant (I don't think it is) here is some more info about the
> models I'm fitting:
>
> I'm fitting sigmoidal models to microbial growth over time. The specific
> model giving me problems at the moment is only one of a whole class of 
> such
> models which I need to work with. I specify it here to illustrate that it 
> is
> 

[R] locfit: max number of predictors?

[Keith Jewell]

Hi All,

In another thread Andy Liaw, who CRAN lists as locfit maintainer; said:

From: "Liaw, Andy" 
To: "Guy Green" ; 
Subject: Re: Alternatives to linear regression with multiple variables
Date: 22 February 2010 17:50

You can try the locfit package, which I believe can handle up to 5
variables.  E.g.,


Looking in the locfit documentation (e.g. 
http://www.stats.bris.ac.uk/R/web/packages/locfit/locfit.pdf) I can't see an 
upper limit on the number of predictors; if it is 5 I'm getting close in one 
of my applications.

Can anyone confirm or deny the existence of a 'crisp' upper limit on the 
number of predictors in locfit?

If it is 5, or thereabouts, can anyone suggest an alternative which can 
handle a few more? (I'm using it for multidimensional interpolation).

Best regards,

Keith Jewell

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

Re: [R] locfit: max number of predictors = 6? How interpolate in 5-10D?

[Keith Jewell]

Thanks for that suggestion

I've investigated a little more using...
y <- rowSums(x) + runif(n)
... just so I had some correlation to play with.

The error I get when it fails is "Invalid what in exvval", which I don't 
understand either!
With n=5e3 it worked with 6 variables but not with 7.

I wasn't sure the error was caused by number of variables rather than 
something else, so I tried with...
n <- 100

I also tried locfit rather than locfit.raw using...
xd <- lapply(1:10, function(x) runif(n))
xd <- as.data.frame(xd)
names(xd) <- paste("x", 1:10, sep="")
y=rowSums(xd)
xd$y <- y
aF <- formula(paste("y ~ lp(",paste(names(xd)[1:6], collapse=","), ")"))
locfit(aF, xd)

Both of these gave the same results, success with 6 variables but not with 
7.

IT APPEARS, the maximum number of predictors is 6, but I don't know locfit 
well, and it may be that other settings would allow more variables.
CAN anyone give a more DEFINITIVE ANSWER?

My current data sets currently reach 5 predictors, and I expect this it 
increase.
 In S-Plus (v6.2.1) I used loess in which "Locally quadratic models may have 
at most 4 predictor variables; locally linear models may have at most 15". 
In R stats::loess allows only "one to four numeric predictors".
I'd assumed (foolishly) that because locfit didn't mention limits, the only 
limits were practical (memory, time,...) - it seems not :-(
I guess I could write something myself, I only need rough interpolation, 
even "straight line" interpolation between nearest neighbours would be OK. 
But at first glance it seems non-trivial with a substantial non-fixed number 
of dimensions (nnclust::nnfind to identify neighbours??), and I don't want 
to re-invent wheels.
Can anyone suggest an ALTERNATIVE route for INTERPOLATION in 5-10 
DIMENSIONS?

Best...
(apologies for capitals, not shouting, just highlighting key points for 
those skimming quickly)

Keith Jewell

"Liaw, Andy"  wrote in message 
news:b10baa7d28d88b45af82813c4a6ffa934ce...@usctmx1157.merck.com...
> Well, I should think there's an obvious (if not elegant) way to test it:
>
> n <- 5e3
> m <- 20
> x <- matrix(runif(n * m), nrow=n)
> y <- rnorm(n)
>
> require(locfit)
> fit <- locfit.raw(x[, 1:10], y)
>
> The code above took a while on my laptop, and ended up giving some error
> I don't understand.  Not sure if the error was caused by insufficient
> sample size, or some inherent limitation.  At least it didn't choke on
> five variables.  However, if all 20 columns of x is used, locfit.raw()
> will choke because it can't compute the dimension of some variable that
> it needs to allocate memory for.
>
> I had vague recollection of reading that "5" is the limit somewhere.
> Unfortunately my copy of Local Regression and Likelihood has been MIA
> for a few years, so I can't check there.  In any case it doesn't seem
> like the number of data points and/or computing power are bigger issue.
>
> Andy
>
>> -Original Message-
>> From: r-help-boun...@r-project.org
>> [mailto:r-help-boun...@r-project.org] On Behalf Of Keith Jewell
>> Sent: Thursday, February 25, 2010 4:11 AM
>> To: r-h...@stat.math.ethz.ch
>> Subject: [R] locfit: max number of predictors?
>>
>> Hi All,
>>
>> In another thread Andy Liaw, who CRAN lists as locfit
>> maintainer; said:
>> 
>> From: "Liaw, Andy" 
>> To: "Guy Green" ; 
>> Subject: Re: Alternatives to linear regression with multiple variables
>> Date: 22 February 2010 17:50
>>
>> You can try the locfit package, which I believe can handle up to 5
>> variables.  E.g.,
>> 
>>
>> Looking in the locfit documentation (e.g.
>> http://www.stats.bris.ac.uk/R/web/packages/locfit/locfit.pdf)
>> I can't see an
>> upper limit on the number of predictors; if it is 5 I'm
>> getting close in one
>> of my applications.
>>
>> Can anyone confirm or deny the existence of a 'crisp' upper
>> limit on the
>> number of predictors in locfit?
>>
>> If it is 5, or thereabouts, can anyone suggest an alternative
>> which can
>> handle a few more? (I'm using it for multidimensional interpolation).
>>
>> Best regards,
>>
>> Keith Jewell
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
> Notice:  This e-mail message, together with any attachme...{{dropped:10}}
>

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

[R] update.packages with UNC library path

[Keith Jewell]

Hi all,

I hit a small snag. Here is my workaround (copied verbatim from my aide 
memoire) in case it helps others. (or anyone knows a better way... ;-)

Best regards,

Keith Jewell.

The site library file is defined (in Renviron.site 
R_LIBS_SITE=//Server02/stats/R/library/%v) via a UNC name something like 
"//Server02/stats/R/library/2.10".

As of now [Feb 2010, R version 2.10.1 (2009-12-14)] the menu Packages|Update 
Packages... [=update.packages(ask='graphics')] fails at the last step when 
it can't copy files to a directory named like that. It works if the site 
library is defined via a mapped drive like "L:\\R\\library\\2.10", but we 
don't want to require all users to have this drive mapped. The workaround is 
thus:

In order to update packages:

a) have an appropriate mapped drive (e.g. Stats on 'server02')

b) start R-Gui from that mapped drive. This will give the site library via 
the UNC and the base library via the mapped drive thus:
> .libPaths()
[1] "//Server02/stats/R/library/2.10" "L:/R/R-Current/library"

c) use .LibPaths(new= ) to add the drive mapped path to the beginning of the 
list
> .libPaths(new=choose.dir()) # navigate to folder on mapped drive
> .libPaths()
[1] "L:\\R\\library\\2.10""//Server02/stats/R/library/2.10" 
"L:/R/R-Current/library"

update.packages(ask='graphics') will now work

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

Re: [R] Error Running TinnR with R

[Keith Jewell]

Many people seem to have trouble defining '.trPaths' which is the set of 
file names which TinnR uses to communicate with R; the user must be able to 
create/write/read these files/folders.

In my Rprofile.site I have the single assingment:
.trPaths <- paste(paste(Sys.getenv("APPDATA"), "\\Tinn-R\\tmp\\", sep=""),
 c("", "search.txt", "objects.txt", "file.r", "selection.r", "block.r", 
"lines.r"), sep="")

At Rstartup this defines a user-specific vector based on the current value 
of the Windows environment variable APPDATA. In my particular case, on this 
particular machine, at the moment, it results in this vector:
> .trPaths
[1] "C:\\Documents and Settings\\jewell\\Application Data\\Tinn-R\\tmp\\"
[2] "C:\\Documents and Settings\\jewell\\Application 
Data\\Tinn-R\\tmp\\search.txt"
[3] "C:\\Documents and Settings\\jewell\\Application 
Data\\Tinn-R\\tmp\\objects.txt"
[4] "C:\\Documents and Settings\\jewell\\Application 
Data\\Tinn-R\\tmp\\file.r"
[5] "C:\\Documents and Settings\\jewell\\Application 
Data\\Tinn-R\\tmp\\selection.r"
[6] "C:\\Documents and Settings\\jewell\\Application 
Data\\Tinn-R\\tmp\\block.r"
[7] "C:\\Documents and Settings\\jewell\\Application 
Data\\Tinn-R\\tmp\\lines.r"

I'm on Windows Server 2003 R2; I don't know if the same construction will 
work on Vista, but I think it probably would, I suggest trying it!

If that doesn't work, I have a few comments on your code which might be 
helpful:
a) you have .trPath <- , not .trPaths <-

b) you are defining a list, not a vector, so .trPath[5] (for example) will 
return a single entry list, not a single character value. That doesn't work 
for me!
> as.list(.trPaths)[5]
[[1]]
[1] "C:\\Documents and Settings\\jewell\\Application 
Data\\Tinn-R\\tmp\\selection.r"
> source(as.list(.trPaths)[5])
Error in source(as.list(.trPaths)[5]) : invalid connection

c) you are hard-coding the file paths, so if they change your TinnR will 
break.

d) I don't think the elements need to be named, but it probably does no 
harm.

If you want to hard-code your file paths, and want the elements named, I 
suggest this might work:
.trPaths <- c('C:/Users/dennis/AppData/Roaming/Tinn-R/tmp/',
  'C:/Users/dennis/AppData/Roaming/Tinn-R/tmp/search.txt',
  'C:/Users/dennis/AppData/Roaming/Tinn-R/tmp/objects.txt',
  'C:/Users/dennis/AppData/Roaming/Tinn-R/tmp/file.r',
  'C:/Users/dennis/AppData/Roaming/Tinn-R/tmp/selection.r',
  'C:/Users/dennis/AppData/Roaming/Tinn-R/tmp/block.r',
  'C:/Users/dennis/AppData/Roaming/Tinn-R/tmp/lines.r'
)
 names(.trPaths) <- c("Tmp", "Search", "Objects", "File", "Selection", 
"Block","Lines")

Hope that helps,

Keith J

"teck-corp"  wrote in message 
news:1268404061254-1590576.p...@n4.nabble.com...
>
> Hi Stephen,
>
> Thanks a lot for your answer. Unfortunately this does not work for me
> neither.
> Could you maybe let me know what is written in you RprofileSite-file now?
>
> Best
> Dennis

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

Re: [R] Find multiple elements in a vector

[Keith Jewell]

Try
which(x %in% y)

HTH

KJ

"Philip Twumasi-Ankrah"  wrote in message 
news:519822.56196...@web39502.mail.mud.yahoo.com...
This should work

which((x==2)|(x==3))

--Quotable Quotes-

A Smile costs Nothing But Rewards Everything
- Anonymous

Happiness is not perfected until it is shared
-Jane Porter


--- On Wed, 7/22/09, Michael Knudsen  wrote:

From: Michael Knudsen 
Subject: [R] Find multiple elements in a vector
To: "r help" 
Date: Wednesday, July 22, 2009, 2:32 PM

Hi,

Given a vector, say

x=sample(0:9,10)
x
[1] 0 6 3 5 1 9 7 4 8 2

I can find the location of an element by

which(x==2)
[1] 10

but what if I want to find the location of more than one number? I could do

c(which(x==2),which(x==3))

but isn't there something more streamlined? My first guess was

y=c(2,3)
which(x==y)
integer(0)

which doesn't work. I haven't found any clue in the R manual.

Thanks!

-- 
Michael Knudsen
micknud...@gmail.com
http://lifeofknudsen.blogspot.com/

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




[[alternative HTML version deleted]]







>

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

Re: [R] What does .[foo] really mean?

[Keith Jewell]

It's a place holder. See http://en.wikipedia.org/wiki/Foo

HTH

KJ

"Simon Tian"  wrote in message 
news:2119d1440908012150j373bfaf5v9d498bded94f...@mail.gmail.com...
> Hi R users,
>
> I really want to know what exactly ".[foo]" means.
>
> Thanks in advance. -Simon
>
> [[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] netlabR package in English

[Keith McMillan]

Dear R users,

 

Is documentation for the netlabR package available in English?

 

If not does anyone know if or when it will be?

 

Regards,

 

Keith

This
 e-mail and any files transmitted with it are confidential and are intended 
solely for the use of the individual or entity to whom it is addressed. If you 
are not the intended recipient or the person responsible for delivering the 
e-mail to the intended recipient, be advised that you have received this e-mail 
in error, and that any use, dissemination, forwarding, printing, or copying of 
this e-mail is strictly prohibited. If you received this e-mail in error, 
please return the e-mail to the sender at Merrick Bank and delete it from your 
computer. Although Merrick Bank attempts to sweep e-mail and attachments for 
viruses, it does not guarantee that either are virus-free and accepts no 
liability for any damage sustained as a result of viruses.

[[alternative HTML version deleted]]

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

Re: [R] how to update R files included as"source" ?

[Keith Jewell]

Hi,

I'm afraid I really don't have time to enter into a dialogue :-{ but this 
fragment from a function of mine might help...
--
  outfile <- sub("Rd$", "html", hfile, ignore.case=TRUE)  # name of 
corresponding html
  out.mod <- file.info(outfile)[,"mtime"]  # if html absent or needs 
updating
  if (is.na(out.mod)||file.info(hfile)[,"mtime"] > out.mod) 
tools::Rd2HTML(hfile, outfile, .Package) # do it
---
The third line (inter alia) checks the modification time of the 'target' 
file and decides whether to process it.

Hope that helps,

Keith J
---
"arnaud chozo"  wrote in message 
news:caebefd51003290824m67bebf96wd61315359ae97...@mail.gmail.com...
> Hi all,
>
> I have a main file main.R in which I include some other R files. For
> example, in main.R I have: source("functions.R").
> When I modify the file "functions.R", I'd like R to take into account the
> changes and to reload the file functions.R when I run main.R
> Is it possible?
>
> Thanks,
> Arnaud
>
> [[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] Kennard-stone or other uniform subset sampling algorithm?

[Woolner, Keith]

Is there any package in R which performs something like the
Kennard-Stone uniform mapping algorithm (a.k.a. uniform subset
selection)?

 

I've used Rseek.org to search the R-help archives and other resources,
and found only the following two mentions.

 

http://tolstoy.newcastle.edu.au/R/e2/help/07/01/.html

https://stat.ethz.ch/pipermail/r-help/2009-January/184856.html

 

Any pointers or suggestions of packages containing alternative
algorithms would be appreciated.

 

Thanks,

Keith

 


[[alternative HTML version deleted]]

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

[R] Problem with NLSstClosestX; and suggested fix

[Keith Jewell]

Problem is demonstrated with this code, intended to find the approximate 'x' 
at which the 'y' is midway between the left and right asymptotes. This 
particular data set returns NA, which is a bit silly!
--
sXY <- structure(list(x = c(0, 24, 27, 48, 51, 72, 75, 96, 99), y = 
c(4.98227,
6.38021, 6.90309, 7.77815, 7.64345, 7.23045, 7.27875, 7.11394,
6.95424)), .Names = c("x", "y"), row.names = c(NA, 9L), class = 
c("sortedXyData",
"data.frame"))
a <- NLSstLfAsymptote(sXY)
d <- NLSstRtAsymptote(sXY)
NLSstClosestX(sXY, (a+d)/2)

I think the problem arises when the target y value is exactly equal to one 
of the y values in sXY and can be fixed by trapping that situation thus:

NLSstClosestX.sortedXyData <-   function (xy, yval)
{
deviations <- xy$y - yval
if (any(deviations==0)) xy$x[match(0, deviations)] else {   # new line 
inserted
  if (any(deviations <= 0)) {
  dev1 <- max(deviations[deviations <= 0])
  lim1 <- xy$x[match(dev1, deviations)]
  if (all(deviations <= 0)) {
  return(lim1)
  }
  }
  if (any(deviations >= 0)) {
  dev2 <- min(deviations[deviations >= 0])
  lim2 <- xy$x[match(dev2, deviations)]
  if (all(deviations >= 0)) {
  return(lim2)
  }
  }
  dev1 <- abs(dev1)
  dev2 <- abs(dev2)
  lim1 + (lim2 - lim1) * dev1/(dev1 + dev2)
   }   # new line inserted
}
---

Comments/corrections welcome.

Keith Jewell

=

> sessionInfo()
R version 2.9.2 (2009-08-24)
i386-pc-mingw32

locale:
LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United 
Kingdom.1252;LC_MONETARY=English_United 
Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

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

other attached packages:
[1] nlme_3.1-93xlsReadWrite_1.3.3 svSocket_0.9-43svMisc_0.9-48 
TinnR_1.0.3R2HTML_1.59-1
[7] Hmisc_3.6-1

loaded via a namespace (and not attached):
[1] cluster_1.12.0  grid_2.9.2  lattice_0.17-25 stats4_2.9.2 
tools_2.9.2 VGAM_0.7-9

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 create a "heatline" -- heatmap in one dimension?

[Woolner, Keith]

Hi,

 

I'm trying to create a "heatmap" that is just a line segment.  That is,
a graphic that shows an interval and how a dependent variable varies
along that interval, with the value of the dependent variable shown by
color rather than on a y-axis.  

 

The image() function produces something close to what I'm envisioning,
but rather than plotting color on a 2-D grid, I need a line where the
color carries the information about the other dimension.  The following
example code gives a visual approximation of what I'm going for (using
image(), but by shrinking the y-axis down to make the result line-like):

 

x <- c(1:100)

y <- abs(40-x)

 

my.col <- heat.colors(10)

par(pin=c(par("pin")[1],.05))

image(x, 1, as.matrix(y), ylim=c(1,1), col=my.col, xaxt="n",xlab="",
yaxt="n", ylab="", fg="transparent")

 

 

I will be generating several of these "heatlines" to plot on a 2-D grid,
which will contain other information on it as well, and so also need to
be able to add the lines to an existing plot.

 

I've read the doc for stripplot(), contour(), image(), and checked the R
Graph Gallery (http://addictedtor.free.fr/graphiques/) to try to find
something similar, but haven't figured out the best way to do this,  I
envision a function similar to segments(), where you can specify where
the line should be drawn, the line width (lwd), and other graphical
parameters, but where you can specify a vector of color values.  

 

Is there a package that has such functionality already, and if not, what
would be most elegant way to approach it?  If it matters, I'm running R
2.9.1 on XP 32-bit.

 

Thanks for any suggestions.

 

--

Keith 

 

 

 

--

Keith Woolner

Manager, Baseball Research & Analysis

Cleveland Indians

216-420-4625

kwool...@indians.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] Vista truncated Sys.getenv("PATH") when and only when running as administrator

[Keith Ponting]

This is rather strange and I suspect is a Windows issue rather than an R
issue, but it does not

seem to be mentioned on the lists so I am logging it in case it gives
anyone else problems.

 

Running Vista Business Version 6.0.6002 Service Pack 2 Build 6002, with
the latest R 2.10.0 patched.

 

A normal 'R' gui process:

 

> p <- Sys.getenv("PATH")

> nchar(p)

PATH 

1199 

> sessionInfo()

R version 2.10.0 Patched (2009-11-23 r50548) 

i386-pc-mingw32 

 

locale:

[1] LC_COLLATE=English_United Kingdom.1252 

[2] LC_CTYPE=English_United Kingdom.1252   

[3] LC_MONETARY=English_United Kingdom.1252

[4] LC_NUMERIC=C   

[5] LC_TIME=English_United Kingdom.1252

 

attached base packages:

[1] stats graphics  grDevices utils datasets  methods   base


 

Running R gui as administrator, however, the same path variable is
truncated:

 

> p <- Sys.getenv("PATH")

> nchar(p)

PATH 

1024 

> sessionInfo()

R version 2.10.0 Patched (2009-11-23 r50548) 

i386-pc-mingw32 

 

locale:

[1] LC_COLLATE=English_United Kingdom.1252 

[2] LC_CTYPE=English_United Kingdom.1252   

[3] LC_MONETARY=English_United Kingdom.1252

[4] LC_NUMERIC=C   

[5] LC_TIME=English_United Kingdom.1252

 

attached base packages:

[1] stats graphics  grDevices utils datasets  methods   base


 

 

This may be related to
http://tolstoy.newcastle.edu.au/R/e6/devel/09/01/0259.html 

 

Keith Ponting

Aurix Ltd, Malvern WR14 3SZ  UK

 


[[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] multidimensional point.in.polygon??

[Keith Jewell]

Hi,

I seek to identify those points in/outside a multidimensional convex hull 
(geometry::convhulln). Any suggestions?

Background just in case I'm going down a really wrong road:

Given an observed data set with one dependent/observed variable (Y) and 
multiple (3 to 10) independent/design variables (X1, X2, ...) I want to 
increase the number of points by interpolating. I'm using expand.grid(X) to 
generate the X points and locfit::predict.locfit to interpolate the Y 
values. No problem so far.

BUT my observed X data set is very far from rectangular, so the set of 
points produced by expand.grid includes many  "extrapolations" which I'd 
want to remove. I'm aware of the problems of defining extrapolation in 
multiple dimensions and am minded to define it as "outside the convex hull", 
hence my question.

In searching r-project.org I came across a thread in which baptiste auguie 
said "one general way to do this would be to compute the convex hull 
(?chull) of the augmented set of points and test if the point belongs to 
it"; an approach I'd considered generalising to multiple points thus (pseudo 
R code)...

  baseHull <- convhulln(baseSet)
  augHull <- convhulln(augSet)
  while (augHull != baseHull)
{   augSet <- augSet [-(augHull !%in% baseHull)]
augHull <- convhulln(augSet)
}

... but this has that horrible loop including a convhulln!

In the (real, typical) test data set I'm using for development baseSet is 5 
columns by 2637 rows while baseHull has only 42 distinct points. If augHull 
has a similarly simple hull, then each time round the loop will remove only 
a few dozen points; I need to remove many thousands.

And (in my naivete) I think there must be a better way of testing for a 
point inside a polygon, I'd have thought convhulln must need to do that 
often!

Thanks in advance for any pointers,

Keith Jewell

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

Re: [R] multidimensional point.in.polygon??

[Keith Jewell]

Hi,

Doing some more reading, I think the problem is easier because the hull is 
convex. Then an algorithm for testing points might be:

a) Define the convex hull as a set of planes (simplexes).
[as returned by convhulln!!]

b) Define one point, i, known to be interior
[e.g. mean of all the points defining the hull]

c) If point x is
i) for every plane, on the same side as i; x is interior
   ii) for every plane, on the same side as i or in the plane; x is in the 
surface
 iii) else x is exterior

So now I need to find the directions of points from multidimensional 
planes.Perhaps I can find the vectors of the perpendiculars from the points 
to the planes (possibly extended) and test for parallel/anti-parallel?

I feel that I'm in the right direction because this uses the structure of a 
convex hull returned by convhulln. But, I still feel I'm re-inventing the 
wheel. Surely this has been done before? Isn't a (the?) major purpose of a 
convex hull to test other points for inclusion?

Perhaps when I get the geometry sorted this will be so easy I'll understand 
why noone has pointed me to an existing R function, but currently I feel I 
and Baptiste are wandering in the dark :-(

Any hints?

Thanks in advance,

Keith Jewell
-
"baptiste auguie"  wrote in message 
news:de4e29f50912040550m71fbffafnfa1ed6e0f4451...@mail.gmail.com...
Hi,

Yet another one of my very naive ideas on the subject: maybe you can
first evaluate the circumscribed and inscribed spheres of the base set
of points (maximum and minimum of their distances to the center of
gravity). Any points within a distance smaller than the infimum is
good, any point further than the supremum is not good. This should be
faster than the calculation of a convex hull for each point. Of course
the usefulness of this first test really depends on how aspherical is
your base convex hull.

I do hope to read a real answer from someone who knows this stuff!

HTH,

baptiste


2009/12/4 Keith Jewell :
> Hi,
>
> I seek to identify those points in/outside a multidimensional convex hull
> (geometry::convhulln). Any suggestions?
>
> Background just in case I'm going down a really wrong road:
>
> Given an observed data set with one dependent/observed variable (Y) and
> multiple (3 to 10) independent/design variables (X1, X2, ...) I want to
> increase the number of points by interpolating. I'm using expand.grid(X) 
> to
> generate the X points and locfit::predict.locfit to interpolate the Y
> values. No problem so far.
>
> BUT my observed X data set is very far from rectangular, so the set of
> points produced by expand.grid includes many "extrapolations" which I'd
> want to remove. I'm aware of the problems of defining extrapolation in
> multiple dimensions and am minded to define it as "outside the convex 
> hull",
> hence my question.
>
> In searching r-project.org I came across a thread in which baptiste auguie
> said "one general way to do this would be to compute the convex hull
> (?chull) of the augmented set of points and test if the point belongs to
> it"; an approach I'd considered generalising to multiple points thus 
> (pseudo
> R code)...
> 
> baseHull <- convhulln(baseSet)
> augHull <- convhulln(augSet)
> while (augHull != baseHull)
> { augSet <- augSet [-(augHull !%in% baseHull)]
> augHull <- convhulln(augSet)
> }
> 
> ... but this has that horrible loop including a convhulln!
>
> In the (real, typical) test data set I'm using for development baseSet is 
> 5
> columns by 2637 rows while baseHull has only 42 distinct points. If 
> augHull
> has a similarly simple hull, then each time round the loop will remove 
> only
> a few dozen points; I need to remove many thousands.
>
> And (in my naivete) I think there must be a better way of testing for a
> point inside a polygon, I'd have thought convhulln must need to do that
> often!
>
> Thanks in advance for any pointers,
>
> Keith Jewell
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] global variables in a function

[Keith Jones]

Y'all,

I would like to have most of the variables in my function to be global  
instead of local.  I have not found any references that tell me now to  
do that.  If I have missed a reference please excuse me and tell me  
what I have missed.


Thanks,

Keith Jones

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] mutlidimensional in.convex.hull (was multidimensional point.in.polygon??)

[Keith Jewell]

Hi All (especially Duncan and Baptiste),

Summary (of lengthy bits below):
I will have a convex hull in multiple (3 to 10) dimensions derived from many 
(thousands) of points by geometry::convhulln.
I will need to categorise other 'test' points as inside/outside that convex 
hull . e.g. given:
--
require(geometry)
ps <- matrix(rnorm(4000),ncol=4) # 'calibration' set
phull <-  convhulln(ps)  # convex hull
xs <- matrix(rnorm(1200),ncol=4)# 'test' set
-
How do I categorise each point (row) in xs as inside/outside(/on) phull???
There is tripack::in.convex.hull but that doesn't handle my dimensionality.

Thanks to Duncan Murdoch for the suggestion (just a few lines down, 
previously made by Baptiste Auguie): of testing a single point thus:
  i) add the test point to the set of points defining the convex hull,
  ii) recalculate the convex hull
  iii) if the test point is part of the new convex hull, then it was outside 
the original

BUT I have many (thousands) of test points, so this would involve very many 
convex hull calculations. My suggestion, immediately below, requires finding 
the signs of perpendicular distances from each test point to each 
multidimensional 'plane' defining the convex hull (NB: phull  is a matrix in 
which each row defines such a 'plane').

Baptiste has found a Matlab implementation 
<http://www.mathworks.com/matlabcentral/fileexchange/10226-inhull> of (what 
looks like) my algorithm. I don't speak Matlab, but this looks non-trivial 
to code in R. I'll do it if I have to, but if it already exists it would be 
nice. If I do have to code it, I'd really appreciate an expression in 
algebra rather than Matlab!

Any pointers will be much appreciated,

Keith Jewell
"Duncan Murdoch"  wrote in message 
news:4b20e1ea.3030...@stats.uwo.ca...
> On 10/12/2009 5:15 AM, Keith Jewell wrote:
>> Hi,
>>
>> Doing some more reading, I think the problem is easier because the hull 
>> is convex. Then an algorithm for testing points might be:
>>
>> a) Define the convex hull as a set of planes (simplexes).
>> [as returned by convhulln!!]
>>
>> b) Define one point, i, known to be interior
>> [e.g. mean of all the points defining the hull]
>>
>> c) If point x is
>> i) for every plane, on the same side as i; x is interior
>>ii) for every plane, on the same side as i or in the plane; x is in 
>> the surface
>>  iii) else x is exterior
>
> That looks like it would work, but wouldn't it be easier to do the 
> following:
>
> Compute the convex hull with the new point added. If the point is 
> exterior, the new point will be part of the hull.  If interior, it won't 
> be.  If it is on the boundary, it's probably unpredictable, but due to 
> rounding error, that's probably true even with a perfect algorithm.
>
> I didn't notice that you said how your original polygon is defined, but if 
> it is defined as a convex hull or in terms of its vertices, the above 
> method would work.  If it's defined some other way, it might be hard.
>
> Duncan Murdoch
>
>
>>
>> So now I need to find the directions of points from multidimensional 
>> planes.Perhaps I can find the vectors of the perpendiculars from the 
>> points to the planes (possibly extended) and test for 
>> parallel/anti-parallel?
>>
>> I feel that I'm in the right direction because this uses the structure of 
>> a convex hull returned by convhulln. But, I still feel I'm re-inventing 
>> the wheel. Surely this has been done before? Isn't a (the?) major purpose 
>> of a convex hull to test other points for inclusion?
>>
>> Perhaps when I get the geometry sorted this will be so easy I'll 
>> understand why noone has pointed me to an existing R function, but 
>> currently I feel I and Baptiste are wandering in the dark :-(
>>
>> Any hints?
>>
>> Thanks in advance,
>>
>> Keith Jewell
>> -
>> "baptiste auguie"  wrote in message 
>> news:de4e29f50912040550m71fbffafnfa1ed6e0f4451...@mail.gmail.com...
>> Hi,
>>
>> Yet another one of my very naive ideas on the subject: maybe you can
>> first evaluate the circumscribed and inscribed spheres of the base set
>> of points (maximum and minimum of their distances to the center of
>> gravity). Any points within a distance smaller than the infimum is
>> good, any point further than the supremum is not good. This should be
>> faster than the calculation of a convex hull for each point. Of course
>> the 

Re: [R] mutlidimensional in.convex.hull (was multidimensionalpoint.in.polygon??)

[Keith Jewell]

S
  n missing  uniqueMean
  11200   0   4   3
0 (2800, 25%), 2 (2800, 25%), 4 (2800, 25%), 6 (2800, 25%)
-
N
  n missing  uniqueMean
  11200   0   4 120
0 (2800, 25%), 80 (2800, 25%), 160 (2800, 25%), 240 (2800, 25%)
-
> cor(xs)
  tpH T S N
t  1.00e+00  4.276263e-23  2.346103e-20 0  0.00e+00
pH 4.276263e-23  1.00e+00 -1.716171e-19 0  0.00e+00
T  2.346103e-20 -1.716171e-19  1.00e+00 0 -3.925415e-21
S  0.00e+00  0.00e+00  0.00e+00 1  0.00e+00
N  0.00e+00  0.00e+00 -3.925415e-21 0  1.00e+00
> phull <-  convhulln(ps)
> phull2 <- convhulln(rbind(ps,xs))
> nrp <- nrow(ps)
> nrx <- nrow(xs)
> outside <- unique(phull2[phull2>nrp])-nrp
> done <- FALSE
> print(length(outside))
[1] 20
> while(!done){
phull3 <- convhulln(rbind(ps,xs[-(outside),]))
also.outside <- (1:nrx)[-outside][unique(phull3[phull3>nrp])-nrp]
outside <- c(outside,also.outside)
print(length(outside))
done <- length(also.outside)==0
}
[1] 97
[1] 243
[1] 395
[1] 501
[1] 571
[1] 642
[1] 718
[1] 790
[1] 836
[1] 876
[1] 927
[1] 973
[1] 1012
[1] 1038
[1] 1060
[1] 1088
[1] 1120
[1] 1148
[1] 1162
[1] 1168
[1] 1170
[1] 1170
> describe(xs[-(outside),])
xs[-(outside), ]
 5  Variables  10030  Observations
-
t
  n missing  uniqueMean .05 .10 .25 .50 .75 
.90 .95
  10030   0  35   120.5   2   4  30 120 172 
240 336
lowest :   0   2   3   4  24, highest: 312 336 360 384 504
-
pH
  n missing  uniqueMean
  10030   0   4   5.781
4.6 (2548, 25%), 5.4 (2559, 26%), 6.2 (2520, 25%), 7 (2403, 24%)
-
T
  n missing  uniqueMean
  10030   0   5   9.584

 258   15   22
Frequency 2099 2232 2174 2024 1501
%   21   22   22   20   15
-
S
  n missing  uniqueMean
  10030   0   4   3.067
0 (2371, 24%), 2 (2512, 25%), 4 (2572, 26%), 6 (2575, 26%)
-
N
  n missing  uniqueMean
  10030   0   4   120.0
0 (2502, 25%), 80 (2514, 25%), 160 (2515, 25%), 240 (2499, 25%)
-
> cor(xs[-(outside),])
   tpH T S N
t   1.00 -0.0252190240 -0.1599673922  0.0260853161 -0.0001862074
pH -0.0252190240  1.00 -0.0308956477 -0.0109761096  0.0006346428
T  -0.1599673922 -0.0308956477  1.00  0.0506125967  0.0005164211
S   0.0260853161 -0.0109761096  0.0506125967  1.00  0.0006132354
N  -0.0001862074  0.0006346428  0.0005164211  0.0006132354  1.00
### 

Baptiste said:
>
>Regarding the "inhull" Matlab code, I came to the opposite conclusion:
>it should be easily ported to R. 1) it is a very short piece of code
>(even more so if one disregards the various checks and handling of
>special cases), with no Matlab-specific objects (only integers,
>booleans, matrices and vectors). 2) The core of the program relies on
>the qhull library, and the same applies to R I think. 3) Matlab and R
>use very similar indexing for matrices and similar linear algebra in
>general.
>
>That said, I'm a bit short on time to give it a go myself. I think the
>open-source Octave could run this code too, so it might help in
>checking the code step-by-step.
>
>
>All the best,
>
>baptiste

I agree the Matlab code didn't look very complicated. The main reason I 
thought it non-trival (for me) to code in R was my ignorance; lines like...
nullsp = null(xyz(tess(i,2:end),:) - repmat(xyz(tess(i,1),:),p-1,1))';
...seem important but I don't know what 'null' and 'repmat' do. I don't even 
understand the algebra being implemented and I'm wary of coding in those 
circumstances.

If I had time, I'd lik

Re: [R] mutlidimensional in.convex.hull (wasmultidimensionalpoint.in.polygon??)

[Keith Jewell]

Hi All,

I couldn't resist doing this the "right" way!
A colleague explained the vector algebra to me (thanks Martin!) and I've 
followed the structure of the Matlab code referenced below, but all errors 
are mine!

I don't pretend to great R expertise, so this code may not be optimal (in 
time or the memory issue addressed by the Matlab code), but it is a lot 
faster than the other algorithm discussed in this thread. These timings are 
on the real example data described in my previous message in this thread.

> system.time(inhull(xs,ps))  # the "right" way
   user  system elapsed
   1.340.071.41
> system.time({phull <-  convhulln(ps) # the other algorithm
+ phull2 <- convhulln(rbind(ps,xs))
+ nrp <- nrow(ps)
+ nrx <- nrow(xs)
+ outside <- unique(phull2[phull2>nrp])-nrp
+ done <- FALSE
+ while(!done){
+ phull3 <- convhulln(rbind(ps,xs[-(outside),]))
+ also.outside <- (1:nrx)[-outside][unique(phull3[phull3>nrp])-nrp]
+ outside <- c(outside,also.outside)
+ done <- length(also.outside)==0
+ }})
   user  system elapsed
  15.090.09   15.22

Although I really must move on now, if anyone has comments, criticisms, 
suggestions for improvement I would be interested.

Enjoy!

Keith Jewell

inhull <- function(testpts, calpts, hull=convhulln(calpts), 
tol=mean(mean(abs(calpts)))*sqrt(.Machine$double.eps)) {
#
# R implementation of the Matlab code by John D'Errico 04 Mar 2006 (Updated 
30 Oct 2006)
# downloaded from 
http://www.mathworks.com/matlabcentral/fileexchange/10226-inhull
# with some modifications and simplifications
#
# Efficient test for points inside a convex hull in n dimensions
#
#% arguments: (input)
#%  testpts - nxp array to test, n data points, in p dimensions
#%   If you have many points to test, it is most efficient to
#%   call this function once with the entire set.
#%
#%  calpts - mxp array of vertices of the convex hull, as used by
#%   convhulln.
#%
#%  hull - (OPTIONAL) tessellation (or triangulation) generated by convhulln
#%   If hull is left empty or not supplied, then it will be
#%   generated.
#%
#%  tol - (OPTIONAL) tolerance on the tests for inclusion in the
#%   convex hull. You can think of tol as the distance a point
#%   may possibly lie outside the hull, and still be perceived
#%   as on the surface of the hull. Because of numerical slop
#%   nothing can ever be done exactly here. I might guess a
#%   semi-intelligent value of tol to be
#%
#% tol = 1.e-13*mean(abs(calpts(:)))
#%
#%   In higher dimensions, the numerical issues of floating
#%   point arithmetic will probably suggest a larger value
#%   of tol.
#%
# In this R implementation default 
tol=mean(mean(abs(calpts)))*sqrt(.Machine$double.eps)
#   DEFAULT: tol = 1e-6
#
# VALUE: Matlab returns a vector of TRUE (inside/on) or FALSE (outside)
#   This R implementation returns an integer vector of length n
#   1 = inside hull
#  -1 = inside hull
#   0 = on hull (to precision indicated by tol)
#
   require(geometry, quietly=TRUE)  # for  convhulln
   require(MASS, quietly=TRUE)  # for Null
# ensure arguments are matrices (not data frames) and get sizes
   calpts <- as.matrix(calpts)
   testpts <- as.matrix(testpts)
   p <- dim(calpts)[2]   # columns in calpts
   cx <- dim(testpts)[1]  # rows in testpts
   nt <- dim(hull)[1]# number of simplexes in hull
# find normal vectors to each simplex
   nrmls <- matrix(NA, nt, p) # predefine each nrml as NA, 
degenerate
   degenflag <- matrix(TRUE, nt, 1)
   for (i in  1:nt) {
nullsp <- t(Null(t(calpts[hull[i,-1],] - 
matrix(calpts[hull[i,1],],p-1,p, byrow=TRUE
if (dim(nullsp)[1] == 1) { nrmls[i,] <- nullsp
   degenflag[i] <- FALSE}}
# Warn of degenerate faces, and remove corresponding normals
   if(length(degenflag[degenflag]) > 0) 
warning(length(degenflag[degenflag])," degenerate faces in convex hull")
   nrmls <- nrmls[!degenflag,]
   nt <- dim(nrmls)[1]
# find center point in hull, and any (1st) point in the plane of each 
simplex
   center = apply(calpts, 2, mean)
   a <- calpts[hull[!degenflag,1],]
# scale normal vectors to unit length and ensure pointing inwards
   nrmls <- nrmls/matrix(apply(nrmls, 1, function(x) sqrt(sum(x^2))), nt, p)
   dp <- sign(apply((matrix(center, nt, p, byrow=TRUE)-a) * nrmls, 1, sum))
   nrmls <- nrmls*matrix(dp, nt, p)
# if  min across all faces of dot((x - a),nrml) is
#  +ve then x is inside hull
#  0   then x is on hull
#  -ve then x is outside hull
# Instead of dot((x - a),nrml)  use dot(x,nrml) - dot(a, nrml)
   aN <- diag(a %*% t(nrmls))
   val <- apply(testpts %*% t(nrmls) - matrix(aN, cx, nt, byrow=TRUE), 
1,min)
# code  values inside 'tol' to zer

Re: [R] mutlidimensional in.convex.hull(wasmultidimensionalpoint.in.polygon??)

[Keith Jewell]

Syntax suggestions implemented.
Inhull's original author consulted.
Function submitted for potential inclusion in geometry package.

Seasons greetings to all.

Keith Jewell
-
"baptiste auguie"  wrote in message 
news:de4e29f50912180238m45c4b7c3g792113d5b6c4c...@mail.gmail.com...
Hi,

Excellent, thanks for doing this!

I had tried the 2D case myself but I was put off by the fact that
Octave's convhulln had a different ordering of the points to R's
geometry package (an improvement to Octave's convhulln was made after
it was ported to R). I'm not sure how you got around this but it's
good to see that the result was worth the effort!

Below are a few minor syntax suggestions,

# p <- dim(calpts)[2]   # columns in calpts
# and other similar lines could be replaced with
ncol(calpts)
nrow(testpts)
nrow(hull)

# length(degenflag[degenflag])
# can probably be written
sum(degenflag)

# center = apply(calpts, 2, mean)
# more efficient
colMeans(calpts)

Would you consider submitting this function to the maintainer of the
geometry package, after checking it's OK with inhull's original
author?

Best regards,

baptiste

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

[R] Problem with expand.grid

[Keith Jewell]

Hi All,

This example code

dDF <- structure(list(y = c(4.75587, 4.8451, 5.04139, 4.85733, 5.20412,
 5.92428, 5.69897, 4.78958, 4, 4), t = c(0, 48, 144, 192, 240,
 312, 360, 0, 48, 144), Batch = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1
 ), T = c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2), pH = c(4.6, 4.6, 4.6,
 4.6, 4.6, 4.6, 4.6, 4.6, 4.6, 4.6), S = c(0, 0, 0, 0, 0, 0, 0,
 0, 0, 0), N = c(0, 0, 0, 0, 0, 0, 0, 80, 80, 80)), .Names = c("y",
 "t", "Batch", "T", "pH", "S", "N"), row.names = c(NA, 10L), class = 
"data.frame")
str(dDF)
expand.grid(dDF)

'hangs' for a while and then gives an error

Error in `[[<-.data.frame`(`*tmp*`, i, value = c(4.75587, 4.8451, 5.04139, 
:
  replacement has 1000 rows, data has 10

In NEWS.R-2.11.0dev I read:
o   The new (in 2.9.0) 'stringsAsFactors' argument to expand.grid()
was not working: it now does work but has default TRUE for
backwards compatibility.

but I don't think that's relevant, I have no factors.

I'm probably being silly. Can anyone point out where?

Best...

Keith Jewell

--please do not edit the information below--

Version:
 platform = i386-pc-mingw32
 arch = i386
 os = mingw32
 system = i386, mingw32
 status = Patched
 major = 2
 minor = 10.1
 year = 2009
 month = 12
 day = 21
 svn rev = 50796
 language = R
 version.string = R version 2.10.1 Patched (2009-12-21 r50796)

Windows Server 2003 x64 (build 3790) Service Pack 2

Locale:
LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United 
Kingdom.1252;LC_MONETARY=English_United 
Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

Search Path:
 .GlobalEnv, package:stats, package:graphics, package:grDevices, 
package:utils, package:datasets, package:methods, Autoloads, package:base

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

Re: [R] Problem with expand.grid

[Keith Jewell]

Just confirming it isn't the bug fixed in 2.11.0dev, and giving an even 
simpler example:

R version 2.11.0 Under development (unstable) (2009-12-20 r50794)

> expand.grid(data.frame(y=1:10, t=1:10))
Error in `[[<-.data.frame`(`*tmp*`, i, value = c(1L, 2L, 3L, 4L, 5L, 6L,  :
  replacement has 100 rows, data has 10

"Keith Jewell"  wrote in message 
news:hgqqja$rk...@ger.gmane.org...
> Hi All,
>
> This example code
> 
> dDF <- structure(list(y = c(4.75587, 4.8451, 5.04139, 4.85733, 5.20412,
> 5.92428, 5.69897, 4.78958, 4, 4), t = c(0, 48, 144, 192, 240,
> 312, 360, 0, 48, 144), Batch = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1
> ), T = c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2), pH = c(4.6, 4.6, 4.6,
> 4.6, 4.6, 4.6, 4.6, 4.6, 4.6, 4.6), S = c(0, 0, 0, 0, 0, 0, 0,
> 0, 0, 0), N = c(0, 0, 0, 0, 0, 0, 0, 80, 80, 80)), .Names = c("y",
> "t", "Batch", "T", "pH", "S", "N"), row.names = c(NA, 10L), class = 
> "data.frame")
> str(dDF)
> expand.grid(dDF)
> 
> 'hangs' for a while and then gives an error
>
> Error in `[[<-.data.frame`(`*tmp*`, i, value = c(4.75587, 4.8451, 5.04139, 
> :
>  replacement has 1000 rows, data has 10
>
> In NEWS.R-2.11.0dev I read:
>o The new (in 2.9.0) 'stringsAsFactors' argument to expand.grid()
> was not working: it now does work but has default TRUE for
> backwards compatibility.
>
> but I don't think that's relevant, I have no factors.
>
> I'm probably being silly. Can anyone point out where?
>
> Best...
>
> Keith Jewell
>
> --please do not edit the information below--
>
> Version:
> platform = i386-pc-mingw32
> arch = i386
> os = mingw32
> system = i386, mingw32
> status = Patched
> major = 2
> minor = 10.1
> year = 2009
> month = 12
> day = 21
> svn rev = 50796
> language = R
> version.string = R version 2.10.1 Patched (2009-12-21 r50796)
>
> Windows Server 2003 x64 (build 3790) Service Pack 2
>
> Locale:
> LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United 
> Kingdom.1252;LC_MONETARY=English_United 
> Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252
>
> Search Path:
> .GlobalEnv, package:stats, package:graphics, package:grDevices, 
> package:utils, package:datasets, package:methods, Autoloads, package:base
>

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

Re: [R] Problem with expand.grid

[Keith Jewell]

Thanks, Prof Ripley (and everyone else who's responding while I type this!).

I thought I was being stupid. What I wanted was:
  expand.grid(lapply(dDF, unique))
which works fine!



Seasonal greetings to all,

Keith Jewell

"Prof Brian Ripley"  wrote in message 
news:alpine.lfd.2.00.0912221619440.20...@gannet.stats.ox.ac.uk...
> 1) Look at the help: a data frame may be a list, but do pass a list such 
> as unclass(dDF) when it says 'list'.
>
> 2) You have 7 columns of 10 items, which gives 10 million rows.  Is that 
> really what you want (especially as some of the columns are constant)? 
> It's an object of ca 500MB.
>
> On Tue, 22 Dec 2009, Keith Jewell wrote:
>
>> Hi All,
>>
>> This example code
>> 
>> dDF <- structure(list(y = c(4.75587, 4.8451, 5.04139, 4.85733, 5.20412,
>> 5.92428, 5.69897, 4.78958, 4, 4), t = c(0, 48, 144, 192, 240,
>> 312, 360, 0, 48, 144), Batch = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1
>> ), T = c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2), pH = c(4.6, 4.6, 4.6,
>> 4.6, 4.6, 4.6, 4.6, 4.6, 4.6, 4.6), S = c(0, 0, 0, 0, 0, 0, 0,
>> 0, 0, 0), N = c(0, 0, 0, 0, 0, 0, 0, 80, 80, 80)), .Names = c("y",
>> "t", "Batch", "T", "pH", "S", "N"), row.names = c(NA, 10L), class =
>> "data.frame")
>> str(dDF)
>> expand.grid(dDF)
>> 
>> 'hangs' for a while and then gives an error
>>
>> Error in `[[<-.data.frame`(`*tmp*`, i, value = c(4.75587, 4.8451, 
>> 5.04139,
>> :
>>  replacement has 1000 rows, data has 10
>>
>> In NEWS.R-2.11.0dev I read:
>>    o The new (in 2.9.0) 'stringsAsFactors' argument to expand.grid()
>> was not working: it now does work but has default TRUE for
>> backwards compatibility.
>>
>> but I don't think that's relevant, I have no factors.
>>
>> I'm probably being silly. Can anyone point out where?
>>
>> Best...
>>
>> Keith Jewell
>>
>> --please do not edit the information below--
>>
>> Version:
>> platform = i386-pc-mingw32
>> arch = i386
>> os = mingw32
>> system = i386, mingw32
>> status = Patched
>> major = 2
>> minor = 10.1
>> year = 2009
>> month = 12
>> day = 21
>> svn rev = 50796
>> language = R
>> version.string = R version 2.10.1 Patched (2009-12-21 r50796)
>>
>> Windows Server 2003 x64 (build 3790) Service Pack 2
>>
>> Locale:
>> LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United
>> Kingdom.1252;LC_MONETARY=English_United
>> Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252
>>
>> Search Path:
>> .GlobalEnv, package:stats, package:graphics, package:grDevices,
>> package:utils, package:datasets, package:methods, Autoloads, package:base
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide 
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
> -- 
> Brian D. Ripley,  rip...@stats.ox.ac.uk
> Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
> University of Oxford, Tel:  +44 1865 272861 (self)
> 1 South Parks Road, +44 1865 272866 (PA)
> Oxford OX1 3TG, UKFax:  +44 1865 272595
>

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

Re: [R] Wilcoxon nonparametric p-values

[Keith Jewell]

I just tried it in Minitab and got
--
Test of median = 8.500 versus median not = 8.500

N for   Wilcoxon Estimated
 N   Test  Statistic  P Median
C1  20 19   80.5  0.573  8.460
-
One tailed gave me closer to the textbook, but still not very close.
---
Test of median = 8.500 versus median < 8.500

N for   Wilcoxon Estimated
 N   Test  Statistic  P Median
C1  20 19   80.5  0.287  8.460
---

I agree with Peter Dalgaard, the book has it wrong (for some value of 
"wrong")

Regards

KJ

"Peter Dalgaard"  wrote in message 
news:4a127d6d.1060...@biostat.ku.dk...
> cvandy wrote:
>> When I use wilcox.test, I get vastly different p-values than the problems
>> from Statistics textbooks.
>> For example:
>> The following problem comes from "Applied Statistics and Probability for
>> Engineers", 2nd Edition, by D. C. Montgomery.  Page736, problem 14.7. 
>> The
>> problem is to compare the sample data with a population median of 8.5. 
>> The
>> book answer is p = 0.25, wilcox.test answer is p = 0.573.
>> I've tried several other similar problems with similar results.  I've 
>> copied
>> the following directly from my workspace.
>
> wilcox.exact (from exactRankTests) gives
>
> > wilcox.exact(x - 8.5)
>
> Exact Wilcoxon signed rank test
>
> data:  x - 8.5
> V = 80.5, p-value = 0.5748
>
> so I'd suspect the textbook. One-sided p-value perhaps? or table 
> limitation (as in "p > .25"). If you want to dig deeper, you'll probably 
> have to check the computations implied by the text.
>
>> Thanks for any help,
>> CHV
>>> x<-c(8.32,8.05,
>>> 8.93,8.65,8.25,8.46,8.52,8.35,8.36,8.41,8.42,8.30,8.71,8.75,8.6,8.83,8.5,8.38,8.29,8.46)
>>> wilcox.test(x,y=NULL,mu=8.5)
>> Wilcoxon signed rank test with continuity correction
>>  data:  x V = 80.5, p-value = 0.573
>> alternative hypothesis: true location is not equal to 8.5 Warning 
>> messages:
>> 1: In wilcox.test.default(x, y = NULL, mu = 8.5) :
>>   cannot compute exact p-value with ties
>> 2: In wilcox.test.default(x, y = NULL, mu = 8.5) :
>>   cannot compute exact p-value with zeroes
>> ? ? Charles H Van deZande
>>  ?
>
>
> -- 
>O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
>   c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
>  (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45) 35327918
> ~~ - (p.dalga...@biostat.ku.dk)  FAX: (+45) 35327907
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] subset POSIXct

[Keith Jones]

Hi,

I have a data frame with two columns: dt and tf.  The dt column is  
datetime and the tf column is a temperature.


   dt  tf
1 2009-06-20 00:53:00 73
2 2009-06-20 01:08:00 73
3 2009-06-20 01:44:00 72
4 2009-06-20 01:53:00 71
5 2009-06-20 02:07:00 72
...

I need a subset of the rows where the minutes are 53.  The hour is  
immaterial.  I can not find a wildcard character to use to select the  
rows.


as.character(wtd$dt) %in% "2009-06-21 ??:53:00" or as.character(wtd 
$dt) %in% "2009-06-21 **:53:00" does not work.


What would you recommend?

Thanks,

Keith Jones, Ph.D.
VTS Pumps

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

Re: [R] Factor to continuous

[Keith Jewell]

I suspect the OP is looking for FAQ 7.10

http://cran.r-project.org/doc/FAQ/R-FAQ.html#How-do-I-convert-factors-to-numeric_003f

HTH

Keith J

"Luc Villandre"  wrote in message 
news:49e495fe.1090...@dms.umontreal.ca...
> Hi,
>
> You'll need to be more specific about the nature of your problem. Could 
> you describe what exactly you're trying to do and why you need to 
> transform factors into continuous variables? Could you also provide a 
> small subsample of the data you're working with? Without such information, 
> I doubt that even the most knowledgeable specialists of R might be able to 
> provide meaningful assistance.
>
> Best of luck,
>
> Luc
>
> joewest wrote:
>> Hi
>>
>> I am really struggling with changing my factors into continuous 
>> variables. There is plenty of information on changing continuous to a 
>> factor which is
>> the opposite of what i need to do.  I would be soo grateful for any 
>> help
>>
>> Thanks
>>
>> Joe
>>
>

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

Re: [R] find the corresponding mean y

[Keith Jewell]

Not very clear what you want, but perhaps...
?sortedXyData
...might help.

hth

KJ
"luciferyan"  wrote in message 
news:1263231740556-1011427.p...@n4.nabble.com...
>
> Hello, I have 49 paired data, x, y.
> I have sampled x (where replacement is true), and find its mean.
> How can I find the corresponding mean y, which is the paired data of above
> sample x?
> Thank you very much,
> Annie
> -- 
> View this message in context: 
> http://n4.nabble.com/find-the-corresponding-mean-y-tp1011427p1011427.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] update.packages on MS Windows with //server/share paths

[Keith Jewell]

Hi,

> update.packages(ask='graphics')

  gives me  multiple warning (one per updated package?) similar to ...

Warning: unable to move temporary installation 
'\\Server02\stats\R\library\2.10\file3de56e0d\locfit' to 
'\\Server02\stats\R\library\2.10\locfit'

The final, updated, folders do not end up where they should be. I can move 
them 'by hand', but it is an inconvenience.

Checking the archives I find 
http://finzi.psych.upenn.edu/Rhelp08/2008-April/160963.html where Prof 
Ripley opines "The issue appears to be that your OS is garbling file names", 
and it surely does seem to be a filename problem - in combination R and 
Windows don't seem to be handling the '//server/share' path construction. I 
have:
> .libPaths()
[1] "//Server02/stats/R/library/2.10" 
"//Server02/stats/R/R-Current/library"

I can work around by mapping a drive letter to the '//server/share', but 
prefer not to (for local reasons).

In CHANGES.R-2.10.1pat I find
  CHANGES IN R VERSION 2.7.2 patched
o dir.create(recursive = TRUE) was not working on //server/share  paths.

In CHANGES.R-2.11.0dev I find
  CHANGES IN R VERSION 2.11.0
NEW FEATURES
o file.rename() can work across volumes (by copy-and-delete)

In another R function I've hit a similar problem and solved it using the 
construction:
  file.path(dirname(x), basename(x))
to convert  '//server/share' to (printed as) 'server/share' which worked 
in that function. In this case:
> file.path(dirname(.libPaths()), basename(.libPaths()))
[1] "Server02/stats/R/library/2.10" 
"Server02/stats/R/R-Current/library"
which looks OK but
> .libPaths(file.path(dirname(.libPaths()), basename(.libPaths(
> .libPaths()
[1] "//Server02/stats/R/library/2.10" 
"//Server02/stats/R/R-Current/library"
doesn't actually change the internal paths.

I don't see the problem in V2.9.2 .
I do see the problem in
   V2.10.1,
   V2.10.1 Patched (2010-01-24 r51030)
   V2.11.0 Under development (unstable) (2010-01-24 r51030)

Grateful for any suggestions,

Keith Jewell
-
Version:
 platform = i386-pc-mingw32
 arch = i386
 os = mingw32
 system = i386, mingw32
 status = Patched
 major = 2
 minor = 10.1
 year = 2010
 month = 01
 day = 24
 svn rev = 51030
 language = R
 version.string = R version 2.10.1 Patched (2010-01-24 r51030)

Windows Server 2003 x64 (build 3790) Service Pack 2

Locale:
LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United 
Kingdom.1252;LC_MONETARY=English_United 
Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

Search Path:
 .GlobalEnv, package:stats, package:graphics, package:grDevices, 
package:utils, package:datasets, package:methods, Autoloads, package:base

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

Re: [R] update.packages on MS Windows with //server/share paths

[Keith Jewell]

Thanks Gabor, but expanding on my remark...

> I can work around by mapping a drive letter to the '//server/share', but
> prefer not to (for local reasons).

R is installed on our network, used by others (with readonly access), 
maintained by me.
In 'Renviron.site' I have...
   R_LIBS_SITE=//Server02/stats/R/library/%v
...which works for everyone else (and for me as long as I don't try 
update.packages).

I (but not everyone else!) routinely have a drive (L:) mapped to 
//Server02/stats so that I can (and currently do) update.packages by
a) altering 'Renviron.site' to R_LIBS_SITE=L://R/library/%v
b) starting Rgui.exe from (e.g.)  L:\R\R-2.10.1\bin
BUT while 'Renviron.site' is altered, users who don't have the drive mapping 
don't get the site library, so I have to be quick and risk a couple of 
complaints.

I've tried putting a '.Renviron' file in my personal home directory 
containing R_LIBS_SITE=L://R/library/%v hoping I could implement a 
'personal' site library, but it doesn't seem to have any effect (not 
unreasonably?). As I outlined, my attempts to alter .libPaths() or 
.Library.site within my R session haven't worked.

Anyway, (although I say this with trepidation) this does look to me like a 
minor bug. It seems to me that (for example)...
> shell.exec(.libPaths()[1])
  and
> shell.exec(file.path(dirname(.libPaths()[1]), basename(.libPaths()[1])))
...should give the same result, especially when (as in this case) the syntax 
of the file spec...
> .libPaths()[1]
[1] "//Server02/stats/R/library/2.10"
... is controlled by R (in this example the former fails with "...not 
found", the latter opens the folder in Windows Explorer).

Best regards,

Keith Jewell
-
"Gabor Grothendieck"  wrote in message 
news:971536df1001260714i430b9ddcm37ffa12642321...@mail.gmail.com...

Try mapping a drive letter to \\Server02 and then use that.  For more,
google for: map network drive

On Tue, Jan 26, 2010 at 10:01 AM, Keith Jewell  
wrote:
> Hi,
>
>> update.packages(ask='graphics')
>
> gives me multiple warning (one per updated package?) similar to ...
>
> Warning: unable to move temporary installation
> '\\Server02\stats\R\library\2.10\file3de56e0d\locfit' to
> '\\Server02\stats\R\library\2.10\locfit'
>
> The final, updated, folders do not end up where they should be. I can move
> them 'by hand', but it is an inconvenience.
>
> Checking the archives I find
> http://finzi.psych.upenn.edu/Rhelp08/2008-April/160963.html where Prof
> Ripley opines "The issue appears to be that your OS is garbling file 
> names",
> and it surely does seem to be a filename problem - in combination R and
> Windows don't seem to be handling the '//server/share' path construction. 
> I
> have:
>> .libPaths()
> [1] "//Server02/stats/R/library/2.10"
> "//Server02/stats/R/R-Current/library"
>
> I can work around by mapping a drive letter to the '//server/share', but
> prefer not to (for local reasons).
>
> In CHANGES.R-2.10.1pat I find
> CHANGES IN R VERSION 2.7.2 patched
> o dir.create(recursive = TRUE) was not working on //server/share paths.
>
> In CHANGES.R-2.11.0dev I find
> CHANGES IN R VERSION 2.11.0
> NEW FEATURES
> o file.rename() can work across volumes (by copy-and-delete)
>
> In another R function I've hit a similar problem and solved it using the
> construction:
> file.path(dirname(x), basename(x))
> to convert '//server/share' to (printed as) 'server/share' which 
> worked
> in that function. In this case:
>> file.path(dirname(.libPaths()), basename(.libPaths()))
> [1] "Server02/stats/R/library/2.10"
> "Server02/stats/R/R-Current/library"
> which looks OK but
>> .libPaths(file.path(dirname(.libPaths()), basename(.libPaths(
>> .libPaths()
> [1] "//Server02/stats/R/library/2.10"
> "//Server02/stats/R/R-Current/library"
> doesn't actually change the internal paths.
>
> I don't see the problem in V2.9.2 .
> I do see the problem in
> V2.10.1,
> V2.10.1 Patched (2010-01-24 r51030)
> V2.11.0 Under development (unstable) (2010-01-24 r51030)
>
> Grateful for any suggestions,
>
> Keith Jewell
> -
> Version:
> platform = i386-pc-mingw32
> arch = i386
> os = mingw32
> system = i386, mingw32
> status = Patched
> major = 2
> minor = 10.1
> year = 2010
> month = 01
> day = 24
> svn rev = 51030
> language = R
> version.string = R version 2.10.1 Patched (2010-01-24 r51030)
>
> Windows Serve

[R] Linear regression and stand deviation at the Linux command line

[Keith Christian]

R List,

Please excuse this ultra-newbie post.
I looked at this page but it's a bit beyond me.
https://www2.kenyon.edu/Depts/Math/hartlaub/Math305%20Fall2011/R.htm

I'm interested in R construct(s) to be entered at the command
line that would output slope, y-intercept, and r-squared values read
from a csv or other filename entered at the command line, and the same
for standard deviation calculations, namely the standard deviation,
variance, and z-scores for every data point in the file.

E.g.
$ ((R function for linear regression here))slope, y-intercept, and
r-squared, other related stats that R seems most capable of
generating.
linear_regression_data.csv file contents (Are line numbers, commas,
etc. needed or no?)
1 20279
2 899
3 24747
4 12564
5 29543

$ ((R function for standard deviation here))standard deviation,
variance, z-scores, other related stats that R seems most capable of
generating.
standard_deviation_data.csv file contents (Are line numbers, commas,
etc. needed or no?)
1 16837
2 9498
3 31389
4 2365
5 17384

Many thanks,

--Keith

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

Re: [R] Linear regression and stand deviation at the Linux command line

[Keith Christian]

Hi Ivan,

Thanks for the suggestions.  Will try them.
 Keith

On Fri, Aug 23, 2024 at 1:57 AM Ivan Krylov  wrote:
>
> В Thu, 22 Aug 2024 13:07:37 -0600
> Keith Christian  пишет:
>
> > I'm interested in R construct(s) to be entered at the command
> > line that would output slope, y-intercept, and r-squared values read
> > from a csv or other filename entered at the command line, and the same
> > for standard deviation calculations, namely the standard deviation,
> > variance, and z-scores for every data point in the file.
>
> If you'd like to script R at the command line, consider the
> commandArgs() function (try entering ?commandArgs at the R prompt).
> This way you can pass a file path to an R process without unsafely
> interpolating it into the R expression itself. These arguments can be
> given to R --args or to Rscript (without the --args).
>
> Also consider the 'littler' scripting-oriented R front-end
> <https://CRAN.R-project.org/package=littler>, which puts the command
> line arguments into the 'argv' variable and has a very convenient -d
> option which loads CSV data from the standard input into a variable
> named 'X'.
>
> > Are line numbers, commas, etc. needed or no?
>
> Depends on how you read it. By default, the function read.table() will
> expect your data to be separated by a mixture of tabs and spaces and
> will recognise a header if the first line contains one less column than
> the rest of the file. Enter ?read.table at the R prompt to see the
> available options (which include read.csv).
>
> Good introductions to R include, well, "An Introduction to R" [1] (also
> available by typing RShowDoc('R-intro') into the R prompt) and "Visual
> Statistics" by Dr. A. Shipunov [2].
>
> Start with functions read.table(), lm(), scale(), sd(), summary(). Use
> str() to look at the structure of a variable: summary(lm(...)) will
> return a named list from which you can extract the values you are
> interested in (see ?summary.lm). When in doubt, call
> help(name_of_the_function).
>
> --
> Best regards,
> Ivan
>
> [1]
> https://cran.r-project.org/doc/manuals/R-intro.html
>
> [2]
> http://web.archive.org/web/20230106210646/http://ashipunov.info/shipunov/school/biol_240/en/visual_statistics.pdf

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

Re: [R] just a small variable-naming question

[Keith Jewell]

In case the OP wanted to append columns, rather than rename existing 
columns:


cbind(myData, var1=NA, var2=NA, var3=NA)
  col1 col2 col3 var1 var2 var3
1123   NA   NA   NA
2234   NA   NA   NA
3345   NA   NA   NA


On 24/08/2013 17:22, arun wrote:

Hi,
You could use ?paste()
colnames(myData)<-paste(colnames(myData),paste0("var",1:3),sep="_")
  myData
#  col1_var1 col2_var2 col3_var3
#1 1 2 3
#2 2 3 4
#3 3 4 5

A.K.



- Original Message -
From: Richard Sherman
To: r-help@r-project.org
Cc:
Sent: Saturday, August 24, 2013 12:09 PM
Subject: [R] just a small variable-naming question

Hi all,

If I have this little data set:


myData<- data.frame(col1 = 1:3, col2 =2:4, col3 = 3:5)
myData

   col1 col2 col3
1123
2234
3345

and I want to add (not replace) the names var1, var2, var3 to the existing 
col1, col2, col3, what is the right way to do that?

Thanks.

-Richard

---
Prof. Richard Sherman
Division of International Studies
Korea University

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] XLSX package + Excel creation question

[Keith Jewell]

I'll skip over the courtesy implications of double posting/pointing to 
stackoverflow.


The stackoverflow thread makes it look as if you need to learn more 
Excel. Do you really not know what an Excel template is?


It sounds as if you want what Excel calls "conditional formatting" which 
you can specify as custom number formats, see 
http://www.ozgrid.com/Excel/CustomFormats.htm.


Excel's help on custom number formats says:

To specify number formats that will be applied only if a number meets a 
condition that you specify, enclose the condition in square brackets. 
The condition consists of a comparison operator (comparison operator: A 
sign that is used in comparison criteria to compare two values. 
Operators include: = Equal to, > Greater than, < Less than, >= Greater 
than or equal to, <= Less than or equal to, and <> Not equal to.) and a 
value. For example, the following format displays numbers that are less 
than or equal to 100 in a red font and numbers that are greater than 100 
in a blue font.

[Red][<=100];[Blue][>100]
--

R package xlsx allows such formats (?DataFormat) as does R package 
XLConnect (?setDataFormat).


HTH

Keith J

On 04/09/2013 09:57, Zsurzsa Laszlo wrote:

http://stackoverflow.com/questions/18511249/excel-cell-coloring-using-xlsx

This is the initial post on stackoverflow. Please look at this maybe I'm
clearer here.

Thank you in advance,

-
- László-András Zsurzsa,-
- Msc. Infromatics, Technical University Munich, Germany -
- Scientific Employee, TUM -
-


On Fri, Aug 30, 2013 at 3:48 PM, jim holtman  wrote:


You can also look at the XLConnect package.
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.


On Thu, Aug 29, 2013 at 9:40 AM, Zsurzsa Laszlo
wrote:

I understand you response but it does not solve the problem. I'am aware
that one can simply color every cell in an excel file by using his own
algorithm.

The question was if I can write my data to a *single* cells and use
different formatting for every piece of data.



-

- László-András Zsurzsa,-
- Msc. Infromatics, Technical University Munich, Germany -
- Scientific Employee, TUM -


-



On Thu, Aug 29, 2013 at 3:36 PM, Rainer Hurling  wrote:


Am 29.08.2013 15:03 (UTC+1) schrieb Zsurzsa Laszlo:

First of all thank you for the quick resposen.

I know I can color and set up every cell. I will take a look again *
CellStyle* but is it possbile for example to write an array to a

single

cell that has different colors for some data. Basically the color

depends

on the data.


As far as I know there is no ready to use functionality to mask groups
of selected cells. You have to write your own function, which selects
the right cells and changes their style with setCellStyle(cell,

cellStyle).


Some hints are given in the examples section of ?CellStyle.







-

- László-András Zsurzsa,

-

- Msc. Infromatics, Technical University Munich, Germany -
- Scientific Employee, TUM

   -





-



On Thu, Aug 29, 2013 at 2:55 PM, Rainer Hurling

wrote:



Am 29.08.2013 12:08 (UTC+1) schrieb Zsurzsa Laszlo:

Dear R users,

I have a question about the xlsx package. It's possible to create

excel

files and color cells and etc.


yes, with package xlsx you can colourize you data sheets, even the
fonts. See for example ?CellStyle .

A good demonstration of the capabilities is on





http://tradeblotter.wordpress.com/2013/05/02/writing-from-r-to-excel-with-xlsx/




My question would be that is it possible to color only some part of

the

data hold in a cell. Let's assume I've got the following data :
167,153,120,100 and I want to color to red everything that is bigger

then

120. How can I achive this using R.

Example file setup with a few lines in attachment. (SEL_MASS column

can

be

used for example)


Attachment missing ...

HTH,
Rainer



Thank you in advance,






-

- László-András Zsurzsa,

  -

- Msc. Infromatics, Technical University Munich, Germany -
- Scientific Employee, TUM

-







--

Re: [R] Problem with converting F to FALSE

[Keith Jewell]

Depending what you're doing with the data, you might want
colClasses=c("factor","numeric")

On 05/09/2013 13:58, Joshua Wiley wrote:

Hi,

You can either manually specify colClasses or the asis argument.  See
?read.csv for more details.

If you just had those two columns, something like:

  read.table(header = TRUE, text = "
  sex group
  F 1
  T 2
  ", colClasses = c("character", "integer"))

Cheers,

Josh


read.csv("file.csv", colClasses = c("character", "integer"))




On Thu, Sep 5, 2013 at 5:44 AM, Venkata Kirankumar
  wrote:

Hi,
I have a peculier problem in R-Project that is when my CSV file have one
column with all values as 'F' the R-Project converting this 'F' to FALSE.
Can some one please suggest how to stop this convertion. Because I want to
use 'F' in my calculations and show it in screen. for example my data is
like

sex  group
F   1
F   2
F   3

but when I use read.csv and load the csv file data is converting it to

sex  group
FALSE   1
FALSE   2
FALSE   3
but i want it as source data like

sex group
F  1
F  2
F  3


Thanks in advance,
D V Kiran Kumar


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

Re: [R] function coverage

[Keith Jewell]

On 14/01/2013 22:25, Hadley Wickham wrote:

I think codetools could do this reasonably well with the walkCode function,
but I've never done it so I don't have sample code, and walkCode is mostly
an internal function.


There are a couple of approaches here:
http://stackoverflow.com/questions/14276728/

Hadley


I've used foodweb in mvbutils for that kind of thing.

HTH

KJ

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

[R] How do I get Garamond font in R?

[Keith Weintraub]

Folks,
  I run R on a early 2009 MacBook Pro running Mountain Lion.

I have a bunch of fonts in my user Library one of which is Garamond.

I have tried the ttf_import function to no avail. I played with this for a 
couple of hours at least and I have gotten nowhere.

Here is a bit of one of my sessions (note that I used file.choose to get down 
to the Font level)

> require(extrafont)
Loading required package: extrafont
Registering fonts with R
Warning messages:
1: In loadfonts("pdf", quiet = TRUE) :
  More than one version of regular/bold/italic found for Apple Braille. 
Skipping setup for this font.
2: In loadfonts("pdf", quiet = TRUE) :
  No regular (non-bold, non-italic) version of Brush Script MT. Skipping setup 
for this font.
3: In loadfonts("postscript", quiet = TRUE) :
  More than one version of regular/bold/italic found for Apple Braille. 
Skipping setup for this font.
4: In loadfonts("postscript", quiet = TRUE) :
  No regular (non-bold, non-italic) version of Brush Script MT. Skipping setup 
for this font.

> fonts()
 [1] ".Keyboard" "Andale Mono"   "Apple Braille" 
"AppleMyungjo" 
 [5] "Arial Black"   "Arial" "Arial Narrow"  
"Arial Rounded MT Bold"
 [9] "Arial Unicode MS"  "Batang""Brush Script MT"   
"Calibri"  
[13] "Cambria"   "Candara"   "Comic Sans MS" 
"Consolas" 
[17] "Constantia""Corbel""Courier New"   
"Georgia"  
[21] "Gulim" "Humor Sans""Impact"
"Jazz" 
[25] "JazzCord"  "JazzPerc"  "JazzText"  
"JazzTextExtended" 
[29] "Khmer Sangam MN"   "Krini" "Krinitsky" 
"Lao Sangam MN"
[33] "Microsoft Sans Serif"  "MS Gothic" "MS Mincho" 
"MS PGothic"   
[37] "MS PMincho""MusicalSymbols""Myanmar Sangam MN" 
"PMingLiU" 
[41] "SimSun""Sinfonia"  "Tahoma"
"Times New Roman"  
[45] "Trebuchet MS"  "Verdana"   "Webdings"  
"Wingdings"
[49] "Wingdings 2"   "Wingdings 3"  

> ttf_import(file.choose())
Scanning ttf files in /Users/kw/Library/Fonts/Garamond ...
Extracting .afm files from .ttf files...
Error in data.frame(fontfile = ttfiles, FontName = "", stringsAsFactors = 
FALSE) : 
  arguments imply differing number of rows: 0, 1

Here is what Font Book has to say about the font:

Garamond
Garamond

PostScript name Garamond
Full name   Garamond
Family  Garamond
Style   Regular
KindTrueType
LanguageAfrikaans, Albanian, Basque, Cornish, Danish, Dutch, 
English, Faroese, French, Galician, German, Icelandic, Indonesian, Irish, 
Italian, Malay, Manx, Norwegian Bokmål, Norwegian Nynorsk, Oromo, Portuguese, 
Somali, Spanish, Swahili, Swedish, Swiss German, Zulu
Script  Latin
Version Version 2.35
Location/Users/kw/Library/Fonts/Garamond
Unique name Monotype - Garamond Regular
ManufacturerMonotype Typography, Inc.
DesignerClaude Garamond
Copyright   Digitized data copyright Monotype Typography, Ltd 
1991-1995. All rights reserved. Monotype Garamond® is a trademark of Monotype 
Typography, Ltd which may be registered in certain jurisdictions.
Trademark   Monotype Garamond® is a trademark of Monotype 
Typography, Ltd which may be registered in certain jurisdictions.
Description Monotype Drawing Office 1922. This typeface is based on 
roman types cut by Jean Jannon in 1615. Jannon followed the designs of Claude 
Garamond which had been cut in the previous century. Garamond's types were, in 
turn, based on those used by Aldus Manutius in 1495 and cut by Francesco 
Griffo. The italic is based on types cut in France circa 1557 by Robert 
Granjon. Garamond is a beautiful typeface with an air of informality which 
looks good in a wide range of applications. It works particularly well in books 
and lengthy text settings.
License NOTIFICATION OF LICENSE AGREEMENT

This typeface is the property of Monotype Typography and its use by you is 
covered under the terms of a license agreement. You have obtained this typeface 
software either directly from Monotype or together with software distributed by 
one of Monotype’s licensees.

This software is a valuable asset of Monotype. Unless you have entered into a 
specific license agreement granting you additional rights, your use of this 
software is limited to your workstation for your own publishing use. You may 
not copy or distribute this software.

If you have any question concerning your rights you should review the license 
agreement you received with the software or contact Monotype for a co

Re: [R] Query on package to use for investment quotes

[Keith Weintraub]

Bruce,
  Have you looked at the Edgar database here: http://www.sec.gov/edgar.shtml

You should be able to get daily mutual fund quotes.

KW

--


[[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] Packages with functionality related to Oil/Gas exploration

[Keith Weintraub]

Folks,

As the subject describes: I would like to know if there are packages that have 
functionality tailored for standard Oil/Gas exploration and monetization.

Thanks,
KW

--


[[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] approxfun values

[Keith Jewell]

Alternatively, with approx() use xout to specify which interpolated 
values you want returned:


approx(dat, xout=dat$V1[is.na(dat$V2)])

KJ

On 14/02/2013 11:43, Rui Barradas wrote:

Hello,

In what follows I've changed your df name to 'dat', to save some
keystrokes.
approxfun returns a function, so if you want just the interpolated
values, apply that function to the x values where y is NA.



dat <- read.table(text = "
V1 V2
1 10 2
2 20 NA
3 30 5
4 40 7
5 50 NA
6 60 NA
7 70 2
8 80 6
9 90 9
10 100 NA
", header = TRUE)

f <- approxfun(dat)
x <- dat$V1[is.na(dat$V2)]
y <- f(x)
y


Hope this helps,

Rui Barradas

Em 14-02-2013 08:43, e-letter escreveu:

Readers,

According to the help '?approxfun', the function can be used to obtain
the interpolated values. The following test was tried:


testinterpolation<-read.csv('test.csv',header=FALSE)
testinterpolation

V1 V2
1 10 2
2 20 NA
3 30 5
4 40 7
5 50 NA
6 60 NA
7 70 2
8 80 6
9 90 9
10 100 NA

testinterpolationvalues<-approxfun(testinterpolation,y=NULL)
testinterpolationvalues

function (v)
.C(C_R_approxfun, as.double(x), as.double(y), as.integer(n),
xout = as.double(v), as.integer(length(v)), as.integer(method),
as.double(yleft), as.double(yright), as.double(f), NAOK = TRUE,
PACKAGE = "stats")$xout



testinterpolationvalues<-approx(testinterpolation,y=NULL)
testinterpolationvalues

$x
[1] 10.0 11.63265 13.26531 14.89796 16.53061 18.16327 19.79592
21.42857
[9] 23.06122 24.69388 26.32653 27.95918 29.59184 31.22449 32.85714
34.48980
[17] 36.12245 37.75510 39.38776 41.02041 42.65306 44.28571 45.91837
47.55102
[25] 49.18367 50.81633 52.44898 54.08163 55.71429 57.34694 58.97959
60.61224
[33] 62.24490 63.87755 65.51020 67.14286 68.77551 70.40816 72.04082
73.67347
[41] 75.30612 76.93878 78.57143 80.20408 81.83673 83.46939 85.10204
86.73469
[49] 88.36735 90.0

$y
[1] 2.00 2.244898 2.489796 2.734694 2.979592 3.224490 3.469388
3.714286
[9] 3.959184 4.204082 4.448980 4.693878 4.938776 5.244898 5.571429
5.897959
[17] 6.224490 6.551020 6.877551 6.829932 6.557823 6.285714 6.013605
5.741497
[25] 5.469388 5.197279 4.925170 4.653061 4.380952 4.108844 3.836735
3.564626
[33] 3.292517 3.020408 2.748299 2.476190 2.204082 2.163265 2.816327
3.469388
[41] 4.122449 4.775510 5.428571 6.061224 6.551020 7.040816 7.530612
8.020408
[49] 8.510204 9.00

How to obtain a vector consisting _only_ of the interpolated values?

It was expected that 'approx' would return both original and
interpolated values (which the above shows) and that 'approxfun' would
not show the original values.

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

[R] What happened to "financial" package?

[Keith Weintraub]

Folks,
  I have been using the cf function from the "financial" package for a few 
years now.

Upon updating my version of R to ... I found that the package no longer exists 
in the main collection of R packages.

Has this package been renamed or merged with another package?

If the functionality is no longer available what is the best way to integrate 
the archived version of the package into my current setup?

Thanks for your time,
KW


--


[[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] What happened to "financial" package?

[Keith Weintraub]

Dr. R,
  Thanks for your help. I had already gotten the archive package and sourced 
the code.

I was just wondering if there was a *best* way to make it a kind of "local" 
package.

I hope this message is in plain text. I thought I changed my settings for this 
awhile ago. If there is still a problem please let me know.

Also I read the digests daily. What is the best way to reply to a message in 
the digest?

Thanks for the help,
KW



> On 26/02/2013 14:45, Keith Weintraub wrote:
>> Folks,
>>   I have been using the cf function from the "financial" package for a few 
>> years now.
> 
> I believe it was never converted for 2.14.0 (it was last touched in 
> 2006, and did not have a NAMESPACE file until CRAN auto-generated one), 
> and the maintainer address is invalid.
> 
>> Upon updating my version of R to ... I found that the package no longer 
>> exists in the main collection of R packages.
>> 
>> Has this package been renamed or merged with another package?
> 
> If so, the CRAN package page would say so.
> 
>> If the functionality is no longer available what is the best way to 
>> integrate the archived version of the package into my current setup?
> 
> Get the package from the archive, extract the file cf.R and source it.
> 
>> Thanks for your time,
>> KW
>> 
>> 
>> --
>> 
>> 
>>  [[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.
>> 
> That does mean you: no HTML.

--

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

[R] string split at xth position

[Keith Weintraub]

Here is another way

require(stringr)

aaa<-paste0("a", 1:20)
bbb<-paste0("b", 101:120)
ab<-paste0(aaa,bbb)
ab


ptrn<-"([ab][[:digit:]]*)"
unlist(str_extract_all(ab, ptrn))



> Hi,
> 
> I have a vector of strings like:
> c("a1b1","a2b2","a1b2") which I want to spilt into two parts like:
> c("a1","a2","a2") and c("b1","b2,"b2"). So there is
> always a first part with a+number and a second part with b+number.
> Unfortunately there is no separator I could use to directly split
> the vectors.. Any idea how to handle such cases?
> 
> /Johannes


--

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

Re: [R] nls: example code throws error

[Keith Jewell]

On 26/04/2013 00:16, Steven LeBlanc wrote:
> Greets,
>
> I'm trying to learn to use nls and was running the example code for 
an exponential model:

>
>  
>
> Perhaps also, a pointer to a comprehensive and correct document that 
details model formulae syntax if someone has one?

>
> Thanks&  Best Regards,
> Steven
>
Others have pointed out that the error is probably from an unclean 
environment.


For model formula syntax, see ?nls
Under Arguments formula, follow the link to ?formula

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

[Keith Weintraub]

Folks,
I am working on a credit card defaults and transition probabilities. For 
example a single credit card account could be in a number of states: 
up-to-date, 30, 60, 90 days in arrears or in default.

* Are there packages in R that do estimation of the transition probabilities 
given historical cohort defaults?
* Any pointers to papers specific to this type of estimation?
* Simulation of future "paths" of defaults.

I know this is not strictly an R question so please feel free to slam me.

Afterwards I would appreciate any pointers you might have.

Thanks much for your time,
KW

--


[[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] Equation of a curve

[Keith Jewell]

On 03/04/2014 16:26, Frances Cheesman wrote:

Hi all,

I have a number of bacterial growth curves I would like to find the
equations for these and then integrate them to find the area under the
curves for me to do stats on later.

Is there any way I can do this in R?

Thanks,

Frances

[[alternative HTML version deleted]]

Responding to the curve fitting question and passing over the 
integration issue...


It is quite common to use nls to fit equations to log(count) v time 
data. You'll have to choose an appropriate model, ideally as a self 
starting nls model. Of those included in the stats package you might 
consider SSfpl, SSgompertz, SSlogis and SSweibull.


But choice of a model is really a microbiological issue and all those 
models might be considered a little passe. Fitting this kind of 
sigmoidal model can be difficult unless the data is good.


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

[Keith Weintraub]

Folks,
  I have a data frame with columns 200401, 200402, ..., 201207, 201208.

These represent years/months. What would be the best way to sum these columns 
by year? What about by quarter?

Thanks for your time,
KW

--


[[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] Aggregate data frame across columns

[Keith Weintraub]

Arun, Jeff, Bert,
   Thanks for your help.

I have put a subset of my data below in mySubset.

I would like to be able to sum the rows by year. In this case the results would 
be the result data.frame below.

How can I automate something like this and how would I do it quarterly if 
necessary.

The data has 9 columns with headers:
   "X200401" "X200402" "X200501" "X200502" "X200503" "X200601" "X200602" 
"X200701" "X200702"

Here is the code I used to create the result: 
  result<-data.frame(rowSums(mySubset[,1:2], na.rm = TRUE), 
rowSums(mySubset[,3:5], na.rm = TRUE), rowSums(mySubset[,6:7], na.rm = TRUE), 
rowSums(mySubset[,8:9], na.rm = TRUE))

Note that my full dataset goes from 200401 through 201208.

I assume there is a way to use one of the apply functions and an inline 
function to call rowSums with na.rm = TRUE.

Thanks again for your time,
KW

__
Two datasets below: result and mySubset.

result<-structure(list(X2004 = c(0.0159924882870401, 0.0914601927232432, 
0.138321748009262, 0.156063783591084, 0.168403383789346, 0.171965759793573, 
0.177147721902362, 0.187522847481161, 0.16541728156, 0.127352907374406, 
0.156908213362621, 0.175500673945803, 0.17516598558791, 0.176671361535308, 
0.174478461658455, 0.157756648535001, 0.180489661678831, 0.189127686535455, 
0.176267288362896, 0.167844722339248, 0.180507725071878, 0.169459551401114, 
0.165939970730443, 0.165709877723436, 0.17229145356651, 0.182795171134028, 
0.166283818929029, 0.15294456192766, 0.166780783496174, 0.181927809974243, 
0.177579619132214, 0.171811823922994, 0.158385247734671, 0.149479196737791, 
0.162477792074099, 0.150845832508427, 0.155104452310268, 0.162456727168325, 
0.155482804148341, 0.138967760361165, 0.146530081719247, 0.157417284793283, 
0.15859793000523, 0.14774834433617, 0.147320895948278, 0.14926677799197, 
0.14171723173142, 0.136909644046266, 0.0683144142254727, 0), 
X2005 = c(0.0314065680219376, 0.163832345277566, 0.265371909518265, 
0.326014428812549, 0.368234576027844, 0.408135729854325, 
0.406609083944697, 0.416411672384248, 0.383445330771284, 
0.303908689700078, 0.279111432968148, 0.299088365053801, 
0.297231620329575, 0.301554936855911, 0.286267190479442, 
0.285101593132908, 0.303617540120288, 0.30912628079129, 0.306740852364357, 
0.288318399163201, 0.290849771380406, 0.284963182810939, 
0.276913327940756, 0.281228989779383, 0.281587206458797, 
0.282575344784775, 0.272750026504263, 0.275012113477194, 
0.275308072336096, 0.2830402877904, 0.278126404864579, 0.270620664127706, 
0.266703443412622, 0.269170883813461, 0.267706882000802, 
0.16400030580057, 0.081430630811921, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0), X2006 = c(0.0263421114412355, 0.189060034450041, 
0.307867112279838, 0.36147798027164, 0.421244020838752, 0.443683746556757, 
0.431455767826625, 0.42139989322102, 0.405361335329633, 0.410239754621491, 
0.295699686523215, 0.286907309706468, 0.279272095208596, 
0.30584849030251, 0.292930451848863, 0.2858086693289, 0.277977449378127, 
0.312603778277765, 0.285799481868903, 0.272257225207343, 
0.264193812603326, 0.260225982968127, 0.26286016640048, 0.246455494757128, 
0.132509032599313, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), X2007 = c(0.0192634102835919, 
0.0821176964052597, 0.125895446649419, 0.166105172058653, 
0.193371112440694, 0.192128664695495, 0.205402741722323, 
0.214804686886319, 0.206944760826543, 0.201633705240255, 
0.224370691628553, 0.213780729123737, 0.116800798312322, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)), .Names = c("X2004", 
"X2005", "X2006", "X2007"), row.names = c("1", "2", "3", "4", 
"5", "6", "7", "8", "9", "10", "65", "66", "67", "68", "69", 
"70", "71", "72", "73", "74", "75", "76", "77", "78", "79", "80", 
"81", "82", "83", "84", "85", "86", "87", "88", "89", "90", "91", 
"92", "93", "94", "95", "96", "97", "98", "99", "100", "101", 
"102", "103", "104"), class = "data.frame")

mySubset<-structure(list(X200401 = c(0.00989521347529282, 0.0373090976259132, 
0.0589205312533588, 0.073423408836654, 0.0684862018994913, 0.0709590032723634, 
0.0659646021314206, 0.066731059852463, 0.0699763372273688, 0.0541142327051023, 
0.0747822118456574, 0.0796152619493638, 0.0898395772710088, 0.0799174340369894, 
0.0858618669044432, 0.0771658556270448, 0.0737537341871075, 0.0920394845122411, 
0.0919450414276911, 0.0766003513927418, 0.0853798302884387, 0.0906216842117783, 
0.0770655693805709, 0.0804880021458966, 0.0793736980924377, 0.0867350210905336, 
0.0844120272331623, 0.0719191481914751, 0.071331576585233, 0.0814113751491231, 
0.089572460115355, 0.078571913014948, 0.080089664942495, 0.0660045583281544, 
0.0757219707900531, 0.0726744688001806, 0.0696490161975653, 0.0771202956076235, 
0.0762013529234614, 0.0699300362070559, 0.0671216348093859, 0.0733447602438439, 
0

[R] Creating yyyymm regexp strings on the fly for aggregation.

[Keith Weintraub]

Folks,
  This question is somewhat related to a previous posting of mine.

I just can't seem to create a generic solution.

Here is a function that I found searching around the internet:
splitIt <- function(x, n) {split(x, sort(rank(x) %% n))}

I use it like so:
> splitIt(1:12, 2)
$`0`
[1] 1 2 3 4 5 6

$`1`
[1]  7  8  9 10 11 12

Or
> splitIt(1:12, 4)
$`0`
[1] 1 2 3

$`1`
[1] 4 5 6

$`2`
[1] 7 8 9

$`3`
[1] 10 11 12

I am splitting 12 months into 6-month or quarterly chunks. I can also use the 
function to create monthly or segments or one big annual segment.

Here is a function that I developed:

groupingStrings<-function(yrs, numSplits) {
 unlist(lapply(yrs, function(x){ paste(x,formatC(numSplits, width = 2, flag = 
0), collapse = "|", sep = "")}))
}

Here is an example of running the function:
groupingStrings(2004:2006, 1:3)
[1] "200401|200402|200403" "200501|200502|200503" "200601|200602|200603"

This would yield first quarter matches for the years 2004 through 2006.

My plan was to use both splitIt and groupingStrings to be able to create 
regexps  all quarters.

In addition I want it to be flexible enough for me to be able to create 
matching regexps for monthly, quarterly, semi-annually and annual regexps.

One more example. Suppose I wanted to look at data semi-annually for 2010 
through 2011. The regexps would be:
 "201001|201002|201003|201004|201005|201006"
 "201007|201008|201009|201010|201011|201012"
 "201101|201102|201103|201104|201105|201106"
 "201107|201108|201109|201110|20|201112"

I hope I have explained my problem clearly.

Thanks much for your time,
KW


--


[[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] Creating yyyymm regexp strings on the fly for aggregation.

[Keith Weintraub]

Thanks,
   I feel like I was so close!!!

KW

--

On Nov 9, 2012, at 4:48 AM, Rui Barradas  wrote:

> Hello,
> 
> Something like this?
> 
> 
> fun <- function(years, periods){
>res <- lapply(splitIt(1:12, periods), function(x) groupingStrings(years, 
> x))
>names(res) <- as.integer(names(res)) + 1L
>res
> }
> 
> fun(2010:2011, 1)
> fun(2010:2011, 2)
> fun(2010:2011, 4)
> fun(2010:2011, 12)
> 
> 
> Hope this helps,
> 
> Rui Barradas
> Em 09-11-2012 03:39, Keith Weintraub escreveu:
>> Folks,
>>   This question is somewhat related to a previous posting of mine.
>> 
>> I just can't seem to create a generic solution.
>> 
>> Here is a function that I found searching around the internet:
>> splitIt <- function(x, n) {split(x, sort(rank(x) %% n))}
>> 
>> I use it like so:
>>> splitIt(1:12, 2)
>> $`0`
>> [1] 1 2 3 4 5 6
>> 
>> $`1`
>> [1]  7  8  9 10 11 12
>> 
>> Or
>>> splitIt(1:12, 4)
>> $`0`
>> [1] 1 2 3
>> 
>> $`1`
>> [1] 4 5 6
>> 
>> $`2`
>> [1] 7 8 9
>> 
>> $`3`
>> [1] 10 11 12
>> 
>> I am splitting 12 months into 6-month or quarterly chunks. I can also use 
>> the function to create monthly or segments or one big annual segment.
>> 
>> Here is a function that I developed:
>> 
>> groupingStrings<-function(yrs, numSplits) {
>>  unlist(lapply(yrs, function(x){ paste(x,formatC(numSplits, width = 2, flag 
>> = 0), collapse = "|", sep = "")}))
>> }
>> 
>> Here is an example of running the function:
>> groupingStrings(2004:2006, 1:3)
>> [1] "200401|200402|200403" "200501|200502|200503" "200601|200602|200603"
>> 
>> This would yield first quarter matches for the years 2004 through 2006.
>> 
>> My plan was to use both splitIt and groupingStrings to be able to create 
>> regexps  all quarters.
>> 
>> In addition I want it to be flexible enough for me to be able to create 
>> matching regexps for monthly, quarterly, semi-annually and annual regexps.
>> 
>> One more example. Suppose I wanted to look at data semi-annually for 2010 
>> through 2011. The regexps would be:
>>  "201001|201002|201003|201004|201005|201006"
>>  "201007|201008|201009|201010|201011|201012"
>>  "201101|201102|201103|201104|201105|201106"
>>  "201107|201108|201109|201110|20|201112"
>> 
>> I hope I have explained my problem clearly.
>> 
>> Thanks much for your time,
>> KW
>> 
>> 
>> --
>> 
>> 
>>  [[alternative HTML version deleted]]
>> 
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
> 


[[alternative HTML version deleted]]

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

[R] Question about VAR (Vector Autoregression) in differences.

[Keith Weintraub]

Folks,
  I have been using the VAR {vars} program to find a fit for the following 
bi-variate time series (subset):

bivariateTS<-structure(c(0.950415958293559, 0.96077848972081, 
0.964348957109053, 
0.967852884998915, 0.967773510751625, 0.970342843688257, 0.97613937178359, 
0.980118627997436, 0.987059493773907, 0.99536830931504, 1.00622672085718, 
1.01198013845981, 1.01866618122606, 1.02039666233818, 1.02400374633675, 
1.0249388841, 1.02865185614599, 1.03337978432753, 1.03244791539977, 
1.03237166203917, 1.02853704067597, 1.0271704031946, 1.02588255626357, 
1.02300209859012, 1.02397946952311, 1.02052043198189, 1.01593679111077, 
1.01535467405129, 1.01368552421158, 1.00833115261092, 1.00495328099247, 
1.00334161454411, 1.00029163432818, 0.999268774926758, 0.998174371988049, 
0.999643729403106, 1.00255061235855, 1.0036328278581, 1.00343133578339, 
1.00020064935273, 0.478881778679413, 0.404641679179503, 0.622804024030457, 
0.474677494404522, 0.433851413612414, 0.544920292447026, 0.561093836992123, 
0.563908924197768, 0.395583533713216, 0.273458857040352, 0.284537535240248, 
0.277701982320208, 0.451280071722745, 0.401744885050984, 0.533760449322162, 
0.488858723259857, 0.426799781286085, 0.503196832449693, 0.526637755320189, 
0.599897302279671, 0.519986312138689, 0.460526054741527, 0.550734747489121, 
0.487644621564283, 0.524566370547446, 0.700614666580712, 0.620626909353966, 
0.661635856684872, 0.631149920530168, 0.668383815717266, 0.646705357249337, 
0.68413824799, 0.650175601510345, 0.679512593649786, 0.67242973935153, 
0.599346651479619, 0.630852735978974, 0.657541680143728, 0.694629910826878, 
0.654572164624051), .Dim = c(40L, 2L), .Dimnames = list(NULL, 
c("a", "b")))

I have used the "usual" techniques to show that both of the series are I(1), 
(integrated of order 1) AND that there is no co-integration. In addition I used 
the VARselect program to estimate the number of lags.

To estimate the relationship between these variables the standard recipe (as 
far as I know) is to take first differences of the individual series and apply 
the VAR program to the new bivariate difference series.

The VAR program seems to derive a reasonable set of coefficients.

My Question:
How do I simulate future paths for the original pair of (undifferenced) series?

Thanks for your time and help,
KW



The rest of this note is just how I have done this up until now...

The method I have been using is to take a simulated path of the differences 
(note that in my case there are 2 simulated paths) say
> x.d
[1]   9  90 900

and then add back the starting value (1 in this case)

> c(1, 1+cumsum(x.d))
[1]1   10  100 1000 

A kind of re-integration of the series.

In a real case I would use the final value(s) of bivariateTS (1.0002006, 
0.6545722 ) in place of the "1" in the cumsum above.

My problem:
The paths that I get make no sense.

Here is the function I use (an alteration of Professor Pfaff's predict.varest) 
to simulate var paths. My apologies if it is bad form to use alter someone 
else's code for other purposes.

# Modified VAR simulation function jigged up from the predict function of vars 
package
simVARpath<-function(object, n.ahead, mult = 1) {
K <- object$K
p <- object$p
obs <- object$obs
type <- object$type
data.all <- object$datamat
ynames <- colnames(object$y)
n.ahead <- as.integer(n.ahead)
Z <- object$datamat[, -c(1:K)]
B <- Bcoef(object)
if (type == "const") {
Zdet <- matrix(rep(1, n.ahead), nrow = n.ahead, ncol = 1)
colnames(Zdet) <- "const"
}
else if (type == "trend") {
trdstart <- nrow(Z) + 1 + p
Zdet <- matrix(seq(trdstart, length = n.ahead), nrow = n.ahead,
ncol = 1)
colnames(Zdet) <- "trend"
}
else if (type == "both") {
trdstart <- nrow(Z) + 1 + p
Zdet <- matrix(c(rep(1, n.ahead), seq(trdstart, length = n.ahead)), 
nrow = n.ahead, ncol = 2)
colnames(Zdet) <- c("const", "trend")
}
else if (type == "none") {
Zdet <- NULL
}
if (!is.null(eval(object$call$season))) {
season <- eval(object$call$season)
seas.names <- paste("sd", 1:(season - 1), sep = "")
cycle <- tail(data.all[, seas.names], season)
seasonal <- as.matrix(cycle, nrow = season, ncol = season - 
1)
if (nrow(seasonal) >= n.ahead) {
seasonal <- as.matrix(cycle[1:n.ahead, ], nrow = n.ahead, 
ncol = season - 1)
}
else {
while (nrow(seasonal) < n.ahead) {
seasonal <- rbind(seasonal, cycle)
}
seasonal <- seasonal[1:n.ahead, ]
}
rownames(seasonal) <- seq(nrow(data.all) + 1, length = n.ahead)
if (!is.null(Zdet)) {
Zdet <- as.matrix(cbind(Zdet, seasonal))
}
else {
Zdet <- as.matrix(seasonal)
}
}
if (!is.null(eval(objec

Re: [R] Effect of each term in the accuracy of Nonlinear multivariate regression fitting equation

[Keith Jewell]

In this context, "linear model" means linear in the _coefficients_ not 
(necessarily) linear in the predictors, so your model:

   JIM ~ z1*A + z2*B + z3*A*B^2 + z4*C*D^3 + z5*A^2*B^2 ...
is a linear model (in z1, z2, ...).

So you don't need to use nls, lm is probably favourite. You can use all 
the techniques around for evaluating linear models; anova.lm might give 
you a start.


KJ

On 27/11/2012 11:40, dsfakianakis wrote:

Dear all,

I have a set of data with 4 inputs (independent variables) and one output
(dependent variable). I want to perform a regression analysis in order to
fit these data to a regression model, however due to the non-linearity of
the model I do not have a clue which equation to use. I am thinking of
starting with a very general equation including ^3 terms and interactions
between the variables however this will lead to a very long equation. Is
there a way to assess the effect of each term to the accuracy of the
regression model in order to discard the terms with the least importance?
Something like a sensitivity analysis of the effect of each term to the
accuracy regression model. I know one possible solution to my problem is
simply 'trial and error' however before going down that road I want to check
if there is an easier way.

e.g. Let's say I have four input variables A B C and D, one output 'JIM' and
let z1, z2, ...  be the coefficients of the terms of the equation.  The
regression will be something like that:

Result = nls(JIM ~ z1*A + z2*B + z3*A*B^2 + z4*C*D^3 + z5*A^2*B^2 ... )

Is there a way to assess the contribution of each term (z1*A, z3*A*B^3 etc)
to the accuracy of the regression model?

Thanks a lot



--
View this message in context: 
http://r.789695.n4.nabble.com/Effect-of-each-term-in-the-accuracy-of-Nonlinear-multivariate-regression-fitting-equation-tp4650949.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] grep for multiple pattern?

[Keith Jewell]

On 13/02/2014 15:51, Marc Schwartz wrote:


On Feb 13, 2014, at 8:43 AM, Rainer M Krug  wrote:


Hi

I want to search for multiple pattern as grep is doing for a single
pattern, but this obviously not work:


grep("an", month.name)

[1] 1

grep("em", month.name)

[1]  9 11 12

grep("eb", month.name)

[1] 2

grep(c("an", "em", "eb"), month.name)

[1] 1
Warning message:
In grep(c("an", "em", "eb"), month.name) :
  argument 'pattern' has length > 1 and only the first element will be used




Is there an equivalent which returns the positions as grep is doing, but
not using the strict full-string matching of match()?

I could obviously do:


unlist( sapply(pat, grep, month.name ) )

an em1 em2 em3  eb
  1   9  11  12   2

but is there a more compact command I am missing?

Thanks,

Rainer



The vertical bar '|' acts as a logical 'or' operator in regex expressions:


grep("an|em|eb", month.name)

[1]  1  2  9 11 12


grep("an|em|eb", month.name, value = TRUE)

[1] "January"   "February"  "September" "November"  "December"


Regards,

Marc Schwartz


and if you want your patterns in a vector
> pat <-c("an", "em", "eb")
> grep(paste(pat, collapse="|"), month.name)
[1]  1  2  9 11 12

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

Re: [R] princomp/prcomp packages not available for 3.0.2

[Keith Jewell]

On 19/02/2014 15:47, Rich Shepard wrote:

   Running 3.0.2 on Slackware here. Tried to install.packages() for both
princomp and prcomp (Principal Components Analysis) but R responded by
telling me neither is available for this version of R.

   Has anyone an idea when either (but especially prcomp) might be
available?

TIA,

Rich

princomp and prcomp are not packages, they are functions in the 'stats' 
package.


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

Re: [R] A vector of normal distributed values with a sum-to-zero constraint

[Keith Jewell]

It seems so simple to me, that I must be missing something.

Subject to Jeff Newmiller's reminder of FAQ 7.31; if the sum is zero 
then the mean is zero and vice versa.


The OP's original attempt of:
-
l <- 100
aux <- rnorm(l,0,0.5)
s <- sum(aux)/l
aux2 <- aux-s
sum(aux2)
-
is equivalent to

  aux2 <- rnorm(l,0,0.5)
  aux2 <- aux2-mean(aux2)

If calculations were exact then aux2 would have mean, and thus sum, 
equal to zero - any difference from zero is attributable entirely to 
machine precision.



On 01/04/2014 15:25, Boris Steipe wrote:

But the result is not Normal. Consider:

set.seed(112358)
N <- 100
x <- rnorm(N-1)
sum(x)

[1] 1.759446   !!!

i.e. you have an outlier at 1.7 sigma, and for larger N...

set.seed(112358)
N <- 1
x <- rnorm(N-1)
sum(x)
[1] -91.19731

B.


On 2014-04-01, at 10:14 AM, jlu...@ria.buffalo.edu wrote:


The sum-to-zero constraint imposes a loss of one degree of freedom.  Of  N 
samples, only (N-1) can be random.   Thus the solution is

N <- 100
x <- rnorm(N-1)
x <- c(x, -sum(x))
sum(x)

[1] -7.199102e-17












Boris Steipe 
Sent by: r-help-boun...@r-project.org
04/01/2014 09:29 AM

To
Marc Marí Dell'Olmo ,
cc
"r-help@r-project.org" 
Subject
Re: [R] A vector of normal distributed values with a sum-to-zero
constraint





Make a copy with opposite sign. This is Normal, symmetric, but no longer random.

  set.seed(112358)
  x <- rnorm(5000, 0, 0.5)
  x <- c(x, -x)
  sum(x)
  hist(x)

B.

On 2014-04-01, at 8:56 AM, Marc Marí Dell'Olmo wrote:


Dear all,

Anyone knows how to generate a vector of Normal distributed values
(for example N(0,0.5)), but with a sum-to-zero constraint??

The sum would be exactly zero, without decimals.

I made some attempts:


l <- 100
aux <- rnorm(l,0,0.5)
s <- sum(aux)/l
aux2 <- aux-s
sum(aux2)

[1] -0.0006131392


aux[1]<- -sum(aux[2:l])
sum(aux)

[1] -0.03530422


but the sum is not exactly zero and not all parameters are N(0,0.5)
distributed...

Perhaps is obvious but I can't find the way to do it..

Thank you very much!

Marc

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


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





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

Re: [R] More Columns than column names Error

[Keith Jewell]

Carl is right.

Going to the nabble post and looking in the source data file 
 I see the 
headings row has 'Material' tab 'Weight...' tab 'Percent'.
Each of the data rows has 1 tab character between the 'Material' and 
'Weight' columns and 3 tab characters between the 'Weight...' and 
'Percent' columns.


On 22/10/2013 14:15, Carl Witthoft wrote:

What is the exact code you are using to try to load this file?
I strongly suspect the problem is a mixture of spaces and multiple tabs in
your text file.



--
View this message in context: 
http://r.789695.n4.nabble.com/More-Columns-than-column-names-Error-tp4678770p4678787.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] R licensing query

[Ed Keith]

Unfortunately this is how things work in the real world. I suspect the reason 
so many people keep getting in trouble for taking classified information home 
is because they can not get any work done on the office computer due to things 
like this. 

Many of the places I've worked have not permuted me to install Vim on my 
computer. I had to use MS Visual C++ editor or MS Notepad for all text editing.

Since I usually get payed by the hour, it just cost them more, and increases my 
income, but I still find it incredibly annoying.

   -EdK

Ed Keith
e_...@yahoo.com

Blog: edkeith.blogspot.com


--- On Thu, 6/17/10, Frank E Harrell Jr  wrote:

> From: Frank E Harrell Jr 
> Date: Thursday, June 17, 2010, 12:11 PM
> Pardon my english but you're working
> for idiots.  I'd look elsewhere if there are other
> options.  IT departments should be here to help get
> things done, not to help prevent good work from being done.
> 
> Frank
> 
> On 06/17/2010 04:28 AM, McAllister, Gina wrote:
> > I have recently started a new job at an NHS hospital
> in Scotland.  Since
> > I took up this post 6 months ago I have had an ongoing
> dispute with the
> > IT secutiry dept. who refuse to install R on my
> computer.  I previously
> > worked in another branch of the NHS where R was widely
> used and yet
> > there is nothing I can say which will persuade the IT
> dept here to even
> > visit the website!  With some help from our head
> of department, they
> > have now agreed to install R but only if they receive
> an email from 'R'
> > ensuring that it is licensed for commercial use, is
> compaitable with
> > Windows XP and will not affect the networked computer
> system here.  My
> > only other option for data anlaysis is Excel, we have
> no money for
> > S-plus or any other stats programme.  Can anyone
> suggest anything or
> > send me a suitable email?
> > 
> > Many thanks,
> > Georgina
> > 
> > 
> 
> -- Frank E Harrell Jr   Professor and
> Chairman        School of Medicine
>                
>      Department of
> Biostatistics   Vanderbilt University
> 
> __
> R-help@r-project.org
> mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Using if statement on function

[Keith Jewell]

You could also consider
isTRUE(all.equal(FUN, mean))

> isTRUE(all.equal(mean, mean))
[1] TRUE
> isTRUE(all.equal(mean, median))
[1] FALSE

HTH

Keith J
"Patrick Burns"  wrote in message 
news:4c28777f.1040...@pburns.seanet.com...
> If I understand the problem properly,
> you want something like this:
>
> function(FUN, ...)
> {
> FunName <- deparse(substitute(FUN))
> if(FunName == "mean") {
> ...
> } else if(FunName == "median") {
> ...
> }
> }
>
> Using 'switch' is an alternative to 'if'.
>
>
> On 28/06/2010 10:50, Etienne Stockhausen wrote:
>> Hello everybody,
>>
>> I'm trying to use a if-statment on a function. For a better
>> understanding I want to present a small example:
>>
>> FUN=mean # could also be median,sd or any other function
>> if (FUN == mean)
>> plot(...)
>> if (FUN == median)
>> plot(...) ...
>>
>> This doesn't work, because FUN is a function. I've already tried to
>> coerce the type of FUN with as.character( ), but that's also not
>> possible. I'm stuck with this task and it is absolutely necessary to
>> give FUN the class of a function.
>> I'm looking forward for any hints, clues or solutions for my problem.
>>
>> So much thanks in advance
>>
>> Etienne
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
> -- 
> Patrick Burns
> pbu...@pburns.seanet.com
> http://www.burns-stat.com
> (home of 'Some hints for the R beginner'
> and 'The R Inferno')
>

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

Re: [R] F# vs. R

[Ed Keith]

It's been a long time since I used Fortran, and I have only dabbled in F#, but 
I do not think translating Fortran (or R) to F# will be easy. F# is basicly a 
functional language (like ML) and a very differant mind set than Fortran (or R).

   -EdK

Ed Keith
e_...@yahoo.com

Blog: edkeith.blogspot.com


--- On Thu, 7/8/10, rkevinbur...@charter.net  wrote:

> From: rkevinbur...@charter.net 
> Subject: Re: [R] F# vs. R
> To: r-help@r-project.org, "Patrick Burns" , 
> serg...@gmail.com
> Date: Thursday, July 8, 2010, 10:16 AM
> True, porting old C and Fortran code
> to C# or F# would be a pain and probably riddled with errors
> but it is not too soon to start looking to see if there is a
> better way. There have been numerous ports of LAPACK, BLAS,
> etc. to C#. Maybe they could be leveraged.
> 
> Maybe just allowing packages to be wrtten in C# or F# would
> be helpful. And remember there is Mono.
> 
> Just my 2 cents.
> 
>  Patrick Burns 
> wrote: 
> > I'd like to hear answers to this as well.
> > A language doesn't have to be a complete
> > replacement to be useful.
> > 
> > F# seems to have some nice features.
> > 
> > Pat
> > 
> > On 07/07/2010 17:54, Sergey Goriatchev wrote:
> > > Hello, Marc
> > >
> > > No, I do not want to validate Cox PH. :-)
> > > I do use R daily, though right now I do not use
> the statistical part that much.
> > >
> > > I just generally wonder if any R-user tried F#
> and his/her opinions.
> > >
> > > Regards,
> > > Sergey
> > >
> > >
> > > On Wed, Jul 7, 2010 at 17:56, Marc Schwartz 
> wrote:
> > >> On Jul 7, 2010, at 10:31 AM, Sergey
> Goriatchev wrote:
> > >>
> > >>> Hello, everyone
> > >>>
> > >>> F# is now public. Compiled code should
> run  faster than R.
> > >>>
> > >>> Anyone has opinion on F# vs. R? Just
> curious
> > >>>
> > >>> Best,
> > >>> S
> > >>
> > >>
> > >> The key time critical parts of R are written
> in compiled C and FORTRAN.
> > >>
> > >> Of course, if you want to take the time to
> code and validate a Cox PH or mixed effects model in F# and
> then run them against R's coxph() or lme()/lmer() functions
> to test the timing, feel free...  :-)
> > >>
> > >> So unless there is a pre-existing library of
> statistical and related functionality for F#, perhaps you
> need to reconsider your query.
> > >>
> > >> Regards,
> > >>
> > >> Marc Schwartz
> > >>
> > >>
> > >
> > >
> > >
> > 
> > -- 
> > Patrick Burns
> > pbu...@pburns.seanet.com
> > http://www.burns-stat.com
> > (home of 'Some hints for the R beginner'
> > and 'The R Inferno')
> > 
> > __
> > R-help@r-project.org
> mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained,
> reproducible code.
> 
> __
> R-help@r-project.org
> mailing list
> https://stat.ethz.ch/mailman/listinfo/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] nls with some coefficients fixed

[Keith Jewell]

I'm using nls to fit a variety of different models. Here I use SSgompertz as 
an example.

I want the ability to fix one (or more) of the coefficients that would 
normally be optimised (e.g. fix b3=0.8).

Examples; based on and using data from example(SSgompertz)
#-
# vanilla call to nls, no coefficients fixed, works fine
nls(density ~ SSgompertz(log(conc), Asym, b2, b3), data = DNase.1)

#-
# naive attempt to fix one parameter. I don't believe the error message
nls(density ~ SSgompertz(log(conc), Asym, b2, b3=0.8), data = DNase.1)
# Error in nlsModel(formula, mf, start, wts) :
#  singular gradient matrix at initial parameter estimates

#
# fix parameter and use explicit start values works fine
nls(density ~ SSgompertz(log(conc), Asym, b2, b3=0.8), data = DNase.1, 
start=list(Asym=3, b2=2))

#
# getInitial returns third coeff, treating value as a name
getInitial(density ~ SSgompertz(log(conc), Asym=a, b2=b, b3=0.8), data = 
DNase.1)
a b   0.8
4.608 2.2713391 0.7164666


I guess the best approach in principle is to change the initial attribute so 
it only estimates and returns values of those coefficients for which 
is.name(mCall[coef]) == TRUE. BUT,that means handling all combinations of 
fixed coefficients so it's a bit of work even for one model. I'm working 
with a large number of models and really don't want to face all that work!

An alternative would be for nls to recognise fixed coefficients and omit 
their values from the start list; the other coefficients might not be 
optimal for the values of the fixed, but it does seem to work and wouldn't 
need all those SSxxx altering. BUT I don't want to get into nls internals.

Any suggestions?

Thanks in advance,

Keith J

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