[R] Indexing within by statement - different coloured lines in abline wanted..
Dear R-list
I'm trying to get each regression line, plotted using abline, to be of a
different colour as the following code illustrates. I'm hoping there is a
simple indexing solution. Many thanks.
## code from here
colours=c(black,red,blue,green,pink)
Mean=500;Sd=10;NosSites=5;Xaxis=seq(1,5,1)
SlopeCoefficient=5;Site=(gl(NosSites,length(Xaxis),labels=1:NosSites))
Predictor=rep(Xaxis,NosSites)
InterceptAdjustment=rnorm(n=NosSites,mean=Xaxis,sd=50)
RandomIntercept=rep(InterceptAdjustment,each=length(Xaxis))
PreResponse=rnorm(n=length(Predictor),
mean=Mean+SlopeCoefficient*1:length(Xaxis),sd=Sd)
Response1=PreResponse+RandomIntercept
#create data frame
Data2=data.frame(Site,Predictor,Mean,SlopeCoefficient,RandomIntercept,Response1)
Data1=data.frame(Site=Data2$Site,Predictor=Data2$Predictor,Response1=Data2$Response1)
#plotting
var=as.numeric(levels(Data1$Site))
par(mfrow=c(1,3))
plot(Response1~Predictor,data=Data1,xlim=c(min(Xaxis),max(Xaxis)),ylim=c(MN,MX),
pch=as.numeric(Site),main=Raw data with linear regresssions by Site)
by(Data1,Data1$Site,function(Site){
par(new=T)
abline(lm(Response1~Predictor,data=Site),col=colours[])#index in here.
})
The Scottish Association for Marine Science (SAMS) is registered in Scotland as
a Company Limited by Guarantee (SC009292) and is a registered charity (9206).
SAMS has an actively trading wholly owned subsidiary company: SAMS Research
Services Ltd a Limited Company (SC224404). All Companies in the group are
registered in Scotland and share a registered office at Scottish Marine
Institute, Oban Argyll PA37 1QA. The content of this message may contain
personal views which are not the views of SAMS unless specifically stated.
Please note that all email traffic is monitored for purposes of security and
spam filtering. As such individual emails may be examined in more detail.
[[alternative HTML version deleted]]
__
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] Power analysis in hierarchical models
Dear All I am attempting some power analyses, based on simulated data. My experimental set up is thus: Bleach: main effect, three levels (control, med, high), Fixed. Temp: main effect, two levels (cold, hot), Fixed. Main effect interactions, six levels (fixed) For each main-effect combination I have three replicates. Within each replicate I can take varying numbers of measurements (response variable = Growth (of marine worms)) but, for this example, assume eight). (I’m interested in changing this to see if the experimental power changes much). Total size = 3 x 2 x 3 x 8 = 144 The script thus far goes: === start of script = library(lme4) #Data frame structure Bleach=rep(c(Control,Med,High),each=48) Temp= rep(rep(c(Cold,Hot),each=24),3) Rep= (rep(rep(rep(c(1,2,3),each=8),2),3)) Ind= (rep(rep(rep(c(1:8),3),2),3))#not required for stats #Fake data (based on pilot studies), only showing a single main effect (bleach) Growth=c( rnorm(48,3.27,0.77),rnorm(48,3.21,0.77),rnorm(48,3.64,1.17)) fake2=data.frame(Bleach,Temp,Rep,Ind,Growth);head(fake2) #generate factor level for lmer as per Crawley, page 649 fake2$rep=fake2$Bleach:fake2$Temp:fake2$Rep#rep is used in the lmer model with(fake2,table(rep))#check that each rep contains 8 measurements # run alternative (?equivalent) models model1=aov(Growth~Bleach*Temp+Error(Bleach*Temp/Rep),data=fake2);summary(model1) model2=lmer(Growth~Bleach*Temp+(1|rep),data=fake2);summary(model2)#note: see above, repRep! end of script == I'd like to get familiar with using lme4 because it is likely that the final results of the experiment will be unbalanced (which precludes the use of aov I think). The df given by model1 seem to make sense. Any guidance on any of the following would be much appreciated: 1. Are model1 and model2 equivalent? 2. For model1 - is the random component correctly specified and is there a (simple) mechanism to get the appropriate F ratios and P values? 3. For model2 - again, is the random component correct (probably not) and why is the random effect (rep) variance and standard deviations so low (zero in most iterations)? 4. For both models - how do I isolate (so I can tabulate and create histograms) the appropriate P and/or t values? (for model2 - the ‘mer’ object doesn’t seem to contain the t values but maybe I’m missing something). Direction to any more generic sources of information regarding power analysis in hierarchical models would be gladly received. Thank you Tom. - Tom Wilding, MSc, PhD, Dip. Stat. Scottish Association for Marine Science, Scottish Marine Institute, OBAN Argyll. PA37 1QA United Kingdom. Phone (+44) (0) 1631 559214 Fax (+44) (0) 1631 559001 +++ The Scottish Association for Marine Science (SAMS) is registered in Scotland as a Company Limited by Guarantee (SC009292) and is a registered charity (9206). SAMS has an actively trading wholly owned subsidiary company: SAMS Research Services Ltd a Limited Company (SC224404). All Companies in the group are registered in Scotland and share a registered office at Scottish Marine Institute, Oban Argyll PA37 1QA. The content of this message may contain personal views which are not the views of SAMS unless specifically stated. Please note that all email traffic is monitored for purposes of security and spam filtering. As such individual emails may be examined in more detail. +++ __ 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] Multifactor boxplots
Dear All I would like to produce interaction boxplots and this seems to work: par(mfrow=c(2,2)) A=sample(rnorm(50,50,10)) B=sample(rnorm(50,100,10)) Test=merge(A,B,by=0)#by=0 where 0 is the row.names TreatA=(gl(2,50,100,labels=c(High,Low))) TreatB=rep(gl(2,25,50,labels=c(High,Low)),2) Newdata=data.frame(TreatA,TreatB,Test) bwplot(x~TreatA:TreatB,data=Newdata) However, I would prefer the X axis labels to be different, such that there are two label rows (TreatA and TreatB) something like this: TreatAHigh HighLow Low TreatBHigh Low High Low any guidance on achieving this much appreciated. Regards Tom [[alternative HTML version deleted]] __ 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] Random selection from a subsample
Dear Mailing List
I have a data set (data4) consisting of a number of factors and a response
variable. I wish to randomly sample from a combination of two of those factors
(GIS_station and Distance_code2) and return a new dataframe containing the
original data structure (i.e. all the columns) but only containing the randomly
selected rows. The number of rows in each combination of GIS_station and
Distance_code2 vary (widely) and some combinations are absent.
This is getting there::
with (data4,{
sub_sample10=by(data4,list(GIS_station,Distance_code2), function(x)
{sample(1:nrow(x),10,replace=T)})
})
but just generates two random numbers from the range 1:nrow(x). It doesn't
return the selected rows, which is what I want.
I'm sure I could this could be done in an elegant manner, using a subscript e.g.
sub_sample10 = data4 [sample (1:nrow (data4), size=10), ]
only somehow combining it with the 'by' statement (e.g. by (data4, list
(GIS_station, Distance_code2)...)) but I cannot get this to work.
Any guidance on this much appreciated.
Thankyou.
__
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.
