Re: [R] Help with writing data to csv

2012-05-13 Thread John
On Sat, 12 May 2012 13:33:19 -0700 (PDT)
DL  wrote:

> I am trying to write data to csv but I am having issues with the
> separations. 
> 
> Basically I have some results that I get with R that I copied and
> pasted into word and then saved as .txt 
> 
> I want to write the results to csv because it's easier to make tables
> in word (all I would have to do is copy and paste into a table,
> instead of typing everything out). 
> 
...
I'm not sure what all you are doing wrong, but the problems are limited
to R.  You don't seem understand, Word, Excel, or how to use csv files
either.

Since you don't actually provide an example what "practice" is, it is
not clear just what the problem might be.  If perchance "Practice" is
the result of running a routine on data and capturing the result, then
you might want to try ?sink(), which, properly instructed, will write
the "results" to a file.  If that is the case, you might also want to
consider what fonts you are using in Word.

However, without providing the list with an example of Practice, we're
left to guess precisely what kind of help you really need.

JWD

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


Re: [R] round up/down to nearest hundredth

2012-05-13 Thread Petr Savicky
On Sat, May 12, 2012 at 01:38:16PM -0700, chuck.01 wrote:
> Hi, 
> How can I round up and down to the nearest hundredth
> 
> example:
> 
> x <- 0.18
> how do I round down to 0.15?
> how do I round down to 0.20?

Hi.

Rounding to 0.15 is not to a nearest hundredth, but to a nearest
multiple of 0.05. Rounding down and up to a nearest multiple
of a unit u can be done as follows

  x <- 0.18
  u <- 0.05
  floor(x/u)*u   # [1] 0.15
  ceiling(x/u)*u # [1] 0.2

Hope this helps.

Petr Savicky.

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


Re: [R] memory growth with rbind and lapply

2012-05-13 Thread Duncan Murdoch

On 12-05-12 10:06 PM, Jack Tanner wrote:

This version of my code makes the R process consume unreasonable amounts of RAM:

   datf<- rbind(lapply(mylist, function(item) {
   with(item, data.frame(col1, col2, col3))
 }))

This version works fine:

   datf<- lapply(mylist, function(item) {
   with(item, data.frame(col1, col2, col3))
 })
   datf<- do.call(rbind, datf)

Is this to be expected?


They are different, so I'd expect differences:

> x <- list(data.frame(a=1), data.frame(a=2),data.frame(a=3))
> rbind(x)
  [,1]   [,2]   [,3]
x List,1 List,1 List,1
> do.call(rbind, x)
  a
1 1
2 2
3 3

The first doesn't really make sense, so the fact that it takes a lot of 
memory may be R trying too hard to do something sensible.


Duncan Murdoch

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


Re: [R] Plotmath bug or my misunderstanding?

2012-05-13 Thread Duncan Murdoch

On 12-05-13 12:27 AM, David Winsemius wrote:


On May 12, 2012, at 6:27 PM, Peter Ehlers wrote:


Hi Bert,

I think the 'cex=' argument applies to the outermost 'atop()'.
It then applies that size specification to each of the two
components of the atop(a,b). If one of the components is itself
another atop(a,b), then the individual parts are sized downward
to produce the required cex for the unit.

As for the 'cex=1:2' specification, only the first value is used.


Isn't there also a 0.67 factor being applied to whatever was the
"outside" size prior to being used for the atop arguments? So the
inner atop arguments get that applied twice if I understand correctly.

text(1,1,labels=expression(atop(
 atop("top level one", "top level two"),
 "another level")
   ) )


I did notice that the cr-linefeed (perhaps unintentionally being
inserted by a mail client) seemed to be "working" at least on my Mac :

