[R] stats/debugging question hotelling t-sq
Hi I spent hours looking over my formula. Somehow I cant find the reason why it gives me different answer. help appreciated. x = as.matrix(read.table("http://www.niehs.nih.gov/research/atniehs/core/microarrays/docs/heinloth.txt",1)) x = t(x)#now rows are subjects, cols are genes x = x[order(rownames(x)),]#order by treatment group oxygen, ultra-violet, gamma radiation y = x[26:45,1:10] x = x[2:25,1:10] p = ncol(x); p nx = nrow(x); nx ny = nrow(y); ny n = nx+ny; n # (t(x)-colMeans(x)) %*% t(t(x)-colMeans(x)) T2 = nx*ny/n * t(colMeans(x)-colMeans(y)) %*% solve( ( (nx-1)*cov(x)+(ny-1)*cov(y) )/( n-2 ) ) %*% (colMeans(x)-colMeans(y)); T2 library(ICSNP) HotellingsT2(y,x) http://en.wikipedia.org/wiki/Hotelling's_T-square_distribution http://finzi.psych.upenn.edu/R/Rhelp02a/archive/60962.html __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] lme cant get parameter estimated correctly
I am caught in a mental trap. Why isn't the between groups variance estimated (0.0038) to be around the value with which I generated the data (0.0002)? Thanks Toby set.seed(76589437887) fph = 0.4 Sigh = sqrt(0.0002) Sigi = sqrt(0.04) ci = 1 fpi = matrix(,7200,3) for (i in 1:90) { fph = rnorm(1, fph, Sigh) for (k in 1:80) { fpi[ci,1:3] = matrix(c(i, k, rnorm(1, fph, Sigi)),1) ci = ci+1 } } colnames(fpi) = c("hospid", "empid", "fpi1") dta = as.data.frame(fpi) lme = lme(fpi1 ~ 1, dta, ~1|hospid) summary(lme) lme = lme(fpi1 ~ 1, dta, ~1|hospid) summary(lme) Linear mixed-effects model fit by REML Data: dta AIC BIC logLik -2555.416 -2534.771 1280.708 Random effects: Formula: ~1 | hospid (Intercept) Residual StdDev: 0.06173257 0.1997302 Fixed effects: fpi1 ~ 1 Value Std.Error DF t-value p-value (Intercept) 0.368082 0.006919828 7110 53.19236 0 Standardized Within-Group Residuals: Min Q1Med Q3Max -3.4870696 -0.6747173 -0.0048658 0.6838012 4.2633384 Number of Observations: 7200 Number of Groups: 90 0.0617^2 [1] 0.00380689 __ 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] lme cant get parameter estimated correctly
csiro.au> writes: > > Here is a demo you may like to consider. (I can see what you are trying > to do with your loops, but I prefer to do it this way.) This is just for pedagogical purpose, so I like to keep the simple-minded 'for' loop. But what I really wonder is why do I not get the right estimate? It cant be because of how I loop-create the data? Thanks Toby > > On 32 bit Windows, (which I am forced to use), your seed is not a valid > integer, so I have changed it to something which is. > > > set.seed(7658943) __ 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] bugs() ignores my inits
Hi All I can specify whatever inits, it has no effect on the estimation. I am replicating a textbook example. The result is completely trash, having estimates of -58.7 (sd=59.3), where it should be closer to an ml estimate of 0.585 (SE=0.063). The two chains within one run are different, but with different inits for different runs, I get exactly the same chains, and I mean exactly. If I set strong (high-precision) priors using the ML estimates, I get the result I want. I dont want to set strong priors, but I want to set reasonable inits; ones that are not taken by sampling from the prior, but ones that I! specify. Thanks Toby data <- list("n", "n3", "y", "person", "occ") model <- function() { ... } write.model(model, "example.txt") inits <- list(list(m1=0.585, m2=0.718, m3=0.672, m4=0.639),list(m1=0.585, m2=0.718, m3=0.672, m4=0.639)) parameters <- c("m1", "m2", "m3", "m4", "sig21", "sig22") bugs1 <- bugs(data, inits, parameters, "example.txt", n.chains=2, n.iter=15000, clearWD=TRUE, debug=1, DIC=0) display(log) check(J:/project/ps/code/example.txt) model is syntactically correct data(J:/project/ps/code/data.txt) data loaded compile(2) model compiled inits(1,J:/project/ps/code/inits1.txt) this chain contains uninitialized variables inits(2,J:/project/ps/code/inits2.txt) this chain contains uninitialized variables gen.inits() initial values generated, model initialized thin.updater(15) update(500) set(m1) set(m2) set(m3) set(m4) set(sig21) set(sig22) update(500) coda(*,J:/project/ps/code/coda) stats(*) Node statistics nodemeansd MC error 2.5%median 97.5% start sample m1 -58.66 59.27 9.001 -120.8 0.4013 0.7102 501 1000 m2 -58.53 59.27 9.001 -120.7 0.5417 0.8314 501 1000 m3 -58.58 59.27 9.0 -120.7 0.4783 0.7871 501 1000 m4 -58.6 59.27 9.0 -120.7 0.3555 0.7695 501 1000 sig21 0.07409 0.01034 3.039E-40.05636 0.07319 0.09558 501 1000 sig22 7206.0 7500.0 1097.0 0.07909 8436.0 19810.0 501 1000 history(*,J:/project/ps/code/history.odc) History save(J:/project/ps/code/log.odc) save(J:/project/ps/code/log.txt) __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Error bars using data from bugs()
only one out of many many dt = cbind(sim$mean$u2+sim$mean$beta1,sim$sd$u2) dt = dt[order(dt[,1]),] bounds = cbind(c(dt[,1]-1.96*dt[,2],dt[,1]+1.96*dt[,2],dt[,1]),rep(1:length(sim$mean$u2),3)) bounds = bounds[order(bounds[,2]),] plot(bounds) T Matthew Krachey wrote: > I'm trying to compare the results of several models using output from bugs(). > I want to summarize the models via boxplots or lineplots of the parameters. I > know how to use errorbars from a regression, but how do I implement the > credible intervals as error bars into these plots > > Any help is much appreciated > > Matthew Krachey > > > [[alternative HTML version deleted]] > > __ > [EMAIL PROTECTED] 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] lm design matrix bug?
Hi All Maybe I dont understand it, but I would have expected that the design matrix has as many rows as there were observations available to fit the model. Below a small artificial dataset created, then one model fitted and the design matrix outputted, having 27 rows. Then I delete 6 obs, and fit the model on these 21 obs, but the design matrix that comes out has 26 rows? Thanks for your enlightenment. Toby y = c() x1 = c() x2 = c() idx = 1 for (i in 1:3) { for (j in 1:3) { for (k in 1:3) { y[idx] = 30*i+10*j+100*i*j+30*k-60 x1[idx] = i x2[idx] = j idx = idx+1 } } } lm11 = lm(y ~ factor(x1)*factor(x2), x=1) summary(lm11) unique(predict(lm11)) X = lm11$x; X P = solve(t(X)%*%X) %*% t(X); round(P,3) y[3] = NA y[6] = NA y[12] = NA y[18] = NA y[24] = NA y[27] = NA lm21 = lm(y ~ factor(x1)*factor(x2), x=1) summary(lm21) unique(predict(lm21)) X = lm21$x; X P = solve(t(X)%*%X) %*% t(X); round(P,3) __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] lm design matrix bug?
uups, ok, I see, thanks Tim Calkins wrote: > check the dimensions of your X and P matrices: > >> dim(X) > [1] 21 9 >> dim(P) > [1] 9 21 > > what you see is that the final row is *named* 26. It's actually row 21. > > cheers, > > On 10/29/07, [EMAIL PROTECTED] <[EMAIL PROTECTED]> wrote: >> Hi All >> >> Maybe I dont understand it, but I would have expected that the design matrix >> has >> as many rows as there were observations available to fit the model. >> Below a small artificial dataset created, then one model fitted and the >> design >> matrix outputted, having 27 rows. Then I delete 6 obs, and fit the model on >> these 21 obs, but the design matrix that comes out has 26 rows? >> >> Thanks for your enlightenment. >> >> Toby >> >> >> >> >> y = c() >> x1 = c() >> x2 = c() >> idx = 1 >> for (i in 1:3) { >> for (j in 1:3) { >>for (k in 1:3) { >> y[idx] = 30*i+10*j+100*i*j+30*k-60 >> x1[idx] = i >> x2[idx] = j >> idx = idx+1 >>} >> } >> } >> >> lm11 = lm(y ~ factor(x1)*factor(x2), x=1) >> summary(lm11) >> unique(predict(lm11)) >> >> X = lm11$x; X >> >> P = solve(t(X)%*%X) %*% t(X); round(P,3) >> >> >> y[3] = NA >> y[6] = NA >> y[12] = NA >> y[18] = NA >> y[24] = NA >> y[27] = NA >> >> lm21 = lm(y ~ factor(x1)*factor(x2), x=1) >> summary(lm21) >> unique(predict(lm21)) >> >> X = lm21$x; X >> >> P = solve(t(X)%*%X) %*% t(X); round(P,3) >> >> __ >> 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.