[R] mb.long {timecourse}

2012-03-29 Thread statfan
How do you identify the genes which are differentially expressed using the
mb.long function?  More specifically, in the fruitfly example (see below),
we begin with an expression matrix containing 2000 genes.  How do I obtain
not only the proportion of differentially expressed genes, but also the
subset of genes which are differentially expressed?

data(fruitfly)
colnames(fruitfly) ## check if arrays are arranged in the default order
gnames - rownames(fruitfly)
assay - rep(c(A, B, C), each = 12)
time.grp - rep(c(1:12), 3)
size - rep(3, nrow(fruitfly))

out1 - mb.long(fruitfly, times=12, reps=size, rep.grp = assay, time.grp =
time.grp)
summary(out1)
plotProfile(out1, type=b, gnames=gnames, legloc=c(2,15),
pch=c(A,B,C), xlab=Hour)

 summary(out1)
Length Class  Mode   
M   72000  -none- numeric
prop1  -none- numeric
nu  1  -none- numeric
Lambda1   121  -none- numeric
percent 2  -none- numeric
size 2000  -none- numeric
con.group  36  -none- numeric
rep.group  36  -none- numeric
time.group 36  -none- numeric
HotellingT2  2000  -none- numeric
pos.HotellingT2  2000  -none- numeric

Thanks to anyone who can shed light on this matter.

--
View this message in context: 
http://r.789695.n4.nabble.com/mb-long-timecourse-tp4515669p4515669.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] limma design matrix

2012-03-25 Thread statfan
I am trying to construct my design matrix needed in the {limma} function
lmfit but am having trouble with the formula I am to specify in the
function model.matrix.  Namely when to I use ~0 + factors (ex 1) vs ~-1 +
factors (ex 2).  Any clarification on this would be greatly appreciated. 
Thanks in advanced.  

Ex 1:
  f - factor(targets$Target, levels=c(RNA1,RNA2,RNA3))
design - model.matrix(~0+f)
[page 409, of Bioinformatics and computational biology solutions using R and
Bioconductor]

and 

