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