


cas<- read.csv(file="cas.csv",header=TRUE,sep=",")

#############

model3<-lme(meansize~year*stock,random=~year|id, data=cas) # the base model, the most simple for my objective.



model3.1<-lme(meansize~year*stock,random=~year|id, data=cas,weights=varFixed(~1/w)) # I want to account for the fact that some gear are more important because they catch more fish. so i want to give more weights to the gears that catch more fish and less weights to the gears that catch less fish.


# compare both models
anova(model3,model3.1)

# explore models
plot(model3, form=resid(., type="p")~fitted(.),abline=0)
plot(model3, form=resid(., type="p")~fitted(.)| stock,abline=0)
plot(model3.1, form=resid(., type="p")~fitted(.),abline=0)
plot(model3.1, form=resid(., type="p")~fitted(.)| stock,abline=0)
qqnorm(model3,~resid(.) | stock,id=0.05,adj=-0.3)
qqnorm(model3.1,~resid(.) | stock,id=0.05,adj=-0.3)

#expand the models to take into account for the fact that different groups(specis) have different variances

#option 1
model3.6<-lme(meansize~year*stock,random=~year|id, data=cas,weights=varIdent(form=~1|stock),
control=list(maxIter=500, msMaxIter=500, niterEM=500))

#option 2
model3.9<-lme(meansize~year*stock,random=~year|id, data=cas,weights=varComb(varFixed(~1/w),varIdent(form=~1|stock)),control=list(maxIter=500, msMaxIter=500, niterEM=500))


#explore the models.
plot(model3.6, form=resid(., type="p")~fitted(.)| stock,abline=0)
plot(model3.9, form=resid(., type="p")~fitted(.)| stock,abline=0)

qqnorm(resid(model3.6,type="p"),cex=1,main="")
qqnorm(resid(model3.9,type="p"),cex=1,main="")

qqnorm(model3.6,~resid(.) | stock,id=0.05,adj=-0.3)
qqnorm(model3.9,~resid(.) | stock,id=0.05,adj=-0.3)

## compare models

anova(model3.6,model3.9)
         Model df      AIC      BIC    logLik
model3.6     1 12 5893.051 5950.667 -2934.526
model3.9     2 12 7599.363 7656.978 -3787.681

# if i compare the four models

anova(model3, model3.1,model3.6,model3.9)
         Model df      AIC      BIC    logLik   Test  L.Ratio p-value
model3       1 10 6532.240 6580.253 -3256.120                        
model3.1     2 10 7876.156 7924.169 -3928.078                        
model3.6     3 12 5893.051 5950.667 -2934.526 2 vs 3 1987.104  <.0001
model3.9     4 12 7599.363 7656.978 -3787.681  