Ex 2:
library(ALL)
pdat - pData(ALL)
design - model.matrix(~-1 + factor(pdat$type) 
[page 237, of Bioinformatics and computational biology solutions using R and
Bioconductor]



--
View this message in context: 
http://r.789695.n4.nabble.com/limma-design-matrix-tp4503651p4503651.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] hypergeometric function in ‘ mvtnorm’

2012-03-19 Thread statfan
Thanks for your advice.  I actually meant to ask about the pmvt for the
distribution function.  Viewing the source code pmvt uses the function
mvt which uses the function probval which sources the fortran code:

Fortran(mvtdst, N = as.integer(n), NU = as.integer(df), 
LOWER = as.double(lower), UPPER = as.double(upper), INFIN =
as.integer(infin), 
CORREL = as.double(corrF), DELTA = as.double(delta), 
MAXPTS = as.integer(x$maxpts), ABSEPS = as.double(x$abseps), 
RELEPS = as.double(x$releps), error = as.double(error), 
value = as.double(value), inform = as.integer(inform), 
PACKAGE = mvtnorm)

I wish to look at how this mvtdst calculates the hypergeometric function
(2_F_1).  Anyway that I can see that?
Thanks 

--
View this message in context: 
http://r.789695.n4.nabble.com/hypergeometric-function-in-mvtnorm-tp4483730p4485277.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] hypergeometric function in ‘ mvtnorm’

2012-03-18 Thread statfan
Is there any way to know how the dmvt function computes the hypergeometric
function needed in the calculation for the density of multivariate t
distribution?  

--
View this message in context: 
http://r.789695.n4.nabble.com/hypergeometric-function-in-mvtnorm-tp4483730p4483730.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] try to silence errors

2011-12-22 Thread statfan
I am trying to use the dmt function in the package {mnormt}.  Throughout my
algorithm, the covariance matrix is sometime calculated to be singular. 
When attempting to calculate the dmt function with a covariance that is not
positive definite, I would like it to return Inf or NaN instead of an error
message.

I have been using the try function, however it is not yeilding the desired
result. (I did not include the values since they are 17 dimensions, all you
need to know is that sigij is singular)

try(out - dmt(y[,j],new$mu[,,4],sigij,ceiling(new$nu[4])),silent=TRUE)
Error in chol.default(x, pivot = FALSE) : 
  the leading minor of order 16 is not positive definite

I have similarly tried the trycatch function (see below for code) however I
am using this within a function made my be  (which i call ll) so I get
errors like
 ll - function(new,y){ 
...
out - tryCatch(dmt(y[,j],new$mu[,,3],sigij,
ceiling(new$nu[3])),error = function(e) -Inf)

}

ll(new,y)
Error in chol.default(x, pivot = FALSE) : 
  the leading minor of order 16 is not positive definite
[1] -Inf

If anyone could shed some light on this problem, it would be greatly
appreciated.


--
View this message in context: 
http://r.789695.n4.nabble.com/try-to-silence-errors-tp4225754p4225754.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] breaking up n object into random groups

2011-12-02 Thread statfan
say n = 100
I want to partition this into 4 random groups wheren n1 + n2 + n3 + n4 = n
and ni is the number of elements in group i.

Thank you for you help

--
View this message in context: 
http://r.789695.n4.nabble.com/breaking-up-n-object-into-random-groups-tp4147476p4147476.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] tryCatch with integration

2011-10-31 Thread statfan
I am trying to compute the approximate numerical integration of the following
expression using the integrate function:


   (integrate(function(x) {log(1+x^2)*(1+x^2)^(-20.04543)},low,Inf)$val)
Error in integrate(function(x) { : the integral is probably divergent


which gives me an error.  If this error happens I want to return the
following instead

 (area(function(x) {-2*log(cos(x))*(cos(x)^(2*e-2))}, atan(low), pi/2))
[1] 0.01045636

In order to accomplish this I tried the following:

 this= tryCatch((integrate(function(x)
 {log(1+x^2)*(1+x^2)^(-e)},low,Inf)$val),error=function(e) e, finally =
 (area(function(x) {-2*log(cos(x))*(cos(x)^(2*e-2))}, atan(low), pi/2)))
 this
simpleError in integrate(function(x) {log(1 + x^2) * (1 + x^2)^(-e)},
low, Inf): the integral is probably divergent

Perhaps I am using tryCatch incorrectly? Any help would be greatly
appreciated.

--
View this message in context: 
http://r.789695.n4.nabble.com/tryCatch-with-integration-tp3954528p3954528.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] symmetric matrix multiplication

2011-10-23 Thread statfan
I have a symmetric matrix B (17x17), and a (17x17) square matrix A.  If do
the following matrix multiplication I SHOULD get a symmetric matrix, however
i don't.  The computation required is:

C = t(A)%*%B%*%A

here are some checks for symmetry
 (max(abs(B - t(B
[1] 0
 C = t(A)%*%B%*%A
 (max(abs(C - t(C
[1] 3.552714e-15

Any help on the matter would be very much appreciated.


--
View this message in context: 
http://r.789695.n4.nabble.com/symmetric-matrix-multiplication-tp3929739p3929739.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] symmetric matrix multiplication

2011-10-23 Thread statfan
Thank you Dan and Ted for these helpful comments. I will implement this
simple force symmetry code you suggested and make sure I familiarize with
this floating-point calculation problem so I can recognize such issues in
the future.

--
View this message in context: 
http://r.789695.n4.nabble.com/symmetric-matrix-multiplication-tp3929739p3930832.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] sampling from the multivariate truncated normal

2011-07-05 Thread statfan

Well, for 0.828324  x[2]  Inf the probablility is roughly 0 hence not 
easy to draw random numbers out there  

Uwe Ligges

How is this probability roughly 0?

--
View this message in context: 
http://r.789695.n4.nabble.com/sampling-from-the-multivariate-truncated-normal-tp3626438p3647039.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] sampling from the multivariate truncated normal

2011-06-26 Thread statfan
I am trying generate a sample for a truncated multivariate normal
distribution via the rtmvnorm function in the  {tmvtnorm} package.

Why does the following produce NaNs?


rtmvnorm(1, mean = rep(0, 2), matrix(c(0.06906084, -0.07463565, -0.07463565,
0.08078086),2),c(-0.4316738,  0.8283240),  c(Inf,Inf), algorithm=gibbsR,
burn.in.samples=100)

Thanks

--
View this message in context: 
http://r.789695.n4.nabble.com/sampling-from-the-multivariate-truncated-normal-tp3626438p3626438.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] rtmvt

2011-05-30 Thread statfan
I want to use the rtmvt from the  {tmvtnorm} package using the gibbs
algorithm but how to i specify the nested function rtmvnorm to use gibbs as
well?

Right now I am using the code:

for (i in 1:g){
for (j in 1:n){
sgamma[,,i,j] = rtmvt(n=50, mean=mu[i,j], sigma[i,j],
df=nu[i], lower=rep(0,2),algorithm=gibbs)
}
}

heres an example of one iteration:
 mu[1,1]
-0.09734357  0.51578628

 sigma[1,1]
  [,1]  [,2]
[1,] 0.4250681 0.0253649
[2,] 0.0253649 0.4250681


when I run this i get 50 errors saying:
Warning messages:
1: In rtmvnorm.rejection(n, mean, sigma, lower, upper, ...) :
  Acceptance rate is very low and rejection sampling becomes inefficient.
Consider using Gibbs sampling.

I have figured out that this is coming from the internal function rtmvnorm
who's default is the rejection algorithm.  Is there any way that I can
specify that to be Gibbs as well? 


Thanks


--
View this message in context: 
http://r.789695.n4.nabble.com/rtmvt-tp3562217p3562217.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] rtmvt

2011-04-10 Thread statfan
I have been using the rtmvt function in the {tmvtnorm} package i'm getting
the warning:

Acceptance rate is very low and rejection sampling becomes inefficient.
Consider using Gibbs sampling.

but i AM specifying the gibbs algorithm!!:

rtmvt(M, mean=q[,,i,j], sigma=((u[i,j] + nu[i])/(p+nu[i]))*delta[,,i],
df=ceiling(nu[i]+p), lower=c(0,0), algorithm=gibbs)

Any ideas why I am getting this warning and how I can fix it?  

--
View this message in context: 
http://r.789695.n4.nabble.com/rtmvt-tp3440751p3440751.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] multivariate t distribution

2011-04-06 Thread statfan
I have been working the the pmt function in the {mnormt} package and which
requires 

S   a positive definite matrix representing the scale matrix of the
distribution, such that S*df/(df-2) is the variance-covariance matrix when
df2; a vector of length 1 is also allowed (in this case, d=1 is set)

is there a way that I can specify the scale covariance matrix instead?  Or
alternatively, how do I convert the scale covariance matrix into this
positive definite S matrix.  Thanks in advanced.

--
View this message in context: 
http://r.789695.n4.nabble.com/multivariate-t-distribution-tp3432257p3432257.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] truncated distributions

2011-04-02 Thread statfan
I am sampling from the truncated multivariate student t distribution rtmvt
in the package {tmvtnorm}. My question is about the mean vector.  Is it
possible to define a mean vector outside of the truncated region? Thank you
in advance for any help.

--
View this message in context: 
http://r.789695.n4.nabble.com/truncated-distributions-tp3422245p3422245.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] truncated distributions

2011-04-02 Thread statfan
The definition of the mean vector is essentially what my question boils
down to.  In the functions details, the author states

We sample x ~ T(mean, Sigma, df) subject to the rectangular truncation
lower = x = upper. Currently, two random number generation methods are
implemented: rejection sampling and the Gibbs Sampler.

So if the mean vector in the rtmvt function is the mean of the parent
distribution's mean (as I hope it is), then it would be acceptable to define
a mean vector outside of the truncated range.  Clarification of this point
would be greatly appreciated.  

--
View this message in context: 
http://r.789695.n4.nabble.com/truncated-distributions-tp3422245p3422434.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] pmt

2011-03-27 Thread statfan
I am working with the pmt function in the {mnormt} package, and i am getting
negative values returned.  the following is an example of one of my outputs:

pmt(x = c(3.024960, -1.010898), mean = c(21.18844, 21.18844), S =
matrix(c(.319,.139,.139,0.319), 2, 2),df = 42)
# -6.585641e-18

Any help on why i'm getting negative numbers would be very much appreciated.

THanks!

--
View this message in context: 
http://r.789695.n4.nabble.com/pmt-tp3409543p3409543.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.