cerana<-c(1,1,1,1,1,1,1,1,1,1,0,1,0,1,1,1,0,1,0,1,1,1,0,1,1,0,1,1,0,0,0,1,0,1,0,0,0,0,1,1,0,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,0,1,0,0,1,1,0,0,0,1,0,1,1,0,0,1,1,1,0,1,1,1,1,0,1,1,1)
dorsata<-c(0,1,1,1,0,0,1,1,0,0,0,0,0,1,0,1,1,1,0,0,1,0,0,0,0,1,0,1,0,0,1,0,0,1,1,1,1,0,1,0,1,1,1,0,1,0,1,0,0,1,1,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0,1, 0,0,0,0,1,0,0,0,0,1,0,0,0,1,0,1,0,1,0,0,1,0)
florea<-c(0,0,0,1,1,1,0,1,0,1,1,1,0,0,1,0,0,0,0,1,1,0,1,1,0,0,0,1,1,1,0,0,0,1,0,0,0,0,1,0,0,1,1,0,1,1,0,0,1,0,0,0,1,1,0,0,0,1,1,1,1,1,1,1,0,0,0,1,1,0,0,0,0,1,1,0,1,0,0,1,0,1,0,0,1,1,0,0,1,1,1,0,0,1,0)
mellifera<-c(0,1,0,1,0,0,0,1,1,1,0,1,1,0,0,0,1,1,1,1,1,1,1,1,0,0,1,0,0,0,0,1,1,1,1,0,0,0,1,1,0,0,1,0,1,1,1,0,1,0,0,1,1,1,0,0,1,0,0,1,0,1,1,1,1,1,1,0,0,1,1,1,0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,1,0,0,0,0,1,1)
data<-cbind(cerana,dorsata,florea,mellifera)
raw_score= sum(rowSums(data[,1:4])^2)
raw_score
n<-nrow(data) # getting the number of plant species
randomisations<-NULL
for(i in 1:100000){
d1<-sample(data[,1],n,replace=F)
d2<-sample(data[,2],n,replace=F)
d3<-sample(data[,3],n,replace=F)
d4<-sample(data[,4],n,replace=F)
shuffled_score<-sum((d1+d2+d3+d4)^2)
randomisations<-c(randomisations,shuffled_score)
}
hist(randomisations)
abline(v=raw_score,lwd=2,col=2)
text(528,12000,"observed",srt=270, col="red")
equals<-length(which(randomisations==raw_score))
below<-length(which(randomisations<raw_score))+0.5*equals
above<-length(which(randomisations>raw_score))+0.5*equals
text(480,16000,paste("less than = ",below,sep=""))
text(560,16000,paste("greater than = ",above,sep=""))
text(560,12000,paste("P = ",round(above/(above+below),4)), col="blue")
psyllids<-c(0:9)
no_Bracon<-c(5,139,153,109,70,28,20,2,3,1)
Bracon<-c(63,57,61,36,20,9,3,0,0,0)
# The proportion containing a Bracon can be calculated by dividing by the total Bracon+no_Bracon
PropBracon<-Bracon/(Bracon+no_Bracon)
plot(psyllids,PropBracon, pch=15, xlab= "Number of psyllids in gall",ylab="Proportion of galls containing a Bracon")
data<- cbind(Bracon, no_Bracon)
model_galls <- glm(data~psyllids, family=binomial) # we can omit “Family=”
summary(model_galls)
xv<-seq(0,9,0.01)
yv<-predict(model_galls,list(psyllids =xv),type="response")
lines(xv,yv,lwd=2,col=4)
nPsyllid<-NULL
Bracon_status<-NULL
for (i in 1:10){
np<-i-1 #because the first row in data corresponds to zero psyllids
if(data[i,1] >=1){
for (j in 1:data[i,1]){
nPsyllid <- c(nPsyllid,np)
Bracon_status<-c(Bracon_status, 0)}}
if(data[i,2] >=1){
for (j in 1:data[i,2]){
nPsyllid <- c(nPsyllid,np)
Bracon_status<-c(Bracon_status, 1)}}
}
xpandedData<-cbind(nPsyllid,Bracon_status)
plot(jitter(xpandedData[,1],1,0.15), jitter(xpandedData[,2],1,0.03), ylab="Bracon in gall", xlab="Total number of psyllids in gall", main="", pch=20, cex=0.3, axes=F)
axis(1,c(0:9),at=c(0:9))
axis(2,c("absent","present"),at=c(0,1))
i
n
View(colon)
a<-1:10
v<
-a°1.7
b<-a^2
plota(a,b)
plot(a,b)