Dear all,
Does anyone have any suggestions on how to block multiple columns of data
one at a time in the midst of an analysis, having specified the blocking
variable?
I am running a random effects model using lmer, and my data set has
multiple columns.
Individuals in the study are grouped into 60 families- which is the
blocking factor.
The random effects are markers (labeled Pxlyy below) that have two levels-
1 (presence) or 0 (absence).
I need to estimate Variance Components for 'peg.no' (family wise, which is
why this is the blocking factor), explained by each PxLyy, and by 'family'
as well.
I have managed to run the lmer model for the first marker column (P1L55),
and have been able to block it according to 'family' as follows:
vc<-read.table("P:\\R\\Testvcomp10.txt",header=T)
attach(vc)
family<-factor(family)
colms<-(vc)[4:13] ## these are the marker columns
lapply(colms,factor)
library(lme4)
for (f in levels(family)) {
print("Family")
print(f)
{
t<- vc [ family == f, ]
try(lmer<-fit(peg.no~1 + (1|P1L55) + (1|family), data=t,
na.action=na.exclude))
print(summary(fit))
}}
However I need to do this for 10 PxLyy columns simultaneously, so that I
have 60 families * 10 markers= 600 runs of the model.
To do this I tried,
for (c in 1:length (colms))
{
for (f in levels (family))
{
print("Marker")
print(c)
print ("Family")
print(f)
{
t<- vc[family == f,] ## this is where I lose the plot
try(fit<-lmer(peg.no~1 + (1|c)+ (1|family), data=t, na.action=na.exclude))
print(summary(fit))
}}}
This however, doesn't run at all after the first two loops.
I cannot figure out how to tell it to block each one of 10 columns in turn.
This works for one marker, but not more in the same run.
Any suggestions please?
Thanks a lot,
Aditi
Data files at:
<http://www.4shared.com/file/131980362/460bdafe/Testvcomp10.html> (excel)
http://www.4shared.com/file/131980512/dc7308b/Testvcomp10.html
(txt)
---------------------
A Singh
aditi.si...@bristol.ac.uk
School of Biological Sciences
University of Bristol
______________________________________________
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.
