[R] Variance component estimation in glmmPQL
Hi all, I've been attempting to fit a logistic glmm using glmmPQL in order to estimate variance components for a score test, where the model is of the form logit(mu) = X*a+ Z1*b1 + Z2*b2. Z1 and Z2 are actually reduced rank square root matrices of the assumed covariance structure (up to a constant) of random effects c1 and c2, respectively, such that b1 ~ N(0,sig.1^2*I) and c1 ~ N(0,sig.1^2*K1) , where K1 = Z1*t(Z1), and c1 = Z1*b1. The model form I've been using is just the following: m<-glmmPQL(y~1,random=list(f=pdBlocked(list(pdIdnot(~Z.1-1),pdIdnot(~Z.2-1 ,family=binomial(link="logit")) I've been extracting the variance components using VarCorr(), but I've noticed that the reported variances associated with my random effects are not even close to the values I get if I evaluate their variances empirically (eg var(random.effects(m.12)). I know that's not how they're actually estimated, but there may be a whole order of magnitude difference in the values. Below is an example under R 2.14 on a linux machine: library(MASS) library(mgcv) library(boot) set.seed(1234) G.1<-matrix(rnorm(5000,0,0.25),nrow=100) G.2<-matrix(rnorm(5000,0,0.25),nrow=100) K.1<-G.1%*%t(G.1) K.2<-G.2%*%t(G.2) Z.1<-mroot(K.1,method="svd") Z.2<-mroot(K.2,method="svd") b.1<-matrix(rnorm(ncol(Z.1),0,0.25),ncol=1) b.2<-matrix(rnorm(ncol(Z.2),0,0.50),ncol=1) p<-inv.logit(Z.1%*%b.1+Z.2%*%b.2) y<-rbinom(100,1,p) f<-rep(1,100) m.fit<-glmmPQL(y~1,random=list(f=pdBlocked(list(pdIdent(~Z.1-1),pdIdent(~Z.2-1, family=binomial(link="logit")) VarCorr(m.fit) var(as.numeric(random.effects(m.fit))[1:ncol(Z.1)]) var(as.numeric(random.effects(m.fit))[-c(1:ncol(Z.1))]) >From the above, VarCorr() results in variance component estimates of sig.1^2 = 0.444 and sig.2^2 = 0.2778, whereas the empirical estimates are sig.1^2 = 0.2060 and sig.1^2 = 0.097. I know variance component estimation in general is a little shaky, but my simulations suggest that VarCorr() is extracting values that are way too large on a consistent basis. I'm largely assuming I'm misinterpreting something here, just can't figure out what. -Nick -- View this message in context: http://r.789695.n4.nabble.com/Variance-component-estimation-in-glmmPQL-tp4650964.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] Difficulty with 'loess' function
Flip date and q in your formula, you've got them backwards from what you've said you're trying to model. armstrwa wrote: > > > Yeah, I did look at the help(loess) page, but I wasn't really sure what to > do with that. I was inputting it as: > >> test<-loess(date ~ q,data.frame(date,q),span=0.5) > > So I guess I'm inputting the formula wrong then? But from looking at the > help(formula) page, I still wasn't sure what it's looking for. > > Thanks. > > Billy > -- View this message in context: http://r.789695.n4.nabble.com/Difficulty-with-loess-function-tp3387537p3388108.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] Matching two vectors
Just use vect_2_id as a subsetting index for vect_1, ie:
vect_2<-vect_1[vect_2_id]
Vincy Pyne wrote:
>
> Dear R helpers
>
> Suppose I have a vector as
>
> vect_1 = c("AAA", "AA", "A", "BBB", "BB", "B", "CCC")
>
> vect_1_id = c(1:length(vect_1))
>
> Through some process I obtain
>
> vect_2_id = c(2, 3, 7), then I need a new vector say vect_2 which will
> give me
>
> vect2 = ("AA", "A", "CCC") i.e. I need the subset of vect_1 as per
> vect_2_id.
>
> Thanking in advance
>
> Regards
>
> Vincy
>
>
>
>
>
> [[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.
>
--
View this message in context:
http://r.789695.n4.nabble.com/Matching-two-vectors-tp3357591p3357603.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] create data set from selection of rows
Try using which() Something like data.frame(dataset[which(dataset[,3]=="text3"),]) e-letter wrote: > > Readers, > > For a data set: > > text1,23,text2,45 > text1,23,text3,78 > text1,23,text3,56 > text1,23,text2,45 > > The following command was entered: > > datasubset<-data.frame(dataset[,3]=="text3") > > The result of > > datasubset > > is > > TRUE > TRUE > > The required result is > > text1,23,text3,78 > text1,23,text3,56 > > What is the correct command to use please? > > Thanks in advance. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > -- View this message in context: http://r.789695.n4.nabble.com/create-data-set-from-selection-of-rows-tp3356866p3356879.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] JAGS/BUGS on gene expression data
64 Bit R w/JAGS seems to be stalling out as well, I ran a test run of 100 iterations and it's been hanging for 8 hours so that doesn't seem to be the solution. I'll take a look at PYMC. That CppBUGS package looks pretty interesting, I'll keep my eye on it. My C programming book arrives today from Amazon...so if all else fails... -- View this message in context: http://r.789695.n4.nabble.com/JAGS-BUGS-on-gene-expression-data-tp3354065p3356723.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] Finding coordinates for maximum of a function
That actually won't work. max(y) will give a value, not a coordinate, so x[max(y)] is definitely not what you want. -- View this message in context: http://r.789695.n4.nabble.com/Finding-coordinates-for-maximum-of-a-function-tp3355369p3356483.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] Finding coordinates for maximum of a function
This is where which.max() comes in handy
n<-length(x)
x.c<-rep(0,n)
for(i in 1:n){
x.c[i]<-which.max[y1]
}
x.c is then a vector of x coordinates for the maximum for columns
y1,y2,...,yn
--
View this message in context:
http://r.789695.n4.nabble.com/Finding-coordinates-for-maximum-of-a-function-tp3355369p3355395.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] code for "permutative" operation
You mean like cumsum()? > a<-1:10 > b<-cumsum(a) > b [1] 1 3 6 10 15 21 28 36 45 55 -Nick Larson -- View this message in context: http://r.789695.n4.nabble.com/code-for-permutative-operation-tp3354839p3354849.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] creating character vector
The reason this doesn't work is because R thinks that in the command as.character(c(first,second)) that first and second are variables that exist within the working environment. Since they don't (I assume), R doesn't know what to do with the command. Using the quotes indicates to R that you're specifying a character string type, instead of calling a variable. As other have said before, if you don't want to type in quotes, you could just read in a text file in the methods they describe. -- View this message in context: http://r.789695.n4.nabble.com/creating-character-vector-tp3353447p3354199.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] JAGS/BUGS on gene expression data
Has anybody had issues running MCMC (either BUGS or JAGS) on data sets of this magnitude (ie 30k x 20-30). I've been trying to run a hierarchical random effects model on expression data but R completely stalls out on jobs run on 32bit R on our server (doesn't respond...then eventually crashes out with an exception fault). Would using 64bit R help with JAGS? (As far as I know there's only a 32bit build of BRugs so BUGS is out) Currently the server doesn't have 64 bit JAGS installed but I've put in a request. It seems that compiling the model is what's crashing R at this point. My last resort is to hardcode everything in C, but I'd like to avoid that. -- View this message in context: http://r.789695.n4.nabble.com/JAGS-BUGS-on-gene-expression-data-tp3354065p3354065.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.