plot(1,1)
text(1,1,labels=expression("another
level")
   )

Whereas it is well known that "\n" will appear as two characters, \n,
and not be interpreted as a special.


I think that only works when the newline (which could equally be \n) is 
within a string; if you want a symbol over another symbol you need the atop.


I.e.

plot(1,1)
text(1,1,labels=expression("another\nlevel"))

is the same as what you typed, but

plot(1,1)
text(1,1,labels=expression(paste("another","\n", "level")))

won't do anything useful.

Duncan Murdoch

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


[R] plot.ts error: cannot plot more than 10 series

2012-05-13 Thread Suhaila Haji Mohd Hussin

Hello.
I'm doing Principal Component Analysis. So I'm trying to plot the new data as 
time series like this:
plot.ts(pr2$x)
But it gives me error: Error in plotts(x = x, y = y, plot.type = plot.type, 
xy.labels = xy.labels,  :  cannot plot more than 10 series as "multiple"
Here's my data dput it's really long so I just put in some essential 
information only for you to see:
structure(list(sdev = c(163.444257295666, 104.206376841326, 97.9593394161277, 
85.13008012602, 78.672940030061, 73.7044913582009, 68.9919940515061, 
66.4914391167912, 63.5364862493486, 53.0067010985736, 50.3850160673695, 
41.0505058921616, 36.6847894205371, 34.4922997539369, 32.8160496561452), 
rotation = structure(c(-0.341624025780871, -0.465361625382171, 
-0.254733743609541, -0.501601736275946, -0.32909803586058, 0.10346905261603, 
0.0826920791661065, 0.04180674828032, 0.0317514777236264, -0.041673145049323, 
0.00552096153263742, 0.302880866160923, 0.0611515661970536, -0.19666860832608, 
-0.0512002313448844, 0.0532192307778435, -0.279691445123637, 
-0.0330375822617488, 0.0296862808186525, 0.00754185914289494, 
-0.0483037204914515, -0.145633004945105, 0.0663131380872606, 
...0.018728860610997, -0.867607192081349, 
-0.00746836023138116, -0.197845760644588, 0.354006015245109, 
-0.0234568336704251, 0.17996661095123, 0.0910352921083386, 0.118413710329641), 
.D!
 im = c(20L, 15L), .Dimnames = list(c("er", "pr", "ar", "l1", "l2", "b1", 
"b2", "h1", "h2", "h3", "p1", "p2", "o1", "o2", "o3", "o4", "o5", "o6", 
"o7", "o8"), c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", 
"PC9", "PC10", "PC11", "PC12", "PC13", "PC14", "PC15"))), center = 
structure(c(94.4986072423398, 94.9637883008357, 71.1838440111421, 
179.793871866295, 180.57938718663, 12.6490250696379, 10.5292479108635, 
10.8690807799443, 38.2479108635098, 137.543175487465, 0.612813370473538, 
45.150417827298, 10.5459610027855, 91.3259052924791, 116.885793871866, 
59.0055710306407, 172.598885793872, 14.3676880222841, 5.0974930362117, 
3.8857938718663), .Names = c("er", "pr", "ar", "l1", "l2", "b1", "b2", 
"h1", "h2", "h3", "p1", "p2", "o1", "o2", "o3", "o4", "o5", "o6", "o7", 
"o8")), scale = FALSE, x = structure(c(-38.7640736739586, 
-85.3987965131485, 279.953119866405, 121.49850702536, 348.387859356463, 
-149.756720052989, 114.1!
 99587790773, -258.129158201887, 245.434330475019, -74.518625315929
, 114.858457403311, 246.989390591038, 264.221327740388, -59.7706080070487, 
-66.9519268984565, -202.267758377027, -61.2743296889704, -63.3534745591336, 
.
403375504553, 41.2280423622485, 78.1688899057297, -17.9229807680497, 
1.892930692659, -0.121464277679548, -10.9161618391865, 9.77739319210259, 
-30.4117535843182), .Dim = c(359L, 15L), .Dimnames = list(c("9", "12", 
"13", "14", "18", "22", "24", "29", "30", "31", "32", "33", "34", "44", 
"49", "53", "54", "58", "59", "60", "62", "66", "69", "70", "75", "76", 
"80", "87", "88", "93", "94", "96", "100", "101", "102", "106", "114", 
"115", "116", "126", "128", "129", "131", "133", "135", "137", "139", 
"141", "142", "144", "147", "151", "154", "157", "158", "160", "161", 
"164", "169", "170", "171", "175", "179", "180", "181", "186", "193", 
"195", "200", "205", "206", "209", "218", "220", "221", "222", "224", 
"225", "226", "228", "231", "232", "233", "235", "239", "243", "244", 
"253", "256", "260", "264", "265", "266", "268", "269", "272", "275", 
"280", .   "1025", 
"1032", "1033"!
 , "1035", "1036", "1038", "1039", "1041", "1042", "1046", "1047", "1048", 
"1051", "1055", "1059", "1060", "1064", "1072", "1076"), c("PC1", "PC2", 
"PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", "PC10", "PC11", "PC12", 
"PC13", "PC14", "PC15", .Names = c("sdev", "rotation", "center", 
"scale", "x"), class = "prcomp")
Please help,Suhaila   
[[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] x axis dates

2012-05-13 Thread Jim Lemon

On 05/13/2012 03:49 AM, pip wrote:

The following data spans out to 2012 and when graphing out the data - R
defaults dates on x axis only shows the years eg 200520072009

My question is how can I get the graph to show more dates on the axis - eg
every qtr or 6 months or every month etc

PRODUCT,DATE,AMOUNT_1,AMOUNT_2
A,3/1/2005,126647400,121942966.2
A,4/1/2005,121942966.2,122244045
A,5/1/2005,122244045,117227480.04
A,6/1/2005,117227480.04,117913720.32
A,7/1/2005,117913720.32,116767296.


Hi Phil,
Try the staxlab function in plotrix.

Jim

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


[R] determine size (width and height) of a graphics file via R - how?

2012-05-13 Thread Mark Heckmann
Hi,

is there a way to determine the size (width, height) of a graphics file saved 
on my hard disk, e.g. a .bmp, via R.
What I want is basically the same information on the dimensions of the graphic 
file that I get from my file browser.  

Thanks
Mark

PS. Why: I use the R2PPT and I need to determine the size of the original 
graphic before adding it to a slide.


Mark Heckmann
Blog: www.markheckmann.de
R-Blog: http://ryouready.wordpress.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] how to calculate risk allele or score allele

2012-05-13 Thread ramakanth reddy
Hello,

In a case control study how to calculate the risk allele or score allele.

Regards

GRR

[[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] categorical analysis - grouping rows

2012-05-13 Thread Rui Barradas


Hello,

If the output of ftable is kept in an object named 'ft', try this.


ftperc <- t(apply(ft, 1, function(x)
if(sum(x)) round(100*x/sum(x)) else rep(0, length(x
attributes(ftperc) <- attributes(ft)
ftperc


(And if you need it to process several datasets, a function with 'ft' as 
the argument and the code above as its body is straightforward to write.)


Hope this helps,

Rui Barradas

Em 13-05-2012 11:00, r-help-requ...@r-project.org escreveu:

Date: Sat, 12 May 2012 14:00:03 -0700 (PDT)
From: cassiorx
To:r-help@r-project.org
Subject: [R] categorical analysis - grouping rows
Message-ID:<1336856403964-4629503.p...@n4.nabble.com>
Content-Type: text/plain; charset=us-ascii

I apologize up front if this has been covered elsewhere - but I can't find
any such question.

I have a data set that contains academic data: term (i.e., semester),
student id, dept, class, success (1=Y, 0=N)

I want to look at dept by term to determine descriptive statistics for
success to failure ratios. The intent being to discover if there are
departments that contribute significantly to the Simpson Paradox, that is,
that make overall success/failure rates undependable.

It's easy to use ftable to get the counts for what I need (row names dept
and success, col name success.  So I get something that looks like this:

  Term  1st   2nd3rd4th5th
dept success
AAA  0   155240163286293
   1   424570349582429
AAB  055  64103   46109
  1   122117145112145
AAC  011 34 4 4
   119   12  23   117

How can I calculate percentages by dept so that I get

AAA 0 27  
  1 73  
AAB 0...

Part of my lack of understanding is that I don't see a way to get the dept
(by term) totals into a data structure that I can use to calculate the
percentages. I can write procedural code to do this but is there some r-way
that would be better?

--
View this message in 
context:http://r.789695.n4.nabble.com/categorical-analysis-grouping-rows-tp4629503.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] Replacing cretin value in a file

2012-05-13 Thread Jonsson
Dear All,

I am trying to replace a value of 528.8933 to - in my file

t<- file("C:\\Users\\Amin\\Desktop\\1999n_Resample11.img", "rb")
e=readBin(t, double(), size=4,n=720*360, signed=TRUE)
e[e != -] <- e[e != -]*0.0099 + 477.65 -273.15
##worked well and values heve been calculated
e[e == 528.8933] <- -## no changes made to
528.8933

any suggestions please


--
View this message in context: 
http://r.789695.n4.nabble.com/Replacing-cretin-value-in-a-file-tp4629840.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] Replacing cretin value in a file

2012-05-13 Thread Jim Holtman
FAQ 7.31

Sent from my iPad

On May 13, 2012, at 7:52, Jonsson  wrote:

> Dear All,
> 
> I am trying to replace a value of 528.8933 to - in my file
> 
> t<- file("C:\\Users\\Amin\\Desktop\\1999n_Resample11.img", "rb")
> e=readBin(t, double(), size=4,n=720*360, signed=TRUE)
> e[e != -] <- e[e != -]*0.0099 + 477.65 -273.15
> ##worked well and values heve been calculated
> e[e == 528.8933] <- -## no changes made to
> 528.8933
> 
> any suggestions please
> 
> 
> --
> View this message in context: 
> http://r.789695.n4.nabble.com/Replacing-cretin-value-in-a-file-tp4629840.html
> Sent from the R help mailing list archive at Nabble.com.
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

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


Re: [R] plot.ts error: cannot plot more than 10 series

2012-05-13 Thread Gabor Grothendieck
On Sun, May 13, 2012 at 5:02 AM, Suhaila Haji Mohd Hussin
 wrote:
>
> Hello.
> I'm doing Principal Component Analysis. So I'm trying to plot the new data as 
> time series like this:
> plot.ts(pr2$x)
> But it gives me error: Error in plotts(x = x, y = y, plot.type = plot.type, 
> xy.labels = xy.labels,  :  cannot plot more than 10 series as "multiple"
> Here's my data dput it's really long so I just put in some essential 
> information only for you to see:
> structure(list(sdev = c(163.444257295666, 104.206376841326, 97.9593394161277, 
> 85.13008012602, 78.672940030061, 73.7044913582009, 68.9919940515061, 
> 66.4914391167912, 63.5364862493486, 53.0067010985736, 50.3850160673695, 
> 41.0505058921616, 36.6847894205371, 34.4922997539369, 32.8160496561452), 
> rotation = structure(c(-0.341624025780871, -0.465361625382171, 
> -0.254733743609541, -0.501601736275946, -0.32909803586058, 0.10346905261603, 
> 0.0826920791661065, 0.04180674828032, 0.0317514777236264, -0.041673145049323, 
> 0.00552096153263742, 0.302880866160923, 0.0611515661970536, 
> -0.19666860832608, -0.0512002313448844, 0.0532192307778435, 
> -0.279691445123637, -0.0330375822617488, 0.0296862808186525, 
> 0.00754185914289494, -0.0483037204914515, -0.145633004945105, 
> 0.0663131380872606, ...0.018728860610997, 
> -0.867607192081349, -0.00746836023138116, -0.197845760644588, 
> 0.354006015245109, -0.0234568336704251, 0.17996661095123, 0.0910352921083386, 
> 0.118413710329641), .D!
>  im = c(20L, 15L), .Dimnames = list(    c("er", "pr", "ar", "l1", "l2", "b1", 
> "b2", "h1", "h2", "h3",     "p1", "p2", "o1", "o2", "o3", "o4", "o5", "o6", 
> "o7", "o8"    ), c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8",    
>  "PC9", "PC10", "PC11", "PC12", "PC13", "PC14", "PC15"))),     center = 
> structure(c(94.4986072423398, 94.9637883008357,     71.1838440111421, 
> 179.793871866295, 180.57938718663, 12.6490250696379,     10.5292479108635, 
> 10.8690807799443, 38.2479108635098, 137.543175487465,     0.612813370473538, 
> 45.150417827298, 10.5459610027855, 91.3259052924791,     116.885793871866, 
> 59.0055710306407, 172.598885793872, 14.3676880222841,     5.0974930362117, 
> 3.8857938718663), .Names = c("er", "pr",     "ar", "l1", "l2", "b1", "b2", 
> "h1", "h2", "h3", "p1", "p2",     "o1", "o2", "o3", "o4", "o5", "o6", "o7", 
> "o8")), scale = FALSE,     x = structure(c(-38.7640736739586, 
> -85.3987965131485, 279.953119866405,     121.49850702536, 348.387859356463, 
> -149.756720052989, 114.1!
>  99587790773,     -258.129158201887, 245.434330475019, -74.518625315929
> , 114.858457403311,     246.989390591038, 264.221327740388, 
> -59.7706080070487, -66.9519268984565,     -202.267758377027, 
> -61.2743296889704, -63.3534745591336, 
> .
> 403375504553, 41.2280423622485,     78.1688899057297, -17.9229807680497, 
> 1.892930692659, -0.121464277679548,     -10.9161618391865, 9.77739319210259, 
> -30.4117535843182), .Dim = c(359L,     15L), .Dimnames = list(c("9", "12", 
> "13", "14", "18", "22",     "24", "29", "30", "31", "32", "33", "34", "44", 
> "49", "53",     "54", "58", "59", "60", "62", "66", "69", "70", "75", "76",   
>   "80", "87", "88", "93", "94", "96", "100", "101", "102",     "106", "114", 
> "115", "116", "126", "128", "129", "131", "133",     "135", "137", "139", 
> "141", "142", "144", "147", "151", "154",     "157", "158", "160", "161", 
> "164", "169", "170", "171", "175",     "179", "180", "181", "186", "193", 
> "195", "200", "205", "206",     "209", "218", "220", "221", "222", "224", 
> "225", "226", "228",     "231", "232", "233", "235", "239", "243", "244", 
> "253", "256",     "260", "264", "265", "266", "268", "269", "272", "275", 
> "280", .   "1025", 
> "1032", "1033"!
>  , "1035", "1036", "1038", "1039", "1041",     "1042", "1046", "1047", 
> "1048", "1051", "1055", "1059", "1060",     "1064", "1072", "1076"), c("PC1", 
> "PC2", "PC3", "PC4", "PC5",     "PC6", "PC7", "PC8", "PC9", "PC10", "PC11", 
> "PC12", "PC13",     "PC14", "PC15", .Names = c("sdev", "rotation", 
> "center", "scale", "x"), class = "prcomp")
> Please help,Suhaila


The dput output seems to have gotten messed up in transmission so lets
use an artificial example:

set.seed(123)
DF <- as.data.frame(matrix(rnorm(600), nc = 12)) # 12 columns
p <- prcomp(DF)

# error - cannot plot more than 10
plot.ts(p$x)

# try it several ways with plot.zoo

library(zoo)
plot(as.zoo(p$x)) # ok

# or - all on one panel

k <- ncol(p$x)
plot(as.zoo(p$x), screen = 1, col = 1:k)

# or - 1st two on 1st panel, 2nd two on 2nd panel, etc.

plot(as.zoo(p$x), screen = seq(0, len = k) %/% 2, col = 1:2)

-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PL

Re: [R] Replacing cretin value in a file

2012-05-13 Thread Jonsson
Dear Jim,

You may forget to add your suggestions because what I see is just my post

--
View this message in context: 
http://r.789695.n4.nabble.com/Replacing-cretin-value-in-a-file-tp4629840p4629843.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] Plotmath bug or my misunderstanding?

2012-05-13 Thread Rui Barradas


Hello,

It seems like a "feature". When trying the example below with more 
atops, the two bottom most lines, and only those two lines, feature 
character expansion relative to the default size and relative to the 
line before last.


plot(1, type="n", xaxt='n', yaxt='n', ann=FALSE, xlim=c(-0.5, 6))
text(1, 1,labels=expression(atop(atop(sigma,"some text"),
"another level")), cex = 2)
text(3, 1, labels=expression(atop(atop(atop(sigma,"some text"),
"another level"), "third")), cex = 2)
text(5, 1, labels=expression(atop(atop(atop(atop(sigma,"some text"),
"another level"), "third"), "4th")), cex = 2)

Anyway, I couldn't find this behavior in the help pages.

Rui Barradas


Em 13-05-2012 11:00,  Bert Gunter escreveu:

Date: Sat, 12 May 2012 14:05:26 -0700
From: Bert Gunter
To:r-help@r-project.org
Subject: [R] Plotmath bug or my misunderstanding?
Message-ID:

Content-Type: text/plain; charset=ISO-8859-1

This is a followup to a recent post on using atop() to obtain
multiline expressions.

My reading of the plotmath docs makes it clear that issuing (in base
graphics) the specification

par(cex = 2)

doubles symbols and regular text in subsequent plotmath expressions.
However, it is unclear to me what specifying cex_within_  the
annotation function using plotmath should do, and the following seems
to want to have it both ways: ignore/obey )or maybe recycle?)

plot(1,type="n", xaxt='n', yaxt='n', ann=FALSE)
  text(1,1,labels=expression(atop(sigma,"some text")),cex = 2)
## obeys the cex specification in symbols and text

HOWEVER

plot(1,type="n", xaxt='n', yaxt='n', ann=FALSE)
  text(1,1,labels=expression(atop(atop(sigma,"some text"),"another
level")),cex = 2)
## ???

For even more fun, try:

plot(1,type="n", xaxt='n', yaxt='n', ann=FALSE)
  text(1,1,labels=expression(atop(atop(sigma,"some text"),"another
level")),cex = 1:2)
##

So I confess to being flummoxed. Enlightenment would be much appreciated.

Cheers,
Bert


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


Re: [R] Replacing cretin value in a file

2012-05-13 Thread David Winsemius


On May 13, 2012, at 8:12 AM, Jonsson wrote:


Dear Jim,

You may forget to add your suggestions because what I see is just my  
post


The suggestion was to read the referenced FAQ.


PLEASE do read the posting guide http://www.R-project.org/posting-guide.html


And while you are referencing the FAQ from the R help pages, you might  
also follow the link and read the Posting Guide. It will explain why  
including context is requested. Most people do not want to use Nabble  
and many of us somewhat annoyed that Nabble users seem to expect us to  
use it.


--

David Winsemius, MD
West Hartford, CT

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


Re: [R] Replacing cretin value in a file

2012-05-13 Thread Jonsson
Dear All, 

 I am trying to replace a value of 528.8933 to - in my file 

 t<- file("C:\\Users\\Amin\\Desktop\\1999n_Resample11.img", "rb") 
e=readBin(t, double(), size=4,n=720*360, signed=TRUE) 
 e[e != -] <- e[e != -]*0.0099 + 477.65 -273.15  

 This code given above would read a binary file . Then It would do some
clculations to  e.  That worked well. I opend the file again in order to
look at the resultsand I found weird numbers ( 528.8933).

I then used this  line  >e[e == 528.8933] <- -to replace
the value of 528.8933  by -. I got no errors but when I looked again at
the values, I found them as they were(528.8933 is still there) 
any idea on how to replace a value by another value?Thanks in advance


--
View this message in context: 
http://r.789695.n4.nabble.com/Replacing-cretin-value-in-a-file-tp4629840p4629846.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] Replacing cretin value in a file

2012-05-13 Thread David Winsemius


On May 13, 2012, at 9:08 AM, Jonsson wrote:


Dear All,

I am trying to replace a value of 528.8933 to - in my file

t<- file("C:\\Users\\Amin\\Desktop\\1999n_Resample11.img", "rb")
e=readBin(t, double(), size=4,n=720*360, signed=TRUE)
e[e != -] <- e[e != -]*0.0099 + 477.65 -273.15

This code given above would read a binary file . Then It would do some
clculations to  e.  That worked well. I opend the file again in  
order to

look at the resultsand I found weird numbers ( 528.8933).

I then used this  line  >e[e == 528.8933] <- -to  
replace
the value of 528.8933  by -. I got no errors but when I looked  
again at

the values, I found them as they were(528.8933 is still there)
any idea on how to replace a value by another value?Thanks in advance


From memory  ... since you have not yet learned to include context.  
You were advised to read the FAQ ... item 7.31.


Perhaps this example will help.

> x <- 7* round(528.8933/7, 6)

> x== 528.8933
[1] FALSE
> x
[1] 528.8933



PLEASE do read the posting guide http://www.R-project.org/posting-guide.html


--

David Winsemius, MD
West Hartford, CT

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


Re: [R] Replacing cretin value in a file

2012-05-13 Thread Berend Hasselman

On 13-05-2012, at 15:08, Jonsson wrote:

> Dear All, 
> 
> I am trying to replace a value of 528.8933 to - in my file 
> 
> t<- file("C:\\Users\\Amin\\Desktop\\1999n_Resample11.img", "rb") 
> e=readBin(t, double(), size=4,n=720*360, signed=TRUE) 
> e[e != -] <- e[e != -]*0.0099 + 477.65 -273.15  
> 
> This code given above would read a binary file . Then It would do some
> clculations to  e.  That worked well. I opend the file again in order to
> look at the resultsand I found weird numbers ( 528.8933).
> 
> I then used this  line  >e[e == 528.8933] <- -to replace
> the value of 528.8933  by -. I got no errors but when I looked again at
> the values, I found them as they were(528.8933 is still there) 
> any idea on how to replace a value by another value?Thanks in advance
> 

Please read the R FAQ section 7.31  "Why doesn't R think these numbers are 
equal?"

http://cran.r-project.org/faqs.html
http://cran.r-project.org/doc/FAQ/R-FAQ.html

Berend

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


Re: [R] Plotmath bug or my misunderstanding?

2012-05-13 Thread Bert Gunter
Peter/David:

1. For some reason, I didn't see Peter's reply on r-help.

2. To Peter: Aha!!
Let me play this back to you. In

 text(1,1,labels=expression(atop(atop(sigma,"some text"),"another
level")),cex = 2)

The (outer) whole atop() specification is allocated twice the amount
of space that would be required for the current font size, cex =1.
Then each of the top and bottom is allocated the amount of space
needed within this cex = 2 space. For the inner atop to fit within
that allocated amount of space, its top and bottom must get smaller,
of course. N'est-ce pas?

This makes sense to me!  My misunderstanding is thinking that cex
applied to the individual text components, not the entire expression.
Enlightenment much appreciated.

3. However, this now begs the question of how to keep all text and
symbols the same size, which is really what the OP really wanted (a
way to emulate multiple lines). I'll fool around with this to see what
I come up with now that I'm enlightened -- or see if you or David come
up with something clever.


Cheers,
Bert



On Sat, May 12, 2012 at 9:27 PM, David Winsemius  wrote:
>
> On May 12, 2012, at 6:27 PM, Peter Ehlers wrote:
>
>> Hi Bert,
>>
>> I think the 'cex=' argument applies to the outermost 'atop()'.
>> It then applies that size specification to each of the two
>> components of the atop(a,b). If one of the components is itself
>> another atop(a,b), then the individual parts are sized downward
>> to produce the required cex for the unit.
>>
>> As for the 'cex=1:2' specification, only the first value is used.
>
>
> Isn't there also a 0.67 factor being applied to whatever was the "outside"
> size prior to being used for the atop arguments? So the inner atop arguments
> get that applied twice if I understand correctly.
>
> text(1,1,labels=expression(atop(
>               atop("top level one", "top level two"),
>               "another level")
>     ) )
>
>
> I did notice that the cr-linefeed (perhaps unintentionally being inserted by
> a mail client) seemed to be "working" at least on my Mac :
>
> plot(1,1)
> text(1,1,labels=expression("another
> level")
>     )
>
> Whereas it is well known that "\n" will appear as two characters, \n, and
> not be interpreted as a special.
> --
> David
>>
>>
>> Cheers,
>> Peter Ehlers
>>
>> On 2012-05-12 14:05, Bert Gunter wrote:
>>>
>>> This is a followup to a recent post on using atop() to obtain
>>> multiline expressions.
>>>
>>> My reading of the plotmath docs makes it clear that issuing (in base
>>> graphics) the specification
>>>
>>> par(cex = 2)
>>>
>>> doubles symbols and regular text in subsequent plotmath expressions.
>>> However, it is unclear to me what specifying cex _within_ the
>>> annotation function using plotmath should do, and the following seems
>>> to want to have it both ways: ignore/obey )or maybe recycle?)
>>>
>>> plot(1,type="n", xaxt='n', yaxt='n', ann=FALSE)
>>>  text(1,1,labels=expression(atop(sigma,"some text")),cex = 2)
>>> ## obeys the cex specification in symbols and text
>>>
>>> HOWEVER
>>>
>>> plot(1,type="n", xaxt='n', yaxt='n', ann=FALSE)
>>>  text(1,1,labels=expression(atop(atop(sigma,"some text"),"another
>>> level")),cex = 2)
>>> ## ???
>>>
>>> For even more fun, try:
>>>
>>> plot(1,type="n", xaxt='n', yaxt='n', ann=FALSE)
>>>  text(1,1,labels=expression(atop(atop(sigma,"some text"),"another 
>>> level")),cex = 1:2)
>>> ##
>>>
>>> So I confess to being flummoxed. Enlightenment would be much appreciated.
>>>
>>> Cheers,
>>> Bert
>>>
>>>
>>>
>>>
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
>
> David Winsemius, MD
> West Hartford, CT
>



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

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

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


Re: [R] set specific contrasts using lapply

2012-05-13 Thread Rui Barradas

Hello,

The value change inside a function's body is lost on exit.
You have changed the value of a copy of the argument passed to the 
function. That copy exists in the function's environment only. And the 
original exists in another environment, the caller's, the parent frame. 
(.GlobalEnv?)

This is valid for contrasts or for any other change you make.
In the case of contrasts, the loop is better, you will always have to 
check wether the object is a factor or logical vector before assigning 
different contrasts levels to it, so you can not use


contrasts(df or list) <- lapply(whatever)

Hope this helps,

Rui Barradas

Em 12-05-2012 11:00, r-help-requ...@r-project.org escreveu:

Date: Fri, 11 May 2012 05:38:59 -0700 (PDT)
From: Frank Paetzold
To:r-help@r-project.org
Subject: [R] set specific contrasts using lapply
Message-ID:<1336739939629-4626289.p...@n4.nabble.com>
Content-Type: text/plain; charset=us-ascii

I have the following data set


>  data

A  B  X1  X2   Y
1 A1 B1 1.1 2.9 1.2
2 A1 B2 1.0 3.2 2.3
3 A2 B1 1.0 3.3 1.6
4 A2 B2 0.5 2.6 3.1


>  sapply(data, class)

 A BX1X2 Y
  "factor"  "factor" "numeric" "numeric" "numeric"

I'd like to set a specific type of contrasts to all the categorical factors
(here A and B).
In a loop, this can be done like this:

for (i in 1:length(data))
{
   if (class(data[[i]]) == 'factor')
   {
 contrasts(data[[i]], length(levels(data[[i]])))<-
contr.treatment(levels(data[[i]]), contrast=FALSE);
   }
}

Can i do this without the loop?

I tried to use the function

set.contrasts<- function(x)
{
   if (class(x) == 'factor')
   {
 contrasts(x, length(levels(x)))<- contr.treatment(levels(x),
contrast=FALSE);
   }
}

in

lapply(data, set.contrasts)

However, it doesn't work.

Any ideas?

Thank You


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


[R] heatmap

2012-05-13 Thread li li
Dear all,
   I have a data set which is 100 \times 100 matrix.
I made a heatmap of the data. The xlim and ylim
are both shown to be c(0,1). I want to make both the
xlim and ylim to be c(0,100). So the heatmap correspond to
the 100 \times 100 locations.
  Thank you.
Hannah

[[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] Plotmath bug or my misunderstanding?

2012-05-13 Thread David Winsemius


On May 13, 2012, at 10:43 AM, Bert Gunter wrote:


Peter/David:

1. For some reason, I didn't see Peter's reply on r-help.

2. To Peter: Aha!!
Let me play this back to you. In

text(1,1,labels=expression(atop(atop(sigma,"some text"),"another
level")),cex = 2)

The (outer) whole atop() specification is allocated twice the amount
of space that would be required for the current font size, cex =1.
Then each of the top and bottom is allocated the amount of space
needed within this cex = 2 space. For the inner atop to fit within
that allocated amount of space, its top and bottom must get smaller,
of course. N'est-ce pas?

This makes sense to me!  My misunderstanding is thinking that cex
applied to the individual text components, not the entire expression.
Enlightenment much appreciated.

3. However, this now begs the question of how to keep all text and
symbols the same size, which is really what the OP really wanted (a
way to emulate multiple lines). I'll fool around with this to see what
I come up with now that I'm enlightened -- or see if you or David come
up with something clever.


Any cleverness I exhibit on these pages is generally ascribable to my  
R-betters. When faced with a request like this, I just fire up Google  
and generally find an answer from Lumley, Grothendieck, Ligges, or in  
this case, Schwartz (from an rhelp post in 2009). This was the first  
hit that wasn't a link to R documentation with a search on "multiple  
lines expression plotmath r"


# values to pull from using bquote
 X= 2.3
Y= 5.6
txtList <- list("Line # 1",
 "Longer line #2",
  bquote(X == .(X)),
 "and",
  bquote(The~value~of~Y == .(Y)))

plot(1,type="n", xaxt='n', yaxt='n', ann=FALSE)
text(1.2, seq(1.2, 1.1, len=5) , labels=do.call(expression, txtList))

I think this is a fairly generally strategy but I'm certainly willing  
to take a crack at further modifications if desired.


--
david.




Cheers,
Bert



On Sat, May 12, 2012 at 9:27 PM, David Winsemius > wrote:


On May 12, 2012, at 6:27 PM, Peter Ehlers wrote:


Hi Bert,

I think the 'cex=' argument applies to the outermost 'atop()'.
It then applies that size specification to each of the two
components of the atop(a,b). If one of the components is itself
another atop(a,b), then the individual parts are sized downward
to produce the required cex for the unit.

As for the 'cex=1:2' specification, only the first value is used.



Isn't there also a 0.67 factor being applied to whatever was the  
"outside"
size prior to being used for the atop arguments? So the inner atop  
arguments

get that applied twice if I understand correctly.

text(1,1,labels=expression(atop(
  atop("top level one", "top level two"),
  "another level")
) )


I did notice that the cr-linefeed (perhaps unintentionally being  
inserted by

a mail client) seemed to be "working" at least on my Mac :

plot(1,1)
text(1,1,labels=expression("another
level")
)

Whereas it is well known that "\n" will appear as two characters,  
\n, and

not be interpreted as a special.
--
David



Cheers,
Peter Ehlers

On 2012-05-12 14:05, Bert Gunter wrote:


This is a followup to a recent post on using atop() to obtain
multiline expressions.

My reading of the plotmath docs makes it clear that issuing (in  
base

graphics) the specification

par(cex = 2)

doubles symbols and regular text in subsequent plotmath  
expressions.

However, it is unclear to me what specifying cex _within_ the
annotation function using plotmath should do, and the following  
seems

to want to have it both ways: ignore/obey )or maybe recycle?)

plot(1,type="n", xaxt='n', yaxt='n', ann=FALSE)
 text(1,1,labels=expression(atop(sigma,"some text")),cex = 2)
## obeys the cex specification in symbols and text

HOWEVER

plot(1,type="n", xaxt='n', yaxt='n', ann=FALSE)
 text(1,1,labels=expression(atop(atop(sigma,"some text"),"another
level")),cex = 2)
## ???

For even more fun, try:

plot(1,type="n", xaxt='n', yaxt='n', ann=FALSE)
 text(1,1,labels=expression(atop(atop(sigma,"some text"),"another  
level")),cex = 1:2)

##

So I confess to being flummoxed. Enlightenment would be much  
appreciated.


Cheers,
Bert






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



David Winsemius, MD
West Hartford, CT





--

Bert Gunter
Genentech Nonclinical Biostatistics

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


David Winsemius, MD
West Hartford, CT

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

Re: [R] Plotmath bug or my misunderstanding?

2012-05-13 Thread P Ehlers

Bert: inline

On 2012-05-13 7:43, Bert Gunter wrote:

Peter/David:

1. For some reason, I didn't see Peter's reply on r-help.

2. To Peter: Aha!!
Let me play this back to you. In

  text(1,1,labels=expression(atop(atop(sigma,"some text"),"another
level")),cex = 2)

The (outer) whole atop() specification is allocated twice the amount
of space that would be required for the current font size, cex =1.
Then each of the top and bottom is allocated the amount of space
needed within this cex = 2 space. For the inner atop to fit within
that allocated amount of space, its top and bottom must get smaller,
of course. N'est-ce pas?


Yes, that's my understanding. With your text() call above, the
size of "another level" is what text(x,y,lab="another text",cex=2)
will print.



This makes sense to me!  My misunderstanding is thinking that cex
applied to the individual text components, not the entire expression.
Enlightenment much appreciated.

3. However, this now begs the question of how to keep all text and
symbols the same size, which is really what the OP really wanted (a
way to emulate multiple lines). I'll fool around with this to see what
I come up with now that I'm enlightened -- or see if you or David come
up with something clever.


I don't think it's possible to nest atop()s and have all members
of the same cex.

Anyway, atop() is really the wrong tool for this sort of thing.
It can be a useful quick fix, but I usually prefer the 'multiple
lines' approach for its greater flexibility.

Peter Ehlers




Cheers,
Bert



On Sat, May 12, 2012 at 9:27 PM, David Winsemius  wrote:


On May 12, 2012, at 6:27 PM, Peter Ehlers wrote:


Hi Bert,

I think the 'cex=' argument applies to the outermost 'atop()'.
It then applies that size specification to each of the two
components of the atop(a,b). If one of the components is itself
another atop(a,b), then the individual parts are sized downward
to produce the required cex for the unit.

As for the 'cex=1:2' specification, only the first value is used.



Isn't there also a 0.67 factor being applied to whatever was the "outside"
size prior to being used for the atop arguments? So the inner atop arguments
get that applied twice if I understand correctly.

text(1,1,labels=expression(atop(
   atop("top level one", "top level two"),
   "another level")
 ) )


I did notice that the cr-linefeed (perhaps unintentionally being inserted by
a mail client) seemed to be "working" at least on my Mac :

plot(1,1)
text(1,1,labels=expression("another
level")
 )

Whereas it is well known that "\n" will appear as two characters, \n, and
not be interpreted as a special.
--
David



Cheers,
Peter Ehlers

On 2012-05-12 14:05, Bert Gunter wrote:


This is a followup to a recent post on using atop() to obtain
multiline expressions.

My reading of the plotmath docs makes it clear that issuing (in base
graphics) the specification

par(cex = 2)

doubles symbols and regular text in subsequent plotmath expressions.
However, it is unclear to me what specifying cex _within_ the
annotation function using plotmath should do, and the following seems
to want to have it both ways: ignore/obey )or maybe recycle?)

plot(1,type="n", xaxt='n', yaxt='n', ann=FALSE)
  text(1,1,labels=expression(atop(sigma,"some text")),cex = 2)
## obeys the cex specification in symbols and text

HOWEVER

plot(1,type="n", xaxt='n', yaxt='n', ann=FALSE)
  text(1,1,labels=expression(atop(atop(sigma,"some text"),"another
level")),cex = 2)
## ???

For even more fun, try:

plot(1,type="n", xaxt='n', yaxt='n', ann=FALSE)
  text(1,1,labels=expression(atop(atop(sigma,"some text"),"another level")),cex 
= 1:2)
##

So I confess to being flummoxed. Enlightenment would be much appreciated.

Cheers,
Bert






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



David Winsemius, MD
West Hartford, CT







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


Re: [R] heatmap

2012-05-13 Thread David Winsemius


On May 13, 2012, at 11:23 AM, li li wrote:


Dear all,
  I have a data set which is 100 \times 100 matrix.
I made a heatmap of the data. The xlim and ylim
are both shown to be c(0,1). I want to make both the
xlim and ylim to be c(0,100). So the heatmap correspond to
the 100 \times 100 locations.


So just suppress the x and y tick marks and then use axis() to label  
as desired.




 Thank you.
   Hannah

[[alternative HTML version deleted]]

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


David Winsemius, MD
West Hartford, CT

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


Re: [R] Replacing cretin value in a file

2012-05-13 Thread Ken Hutchison
Perhaps you are in the first circle of the R Inferno?
   Ken Hutchison



On May 13, 2012, at 9:31 AM, Berend Hasselman  wrote:

> 
> On 13-05-2012, at 15:08, Jonsson wrote:
> 
>> Dear All, 
>> 
>> I am trying to replace a value of 528.8933 to - in my file 
>> 
>> t<- file("C:\\Users\\Amin\\Desktop\\1999n_Resample11.img", "rb") 
>> e=readBin(t, double(), size=4,n=720*360, signed=TRUE) 
>> e[e != -] <- e[e != -]*0.0099 + 477.65 -273.15  
>> 
>> This code given above would read a binary file . Then It would do some
>> clculations to  e.  That worked well. I opend the file again in order to
>> look at the resultsand I found weird numbers ( 528.8933).
>> 
>> I then used this  line  >e[e == 528.8933] <- -to replace
>> the value of 528.8933  by -. I got no errors but when I looked again at
>> the values, I found them as they were(528.8933 is still there) 
>> any idea on how to replace a value by another value?Thanks in advance
>> 
> 
> Please read the R FAQ section 7.31  "Why doesn't R think these numbers are 
> equal?"
> 
> http://cran.r-project.org/faqs.html
> http://cran.r-project.org/doc/FAQ/R-FAQ.html
> 
> Berend
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

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


[R] Residuals from binomial glm

2012-05-13 Thread Petar Milin
Dear all!
Sometimes, when using plain lm(), if I have highly correlated numeric
predictors (covariates), I firstly ran lm() to get residuals of one from
the other. This way, they became uncorrelated (orthogonal), and I can
check if the residualized one contributes to prediction, over and above
the other. I do that if and only if there is "theoretical" sense/ground
for such step.
Now, I would need to run binomial glm, but I think that those residuals
may "trick" me. Logic of glm's various types of residuals tells me that
the "response" should be a match to those which I would get from simple
lm(). However, I do not want to wander blindly into this.
Please, does anyone have experience with using glm.residuals on the
right-hand side of main lm() model? What to use, how and why?
Any advice on literature?

Best,
Petar

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


Re: [R] plot.ts error: cannot plot more than 10 series

2012-05-13 Thread Suhaila Haji Mohd Hussin

Thank you!
Suhaila

> From: ggrothendi...@gmail.com
> Date: Sun, 13 May 2012 07:58:19 -0400
> Subject: Re: [R] plot.ts error: cannot plot more than 10 series
> To: bell_beaut...@hotmail.com
> CC: r-help@r-project.org
> 
> On Sun, May 13, 2012 at 5:02 AM, Suhaila Haji Mohd Hussin
>  wrote:
> >
> > Hello.
> > I'm doing Principal Component Analysis. So I'm trying to plot the new data 
> > as time series like this:
> > plot.ts(pr2$x)
> > But it gives me error: Error in plotts(x = x, y = y, plot.type = plot.type, 
> > xy.labels = xy.labels,  :  cannot plot more than 10 series as "multiple"
> > Here's my data dput it's really long so I just put in some essential 
> > information only for you to see:
> > structure(list(sdev = c(163.444257295666, 104.206376841326, 
> > 97.9593394161277, 85.13008012602, 78.672940030061, 73.7044913582009, 
> > 68.9919940515061, 66.4914391167912, 63.5364862493486, 53.0067010985736, 
> > 50.3850160673695, 41.0505058921616, 36.6847894205371, 34.4922997539369, 
> > 32.8160496561452), rotation = structure(c(-0.341624025780871, 
> > -0.465361625382171, -0.254733743609541, -0.501601736275946, 
> > -0.32909803586058, 0.10346905261603, 0.0826920791661065, 0.04180674828032, 
> > 0.0317514777236264, -0.041673145049323, 0.00552096153263742, 
> > 0.302880866160923, 0.0611515661970536, -0.19666860832608, 
> > -0.0512002313448844, 0.0532192307778435, -0.279691445123637, 
> > -0.0330375822617488, 0.0296862808186525, 0.00754185914289494, 
> > -0.0483037204914515, -0.145633004945105, 0.0663131380872606, 
> > ...0.018728860610997, -0.867607192081349, 
> > -0.00746836023138116, -0.197845760644588, 0.354006015245109, 
> > -0.0234568336704251, 0.17996661095123, 0.0910352921083386, 
> > 0.118413710329641)!
 , .D!
> >  im = c(20L, 15L), .Dimnames = list(c("er", "pr", "ar", "l1", "l2", 
> > "b1", "b2", "h1", "h2", "h3", "p1", "p2", "o1", "o2", "o3", "o4", "o5", 
> > "o6", "o7", "o8"), c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", 
> > "PC8", "PC9", "PC10", "PC11", "PC12", "PC13", "PC14", "PC15"))), 
> > center = structure(c(94.4986072423398, 94.9637883008357, 
> > 71.1838440111421, 179.793871866295, 180.57938718663, 12.6490250696379, 
> > 10.5292479108635, 10.8690807799443, 38.2479108635098, 137.543175487465, 
> > 0.612813370473538, 45.150417827298, 10.5459610027855, 91.3259052924791, 
> > 116.885793871866, 59.0055710306407, 172.598885793872, 14.3676880222841, 
> > 5.0974930362117, 3.8857938718663), .Names = c("er", "pr", "ar", "l1", 
> > "l2", "b1", "b2", "h1", "h2", "h3", "p1", "p2", "o1", "o2", "o3", "o4", 
> > "o5", "o6", "o7", "o8")), scale = FALSE, x = 
> > structure(c(-38.7640736739586, -85.3987965131485, 279.953119866405, 
> > 121.49850702536, 348.387859356463, -149.756720052989, 1!
 14.1!
> >  99587790773, -258.129158201887, 245.434330475019, -74.518625315929
> > , 114.858457403311, 246.989390591038, 264.221327740388, 
> > -59.7706080070487, -66.9519268984565, -202.267758377027, 
> > -61.2743296889704, -63.3534745591336, 
> > .
> > 403375504553, 41.2280423622485, 78.1688899057297, -17.9229807680497, 
> > 1.892930692659, -0.121464277679548, -10.9161618391865, 
> > 9.77739319210259, -30.4117535843182), .Dim = c(359L, 15L), .Dimnames = 
> > list(c("9", "12", "13", "14", "18", "22", "24", "29", "30", "31", "32", 
> > "33", "34", "44", "49", "53", "54", "58", "59", "60", "62", "66", "69", 
> > "70", "75", "76", "80", "87", "88", "93", "94", "96", "100", "101", 
> > "102", "106", "114", "115", "116", "126", "128", "129", "131", "133",   
> >   "135", "137", "139", "141", "142", "144", "147", "151", "154", "157", 
> > "158", "160", "161", "164", "169", "170", "171", "175", "179", "180", 
> > "181", "186", "193", "195", "200", "205", "206", "209", "218", "220", 
> > "221", "222", "224", "225", "226", "228", "231", "232", "233", "235", 
> > "239", "243", "244", "253", "256", "260", "264", "265", "266", "268", 
> > "269", "272", "275", "280", 
> > .   "1025", "1032", 
> > "1!
 033"!
> >  , "1035", "1036", "1038", "1039", "1041", "1042", "1046", "1047", 
> > "1048", "1051", "1055", "1059", "1060", "1064", "1072", "1076"), 
> > c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", 
> > "PC10", "PC11", "PC12", "PC13", "PC14", "PC15", .Names = c("sdev", 
> > "rotation", "center", "scale", "x"), class = "prcomp")
> > Please help,Suhaila
> 
> 
> The dput output seems to have gotten messed up in transmission so lets
> use an artificial example:
> 
> set.seed(123)
> DF <- as.data.frame(matrix(rnorm(600), nc = 12)) # 12 columns
> p <- prcomp(DF)
> 
> # error - cannot plot more than 10
> plot.ts(p$x)
> 
> # try it several ways with plot.zoo
> 
> library(zoo)
> plot(as.zoo(p$x)) # ok
> 
> # or - all on one panel
> 
> k <- ncol(p$x)
> p

[R] Couldn't plot all PCA new attributes on K-Means Clustering

2012-05-13 Thread Suhaila Haji Mohd Hussin

Hello.
I'm trying to plot clusters for all possible combination of attributes.
So I typed in 
plot(pr2NewMedicalData, col = kmpr2NewMedicalData$cluster)
The problem is it only plots PC2 axis against PC1 axis. There should be PC1 to 
PC10. I tried it on normal data (not PCA) and it plots all possible combination 
of attributes.
Here's the data before applying K-means clustering. It's too long so I just 
show you some essential information;
structure(c(8.18794213444232, 108.09939741361, -129.331704413972, 
245.596491781615, 284.97336181805, 46.2393277833597, 102.670149102092, 
-56.0868873906073, 78.7111917257909, 
17.2805919727568, 8.21955994418973, -52.8841201560606), .Dim = c(1076L, 10L), 
.Dimnames = list(NULL, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", 
"PC8", "PC9", "PC10")))
And here's the data after applying K-means clustering. It's too long as well so 
I just show you some essential information only;
structure(list(cluster = c(2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 
1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 
.1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 
2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L), 
centers = structure(c(202.683498408051, -91.234890118988, -6.32822645738145, 
2.84855476653021, 7.14162520784024, -3.21469382670974, 9.13147239059093, 
-4.11039323242235, -2.71400953442991, 1.22167005997248, 2.69685863804327, 
-1.21394984515696, -1.19187335000922, 0.536503637335695, 0.509142465076951, 
-0.229182726867495, 4.82054598771028, -2.16989536374021, 4.66753734226828, 
-2.10102085218006), .Dim = c(2L, 10L), .Dimnames = list(c("1", "2"), c("PC1", 
"PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", "PC10"))), totss = 
86060500.4826787, withinss = c(21892651.1737651, 44156468.7023972), 
tot.withinss = 66049119.8761623, betweenss = 20011380.6065164, size = 
c(334L, 742L)), .Names = c("cluster", !
 "centers", "totss", "withinss", "tot.withinss", "betweenss", "size"), class = 
"kmeans")
Please help,Suhaila   
[[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] Plotmath bug or my misunderstanding?

2012-05-13 Thread Bert Gunter
Peter/David:

Yes, I think this is the general approach for multiple lines in
plotmath expressions. For future reference, a simple canonical example
is:

exl <- list(quote(sigma), "tau",quote(beta))
plot(0:1,0:1,type="n",axes=FALSE)
text(do.call(expression,exl),x=.5,y=seq(.4,.6,by=.1),cex=2)

The points to note are:

1. The plotmath expressions, but not quoted character strings,  in the
list have to be protected from evaluation by quote() (or bquote if
variable values are to be substituted in them);

2. The individual components/lines have to be positioned manually via
x and y; there is no newline ("\n") capability, since plotmath is
"drawing" the expressions.

3. atop()  is NOT the right way to do this.

Thanks to you both for clarifying this for me.

Best,
Bert

On Sun, May 13, 2012 at 9:45 AM, David Winsemius  wrote:
>
> On May 13, 2012, at 10:43 AM, Bert Gunter wrote:
>
>> Peter/David:
>>
>> 1. For some reason, I didn't see Peter's reply on r-help.
>>
>> 2. To Peter: Aha!!
>> Let me play this back to you. In
>>
>> text(1,1,labels=expression(atop(atop(sigma,"some text"),"another
>> level")),cex = 2)
>>
>> The (outer) whole atop() specification is allocated twice the amount
>> of space that would be required for the current font size, cex =1.
>> Then each of the top and bottom is allocated the amount of space
>> needed within this cex = 2 space. For the inner atop to fit within
>> that allocated amount of space, its top and bottom must get smaller,
>> of course. N'est-ce pas?
>>
>> This makes sense to me!  My misunderstanding is thinking that cex
>> applied to the individual text components, not the entire expression.
>> Enlightenment much appreciated.
>>
>> 3. However, this now begs the question of how to keep all text and
>> symbols the same size, which is really what the OP really wanted (a
>> way to emulate multiple lines). I'll fool around with this to see what
>> I come up with now that I'm enlightened -- or see if you or David come
>> up with something clever.
>
>
> Any cleverness I exhibit on these pages is generally ascribable to my
> R-betters. When faced with a request like this, I just fire up Google and
> generally find an answer from Lumley, Grothendieck, Ligges, or in this case,
> Schwartz (from an rhelp post in 2009). This was the first hit that wasn't a
> link to R documentation with a search on "multiple lines expression plotmath
> r"
>
> # values to pull from using bquote
>  X= 2.3
> Y= 5.6
> txtList <- list("Line # 1",
>     "Longer line #2",
>      bquote(X == .(X)),
>     "and",
>      bquote(The~value~of~Y == .(Y)))
>
> plot(1,type="n", xaxt='n', yaxt='n', ann=FALSE)
> text(1.2, seq(1.2, 1.1, len=5) , labels=do.call(expression, txtList))
>
> I think this is a fairly generally strategy but I'm certainly willing to
> take a crack at further modifications if desired.
>
> --
> david.
>
>>
>>
>> Cheers,
>> Bert
>>
>>
>>
>> On Sat, May 12, 2012 at 9:27 PM, David Winsemius 
>> wrote:
>>>
>>>
>>> On May 12, 2012, at 6:27 PM, Peter Ehlers wrote:
>>>
 Hi Bert,

 I think the 'cex=' argument applies to the outermost 'atop()'.
 It then applies that size specification to each of the two
 components of the atop(a,b). If one of the components is itself
 another atop(a,b), then the individual parts are sized downward
 to produce the required cex for the unit.

 As for the 'cex=1:2' specification, only the first value is used.
>>>
>>>
>>>
>>> Isn't there also a 0.67 factor being applied to whatever was the
>>> "outside"
>>> size prior to being used for the atop arguments? So the inner atop
>>> arguments
>>> get that applied twice if I understand correctly.
>>>
>>> text(1,1,labels=expression(atop(
>>>              atop("top level one", "top level two"),
>>>              "another level")
>>>    ) )
>>>
>>>
>>> I did notice that the cr-linefeed (perhaps unintentionally being inserted
>>> by
>>> a mail client) seemed to be "working" at least on my Mac :
>>>
>>> plot(1,1)
>>> text(1,1,labels=expression("another
>>> level")
>>>    )
>>>
>>> Whereas it is well known that "\n" will appear as two characters, \n, and
>>> not be interpreted as a special.
>>> --
>>> David



 Cheers,
 Peter Ehlers

 On 2012-05-12 14:05, Bert Gunter wrote:
>
>
> This is a followup to a recent post on using atop() to obtain
> multiline expressions.
>
> My reading of the plotmath docs makes it clear that issuing (in base
> graphics) the specification
>
> par(cex = 2)
>
> doubles symbols and regular text in subsequent plotmath expressions.
> However, it is unclear to me what specifying cex _within_ the
> annotation function using plotmath should do, and the following seems
> to want to have it both ways: ignore/obey )or maybe recycle?)
>
> plot(1,type="n", xaxt='n', yaxt='n', ann=FALSE)
>  text(1,1,labels=expression(atop(sigma,"some text")),cex = 2)
> ## obeys the cex specification in symbols and t

Re: [R] Couldn't plot all PCA new attributes on K-Means Clustering

2012-05-13 Thread Suhaila Haji Mohd Hussin

Please nevermind as I've just solved this by using scatterplot matrices.
Suhaila

> From: bell_beaut...@hotmail.com
> To: r-help@r-project.org
> Date: Mon, 14 May 2012 06:28:06 +1200
> Subject: [R] Couldn't plot all PCA new attributes on K-Means Clustering
> 
> 
> Hello.
> I'm trying to plot clusters for all possible combination of attributes.
> So I typed in 
> plot(pr2NewMedicalData, col = kmpr2NewMedicalData$cluster)
> The problem is it only plots PC2 axis against PC1 axis. There should be PC1 
> to PC10. I tried it on normal data (not PCA) and it plots all possible 
> combination of attributes.
> Here's the data before applying K-means clustering. It's too long so I just 
> show you some essential information;
> structure(c(8.18794213444232, 108.09939741361, -129.331704413972, 
> 245.596491781615, 284.97336181805, 46.2393277833597, 102.670149102092, 
> -56.0868873906073, 78.7111917257909, 
> 17.2805919727568, 8.21955994418973, -52.8841201560606), .Dim = c(1076L, 10L), 
> .Dimnames = list(NULL, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", 
> "PC8", "PC9", "PC10")))
> And here's the data after applying K-means clustering. It's too long as well 
> so I just show you some essential information only;
> structure(list(cluster = c(2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 
> 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 
> .1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 
> 2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L), 
> centers = structure(c(202.683498408051, -91.234890118988, -6.32822645738145, 
> 2.84855476653021, 7.14162520784024, -3.21469382670974, 9.13147239059093, 
> -4.11039323242235, -2.71400953442991, 1.22167005997248, 2.69685863804327, 
> -1.21394984515696, -1.19187335000922, 0.536503637335695, 0.509142465076951, 
> -0.229182726867495, 4.82054598771028, -2.16989536374021, 4.66753734226828, 
> -2.10102085218006), .Dim = c(2L, 10L), .Dimnames = list(c("1", "2"), c("PC1", 
> "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", "PC10"))), totss = 
> 86060500.4826787, withinss = c(21892651.1737651, 44156468.7023972), 
> tot.withinss = 66049119.8761623, betweenss = 20011380.6065164, size = 
> c(334L, 742L)), .Names = c("cluster"!
 , !
>  "centers", "totss", "withinss", "tot.withinss", "betweenss", "size"), class 
> = "kmeans")
> Please help,Suhaila 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
  
[[alternative HTML version deleted]]

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


Re: [R] Why can we combine design matrix and data-frame in R?

2012-05-13 Thread Michael
Thank you!

But the line you cited was about "response" being a matrix, which is not
our case.

And also I have checked:

Any more thoughts?

Thank you!

> names(tmp)

[1] "School" "Minority" "Sex" "SES" "MathAch" "MEANSES.x"

[7] "Size" "Sector" "PRACAD" "DISCLIM" "HIMINTY" "MEANSES.y"

[13] "Intercept" "y" "X"

> names(MathScores)

[1] "School" "Minority" "Sex" "SES" "MathAch" "MEANSES.x"

[7] "Size" "Sector" "PRACAD" "DISCLIM" "HIMINTY" "MEANSES.y"


> class(tmp)

[1] "data.frame"

> head(tmp)

School Minority Sex SES MathAch MEANSES.x Size Sector PRACAD

1 1224 No Female -1.5280 5.876 -0.428 842 Public 0.35

2 1224 No Female -0.83155757 19.708 -0.428 842 Public 0.35

3 1224 No Male -0.91452283 20.349 -0.428 842 Public 0.35

4 1224 No Male -1.3360 8.781 -0.428 842 Public 0.35

5 1224 No Male -0.35329874 17.898 -0.428 842 Public 0.35

6 1224 No Male 0.05388877 4.583 -0.428 842 Public 0.35

DISCLIM HIMINTY MEANSES.y Intercept y X.(Intercept) X.SES

1 1.597 0 -0.428 1.00 5.87600 1. -1.5280

2 1.597 0 -0.428 1.414214 27.87132 1.41421356 -0.83155757

3 1.597 0 -0.428 1.732051 35.24550 1.73205081 -0.91452283

4 1.597 0 -0.428 2.00 17.56200 2. -1.3360

5 1.597 0 -0.428 2.236068 40.02114 2.23606798 -0.35329874

6 1.597 0 -0.428 2.449490 11.22601 2.44948974 0.05388877

X.factor(Sector)Public

1 1.

2 1.41421356

3 1.73205081

4 2.

5 2.23606798

6 2.44948974


On Sat, May 12, 2012 at 6:47 AM, S Ellison  wrote:

> >I am trying to understand why this line works:
> >
> > lm1x = lm(y~X-1, tmp)
>
> Well, I would not normally define a data frame element as a matrix myself
>  (though I might well define a list element as one). But specifying a
> matrix as the terms part of an lm is documented in lm's details:
> "If response is a matrix a linear model is fitted separately by
> least-squares to each column of the matrix"
>
> So _something_ will happen.
>
> Whether the something is useful depends on the intent.
>
> > Here it seems that I was combining the design matrix and the data
> frame...
> Did you inspect tmp after adding the design matrix? Was it an odd looking
> data frame or a list?
> What seems to have been done is that the design matrix has been added to a
> list. I wouldn't normally do that if tmp is a data frame, and r would not
> do so unless the lengths all matched.  But a list should be ok. And lm
> takes a list or environment as its data argument, so a list of things will
> work even if they are different types. In other words tmp is just a ragbag
> of things, each of which lm understands.
> ***
> This email and any attachments are confidential. Any u...{{dropped:11}}

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


Re: [R] Why can we combine design matrix and data-frame in R?

2012-05-13 Thread Michael
 It's very interesting that the names of tmp has only 14 items, and the
dim(tmp) shows 14 columns too, but actually the number of columns in tmp
should be 16...

This is very strange...

> dim(X)

[1] 7185 3

> dim(MathScores)

[1] 7185 12

> dim(tmp)

[1] 7185 14


> names(tmp)

[1] "School" "Minority" "Sex" "SES" "MathAch" "MEANSES.x"

[7] "Size" "Sector" "PRACAD" "DISCLIM" "HIMINTY" "MEANSES.y"

[13] "y" "X"




On Sun, May 13, 2012 at 3:30 PM, Michael  wrote:

> Thank you!
>
> But the line you cited was about "response" being a matrix, which is not
> our case.
>
> And also I have checked:
>
> Any more thoughts?
>
> Thank you!
>
> >
> names(tmp)
>
> [1] "School" "Minority" "Sex" "SES" "MathAch" "MEANSES.x"
>
> [7] "Size" "Sector" "PRACAD" "DISCLIM" "HIMINTY" "MEANSES.y"
>
> [13] "Intercept" "y" "X"
>
> >
> names(MathScores)
>
> [1] "School" "Minority" "Sex" "SES" "MathAch" "MEANSES.x"
>
> [7] "Size" "Sector" "PRACAD" "DISCLIM" "HIMINTY" "MEANSES.y"
>
>
> >
> class(tmp)
>
> [1] "data.frame"
>
> >
> head(tmp)
>
> School Minority Sex SES MathAch MEANSES.x Size Sector PRACAD
>
> 1 1224 No Female -1.5280 5.876 -0.428 842 Public 0.35
>
> 2 1224 No Female -0.83155757 19.708 -0.428 842 Public 0.35
>
> 3 1224 No Male -0.91452283 20.349 -0.428 842 Public 0.35
>
> 4 1224 No Male -1.3360 8.781 -0.428 842 Public 0.35
>
> 5 1224 No Male -0.35329874 17.898 -0.428 842 Public 0.35
>
> 6 1224 No Male 0.05388877 4.583 -0.428 842 Public 0.35
>
> DISCLIM HIMINTY MEANSES.y Intercept y X.(Intercept) X.SES
>
> 1 1.597 0 -0.428 1.00 5.87600 1. -1.5280
>
> 2 1.597 0 -0.428 1.414214 27.87132 1.41421356 -0.83155757
>
> 3 1.597 0 -0.428 1.732051 35.24550 1.73205081 -0.91452283
>
> 4 1.597 0 -0.428 2.00 17.56200 2. -1.3360
>
> 5 1.597 0 -0.428 2.236068 40.02114 2.23606798 -0.35329874
>
> 6 1.597 0 -0.428 2.449490 11.22601 2.44948974 0.05388877
>
> X.factor(Sector)Public
>
> 1 1.
>
> 2 1.41421356
>
> 3 1.73205081
>
> 4 2.
>
> 5 2.23606798
>
> 6 2.44948974
>
>
>  On Sat, May 12, 2012 at 6:47 AM, S Ellison wrote:
>
>> >I am trying to understand why this line works:
>> >
>> > lm1x = lm(y~X-1, tmp)
>>
>> Well, I would not normally define a data frame element as a matrix myself
>>  (though I might well define a list element as one). But specifying a
>> matrix as the terms part of an lm is documented in lm's details:
>> "If response is a matrix a linear model is fitted separately by
>> least-squares to each column of the matrix"
>>
>> So _something_ will happen.
>>
>> Whether the something is useful depends on the intent.
>>
>> > Here it seems that I was combining the design matrix and the data
>> frame...
>> Did you inspect tmp after adding the design matrix? Was it an odd looking
>> data frame or a list?
>> What seems to have been done is that the design matrix has been added to
>> a list. I wouldn't normally do that if tmp is a data frame, and r would not
>> do so unless the lengths all matched.  But a list should be ok. And lm
>> takes a list or environment as its data argument, so a list of things will
>> work even if they are different types. In other words tmp is just a ragbag
>> of things, each of which lm understands.
>> ***
>> This email and any attachments are confidential. Any use, copying or
>> disclosure other than by the intended recipient is unauthorised. If
>> you have received this message in error, please notify the sender
>> immediately via +44(0)20 8943 7000 or notify postmas...@lgcgroup.com
>> and delete this message and any copies from your computer and network.
>> LGC Limited. Registered in England 2991879.
>> Registered office: Queens Road, Teddington, Middlesex, TW11 0LY, 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] Discrete choice model maximum likelihood estimation

2012-05-13 Thread infinitehorizon
Hello,

I am new to R and I am trying to estimate a discrete model with three
choices. I am stuck at a point and cannot find a solution.

I have probability functions for occurrence of these choices, and then I
build the likelihood functions associated to these choices and finally I
build the general log-likelihood function.

There are four parameters in the model, three of them are associated to
three discrete choices I mentioned, and one of them is for a binary variable
in the data (t).  There are also latent variables but I didn't put them in
this question because if I figure out how to do this, I will be able to add
them as well.

I am not familiar with the syntax I have to write in the likelihood
functions, so I really doubt that they are true.  Below I simplify the
problem and provide the code I've written:

# Probabilities for discrete choices for a=3, a=2 and a=1 respectively
P3  <- function(b3,b,t) {
P   <- exp(b3+b*(t==1))/(1-exp(b3+b*(t==1)))
return(P)   
}
P2  <- function(b2,b,t) {
P   <- exp(b2+b*(t==1))/(1-exp(b2+b*(t==1)))
return(P)   
}
P1  <- function(b1,b,t) {
P   <- exp(b1+b*(t==1))/(1-exp(b1+b*(t==1)))
return(P)   
}

# Likelihood functions for discrete choices for a=3, a=2 and a=1
respectively

L3  <- function(b1,b2,b3,b,t) {
P11 <- P1(b1,b,t)
P22 <- P2(b2,b,t)
P33 <- P3(b3,b,t)

L3l <- (P11=1)*(P22=1)*(P33=1)
return(L3l)
}

L2  <- function(b1,b2,b3,b,t) {
P11 <- P1(b1,b,t)
P22 <- P2(b2,b,t)
P33 <- P3(b3,b,t)

L2l <- (P11=1)*(P22=1)*(P33=0)
return(L2l)
}

L1  <- function(b1,b2,b,t) {
P11 <- P1(b1,b,t)
P22 <- P2(b2,b,t)

L1l <- (P11=1)*(P22=0)
return(L1l)
}

# Log-likelihood function

llfn<- function(par,a,t) {

b1  <- par[1]
b2  <- par[2]
b3  <- par[3]
b   <- par[4]

lL1 <- log(L1(b1,b2,b,t))
lL2 <- log(L2(b1,b2,b3,b,t))
lL3 <- log(L3(b1,b2,b3,b,t))

llfn<- (a==1)*lL1+(a==2)*lL2+(a==3)*lL3
}
est <- optim(par,llfn, method = c("CG"),control=list(trace=2,maxit=2000),
hessian=TRUE)

And when I run this code I get "cannot coerce type 'closure' to vector of
type 'double'" error.
I will really appreciate your help. Thanks,



--
View this message in context: 
http://r.789695.n4.nabble.com/Discrete-choice-model-maximum-likelihood-estimation-tp4629877.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] CRAN (and crantastic) updates this week

2012-05-13 Thread Crantastic
CRAN (and crantastic) updates this week

New packages


* AGD (0.30)
  Maintainer: Stef van Buuren
  Author(s): Stef van Buuren 
  License: GPL-2 | GPL-3
  http://crantastic.org/packages/AGD

  Tools for NIHES course EP18 'Analysis of Growth Data', May 22-23 2012,
  Rotterdam.

* bigml (0.1-1)
  Maintainer: Justin Donaldson
  Author(s): Justin Donaldson [aut, cre]
  License: LGPL-3
  http://crantastic.org/packages/bigml

  The bigml package contains bindings for the BigML API. The package
  includes methods that provide straightforward access to basic API
  functionality, as well as methods that accommodate idiomatic R
  datatypes and concepts.

* BVS (4.12.0)
  Maintainer: Melanie Quintana
  Author(s): Melanie Quintana 
  License: Unlimited
  http://crantastic.org/packages/BVS

  The functions in this package focus on analyzing case-control
  association studies involving a group of genetic variants.  In
  particular, we are interested in modeling the outcome variable as a
  function of a multivariate genetic profile using Bayesian model
  uncertainty and variable selection techniques.  The package
  incorporates functions to analyze data sets involving common
  variants as well as extensions to model rare variants via the
  Bayesian Risk Index (BRI) as well as haplotypes.  Finally, the
  package also allows the incorporation of external biological
  information to inform the marginal inclusion probabilities via the
  iBMU.

* gsbDesign (0.95)
  Maintainer: Unknown
  Author(s): Florian Gerber, Thomas Gsponer
  License: GPL-3
  http://crantastic.org/packages/gsbDesign

  Group Sequential Operating Characteristics for Clinical, Bayesian
  two-arm Trials with known Sigma and Normal Endpoints.

* marqLevAlg (1.0)
  Maintainer: Melanie Prague
  Author(s): D. Commenges , M. Prague
  and Amadou Diakite
  License: GPL (>= 2.0)
  http://crantastic.org/packages/marqLevAlg

  This algorithm provides a numerical solution to the problem of
  minimizing a function. This is more efficient than the
  Gauss-Newton-like algorithm when starting from points vey far from
  the final minimum. A new convergence test is implemented (RDM) in
  addition to the usual stopping criterion : stopping rule is when the
  gradients are small enough in the parameters metric (GH-1G).

* p2distance (1.0.1)
  Maintainer: A.J. Perez-Luque
  Author(s): A.J. Perez-Luque; R. Moreno; R. Perez-Perez and F.J. Bonet
  License: GPL
  http://crantastic.org/packages/p2distance

  The welfare's synthetic indicator provides an ideal tool for measuring
  multi-dimensional concepts such as welfare, development, living
  standards, etc. It enables information from the various indicators
  to be aggregated into a single synthetic measure.

* pcaPA (1.0)
  Maintainer: Carlos A. Arias
  Author(s): Carlos A. Arias  and Víctor H. Cervantes
 .
  License: GPL (>= 2)
  http://crantastic.org/packages/pcaPA

  A set of functions to perform parallel analysis for principal
  components analysis intended mainly for large data sets. It performs
  a parallel analysis of continuous, ordered (including
  dichotomous/binary as a special case) or mixed type of data
  associated with a principal components analysis. Polychoric
  correlations among ordered variables, Pearson correlations among
  continuous variables and polyserial correlation between mixed type
  variables (one ordered and one continuous) are used. Whenever the
  use of polyserial or polychoric correlations yields a non positive
  definite correlation matrix, the resulting matrix is transformed
  into the nearest positive definite matrix. For details on the GSL
  installation in windows please visit the following link:
  http://wiki.rglab.org/index.php?title=Public:Installation_of_GSL_from
  Windows_binary

* phia (0.0-2)
  Maintainer: Helios De Rosario-Martinez
  Author(s): Helios De Rosario-Martinez
  License: GPL (>= 3)
  http://crantastic.org/packages/phia

  Analysis of terms in linear and generalized linear models, on the
  basis of multiple comparisons of factor contrasts. Specially suited
  for the analysis of interaction terms.

* plmDE (1.0)
  Maintainer: Jonas
  Author(s): Jonas Mueller
  License: GPL (>= 2)
  http://crantastic.org/packages/plmDE

  A set of tools for identifying genes whose differential expression is
  associated with measurements of other covariates on a continuous
  scale.  These methods rely on generalized additive partially linear
  models which can be fitted efficiently using a B-spline basis
  approximation.  Still under development: methods for interfacing
  with objects extending the eSet class and a function to pass linear
  models in edgeR and DEseq format.

* polywog (0.1-0)
  Maintainer: Brenton Kenkel
  Author(s): Brenton Kenkel and Curtis S. Signorino
  License: GPL (>= 2)
  http://crantastic.org/packages/polywog

  Routines for flexible functional form estimation via basis regression,
  with model selection via the adaptive LASSO or 

Re: [R] Discrete choice model maximum likelihood estimation

2012-05-13 Thread Rui Barradas
Hello,

There are several issues with your code.

1. The error message. Don't use 'par' as a variable name, it's already an R
function, tyo get or set graphics parameters.
Call it something else, say, 'param'.
This is what causes the error. You must pass initial values to optim, but
the variable you're passing doesn't exist, you haven't created it so R finds
an object with that name, the graphics parameters function.
Avoid the confusion.
And create 'param' with as many values as expected by llfn before the call.

2. 't' is also a parameter. Take it out of llfn formals and put it in
'param'. Then, inside llfn's body,

t <- param[5]

3. It still won't work. llfn will not be passed a value for 'a', for the
same reason it can't find 't'.

4. Then, look at L3 and the others. The line just before return.

L3l <- (P11=1)*(P22=1)*(P33=1)

After computing P11, etc, you're discarding those values and assigning 1 to
each of them.
Your likelihood functions just became constants...
And if this is a typo, if you meant P11 == 1, etc, it's even worse. You
can't expect that ratios of exponentials to be equal to that one real value.

Points 1-3 are workable but this last one means you have to revise your
likelihood.
Good luck.

Hope this helps,

Rui Barradas

infinitehorizon wrote
> 
> Hello,
> 
> I am new to R and I am trying to estimate a discrete model with three
> choices. I am stuck at a point and cannot find a solution.
> 
> I have probability functions for occurrence of these choices, and then I
> build the likelihood functions associated to these choices and finally I
> build the general log-likelihood function.
> 
> There are four parameters in the model, three of them are associated to
> three discrete choices I mentioned, and one of them is for a binary
> variable in the data (t).  There are also latent variables but I didn't
> put them in this question because if I figure out how to do this, I will
> be able to add them as well.
> 
> I am not familiar with the syntax I have to write in the likelihood
> functions, so I really doubt that they are true.  Below I simplify the
> problem and provide the code I've written:
> 
> # Probabilities for discrete choices for a=3, a=2 and a=1 respectively
> P3<- function(b3,b,t) {
> P <- exp(b3+b*(t==1))/(1-exp(b3+b*(t==1)))
> return(P) 
> }
> P2<- function(b2,b,t) {
> P <- exp(b2+b*(t==1))/(1-exp(b2+b*(t==1)))
> return(P) 
> }
> P1<- function(b1,b,t) {
> P <- exp(b1+b*(t==1))/(1-exp(b1+b*(t==1)))
> return(P) 
> }
> 
> # Likelihood functions for discrete choices for a=3, a=2 and a=1
> respectively
> 
> L3<- function(b1,b2,b3,b,t) {
> P11   <- P1(b1,b,t)
> P22   <- P2(b2,b,t)
> P33   <- P3(b3,b,t)
> 
> L3l   <- (P11=1)*(P22=1)*(P33=1)
> return(L3l)
> }
> 
> L2<- function(b1,b2,b3,b,t) {
> P11   <- P1(b1,b,t)
> P22   <- P2(b2,b,t)
> P33   <- P3(b3,b,t)
> 
> L2l   <- (P11=1)*(P22=1)*(P33=0)
> return(L2l)
> }
> 
> L1<- function(b1,b2,b,t) {
> P11   <- P1(b1,b,t)
> P22   <- P2(b2,b,t)
> 
> L1l   <- (P11=1)*(P22=0)
> return(L1l)
> }
> 
> # Log-likelihood function
> 
> llfn  <- function(par,a,t) {
> 
> b1<- par[1]
> b2<- par[2]
> b3<- par[3]
> b <- par[4]
> 
> lL1   <- log(L1(b1,b2,b,t))
> lL2   <- log(L2(b1,b2,b3,b,t))
> lL3   <- log(L3(b1,b2,b3,b,t))
> 
> llfn  <- (a==1)*lL1+(a==2)*lL2+(a==3)*lL3
> }
> est   <- optim(par,llfn, method = c("CG"),control=list(trace=2,maxit=2000),
> hessian=TRUE)
> 
> And when I run this code I get "cannot coerce type 'closure' to vector of
> type 'double'" error.
> I will really appreciate your help. Thanks,
> 


--
View this message in context: 
http://r.789695.n4.nabble.com/Discrete-choice-model-maximum-likelihood-estimation-tp4629877p4629880.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] Why can we combine design matrix and data-frame in R?

2012-05-13 Thread S Ellison
> But the line you cited was about "response" being a matrix, which is not our 
> case.
Yes, you're right; I picked the wrong thing to cite.
The only documentation I found about lm accepting a matrix in the predictors is 
a one-line statement in "Introduction to R" which says "term_i
is either

a vector or matrix expression, or 1,
a factor, or
a formula expression consisting of factors, vectors or matrices 
connected by formula operators. "

Not the most informative documentation. But Peter Dalgaard is a most 
authoritative source!

>And also I have checked:
>
>Any more thoughts?

Data frames are odd things; a column need not contain only a vector if the 
number of rows is OK. I am half surprised that including a matrix in one works. 
But the gods of R are powerful and their magic is strong. Here, names(tmp) is 
showing that the data frame has one element called X (in effect, the whole 
matrix is regarded as one element of the data frame), but on display the magic 
has expanded X to show all the columns of X.

This is the main reason I generally keep to simple things in data frames; 
complicated things make it less easy to predict behaviour.



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

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


Re: [R] Why can we combine design matrix and data-frame in R?

2012-05-13 Thread Luna
Thanks!

Do you think if the correctness of the such results could be generalized to
other future cases?




On Sun, May 13, 2012 at 7:10 PM, S Ellison  wrote:

> > But the line you cited was about "response" being a matrix, which is not
> our case.
> Yes, you're right; I picked the wrong thing to cite.
> The only documentation I found about lm accepting a matrix in the
> predictors is a one-line statement in "Introduction to R" which says "term_i
>is either
>
>a vector or matrix expression, or 1,
>a factor, or
>a formula expression consisting of factors, vectors or matrices
> connected by formula operators. "
>
> Not the most informative documentation. But Peter Dalgaard is a most
> authoritative source!
>
> >And also I have checked:
> >
> >Any more thoughts?
>
> Data frames are odd things; a column need not contain only a vector if the
> number of rows is OK. I am half surprised that including a matrix in one
> works. But the gods of R are powerful and their magic is strong. Here,
> names(tmp) is showing that the data frame has one element called X (in
> effect, the whole matrix is regarded as one element of the data frame), but
> on display the magic has expanded X to show all the columns of X.
>
> This is the main reason I generally keep to simple things in data frames;
> complicated things make it less easy to predict behaviour.
>
>
>
> ***
> This email and any attachments are confidential. Any u...{{dropped:13}}

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


[R] Post stratification weights in survey package in R

2012-05-13 Thread Ruijie
Hi all,

I have data collected from a survey administered on a subset of the
population. I also have the population proportions of variables such as
gender, race and housing type. I would like to combine the weights from
each separate cross tab (of gender, race and housing type) such that the
weighted proportions of my survey data matches that of the population.

I have tried the following:

library(survey)


gender.population <-
read.table("http://dl.dropbox.com/u/822467/Gender.csv";, header = TRUE,
sep = ",")

housing.population <-
read.table("http://dl.dropbox.com/u/822467/Housing.csv";, header =
TRUE, sep = ",")

race.population <-
read.table("http://dl.dropbox.com/u/822467/Race.csv";, header = TRUE,
sep = ",")

survey.sample <-
read.table("http://dl.dropbox.com/u/822467/survey.sample.csv";, header
= TRUE, sep = ",")

survey.object.sample <- svydesign(id = ~1, data = survey.sample)

survey.object.sample.weighted <- rake(survey.object.sample,
list(~gender, ~housing, ~race), list(gender.population,
housing.population, race.population))

str(survey.object.sample.weighted$postStrata)

I see from survey.object.sample.weighted$postStrata that weights have been
assigned separately for each of the variable. My question is: Is it
possible to get 1 weight for each subject instead of 3 weights as shown in
the package?
Regards,
Ruijie (RJ)

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