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)
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))
yv<-predict(model_galls,list(psyllids =xv),type="response")
lines(xv,yv,lwd=2,col=4)