Hi,
I have a quick statistics-related question. I would like to perform a
comparison of means across treatments using a linear model framework.
However, my 'baseline' treatment has many more observations than either of
the other two treatments and the variance within each treatment is different.
Furthermore, each observation is assigned a weight-- where weights within a
treatment sum to 1.
Can I safely use lm() to check differences in means, under these
circumstances? I have posted an example below.
# generate some data
# note lack of balance
x.1 <- rnorm(n=100, mean=0, sd=1)
x.2 <- rnorm(n=10, mean=4, sd=3)
x.3 <- rnorm(n=10, mean=0, sd=1.5)
# generate some weights that sum to 1 for each treatment
# c/o this post:
# https://stat.ethz.ch/pipermail/r-help/2008-July/167044.html
#
w.1 <- diff(c(0,sort(sample(seq(1,99,1),99,replace=T)),100)) / 100
w.2 <- diff(c(0,sort(sample(seq(1,99,1),9,replace=T)),100)) / 100
w.3 <- diff(c(0,sort(sample(seq(1,99,1),9,replace=T)),100)) / 100
# generate treatment labels:
trt <- factor(c(rep('baseline', 100), rep(c('data1','data2'), each=10)))
# combine into dataframe
d <- data.frame(values=c(x.1,x.2,x.3), wts=c(w.1,w.2,w.3), treatment=trt)
# compute means and wt.means
d.means <- sapply(by(d, d$treatment, function(i) mean(i$value)), '[')
d.wt.means <- sapply(by(d, d$treatment, function(i) weighted.mean(i$value,
w=i$wts)), '[')
# quick check:
boxplot(values~ treatment, data=d, varwidth=TRUE, border=grey(0.5),
main='symbol size is proportional to observation weight')
points(1:3, d.means, col='green', pch=15, cex=1.2)
points(1:3, d.wt.means, col='blue', pch=15, cex=1.2)
points(jitter(as.numeric(d$treatment)), d$values, col='red',
cex=sqrt(d$wts*100))
legend('topleft', legend=c('original data','mean','wt.mean'), pch=c(1,15,15),
col=c('red','green','blue'))
# comparison of means (without weights):
# effects are equal to treatment means
summary(lm(values ~ treatment, data=d))
# comparison with means
# effects are equal to treatment weighted-means
summary(lm(values ~ treatment, data=d, weights=wts))
Thanks,
Dylan