[R-sig-ME] Testing differences in measurement variance

Mike Lawrence Mike.Lawrence at dal.ca
Fri Aug 13 00:23:57 CEST 2010


On Thu, Aug 12, 2010 at 6:48 PM, Ben Bolker <bbolker at gmail.com> wrote:
> Var(group i) = Var(process error) + Var(measurement error, group i)


Ah, yes, I think the assumption of equal "measurement error" across
groups but unequal "process error" is appropriate for the question I
raised. To be as clear as possible, I've included below an example of
the type of data I'd anticipate encountering. There are two groups of
individuals, both groups have the same measurement error but have
their own unique process error added onto this measurement error.
Individuals within each group will of course not manifest precisely
the group's base variability but vary from the base variability by
some random amount.

set.seed(1)

measurement_error = 1
group_A_base_sd = measurement_error + 1
group_B_base_sd = measurement_error + 2
within_group_sd_of_sds = .1

n_per_group = 10
obs_per_id = 20

temp = data.frame(
	id = 1:(n_per_group*2)
	, group = rep(c('A','B'))
)

#generate example data
library(plyr) #to avoid loops (for coding convenience only)
obs_data = ddply(
	.data = temp
	, .variables = .(id,group)
	, .fun = function(x){
		#generate a unique sd for this individual
		# based on their group's sd plus some
		# within-group variability
		id_sd = ifelse(
			x$group=='A'
			, rnorm(
				1
				, group_A_base_sd
				, within_group_sd_of_sds
			)
			, rnorm(
				1
				, group_B_base_sd
				, within_group_sd_of_sds
			)
		)
		#generate data points with the above generated
		# variability
		to_return = data.frame(
			obs_num = 1:obs_per_id
			, measurement = rnorm(obs_per_id,0,id_sd)
		)
		return(to_return)
	}
)

#first step of an anova-based approach:
# compute SDs within each Ss
obs_sds = ddply(
	.data = obs_data
	, .variables = .(id,group)
	, .fun = function(x){
		to_return = data.frame(
			obs_sd = sd(x$measurement)
		)
	}
)

#second step of an anova-based approach:
# compute the anova on the SDs
summary(
	aov(
		formula = obs_sd~group
		, data = obs_sds
	)
)




-- 
Mike Lawrence
Graduate Student
Department of Psychology
Dalhousie University

Looking to arrange a meeting? Check my public calendar:
http://tr.im/mikes_public_calendar

~ Certainty is folly... I think. ~




More information about the R-sig-mixed-models mailing list