plots.resid.regression.lines.R 4.6 KB
Newer Older
1
2
3
4
5
6
7
##### SCRIPT TO TEST plots resuls with regression lines
source("R/analysis/lmer.run.R")

output.dir.lmer <- function (model, trait, set, ecoregion,type.filling) {
  file.path("output/lmer", set,ecoregion,trait,type.filling,model)
}

8
9
10
11
12
13
14
15
16
17
18
19
20
21
fun.plot.colors.density <- function(x,y,...){
mat   <- cbind(x,y)
data <- as.data.frame(mat)
colors.dens  <- densCols(mat)
plot(x,y, col=colors.dens, pch=20,...)
}

fun.points.colors.density <- function(x,y,...){
mat   <- cbind(x,y)
data <- as.data.frame(mat)
colors.dens  <- densCols(mat)
points(x,y, col=colors.dens, pch=20,...)
}

Georges Kunstler's avatar
Georges Kunstler committed
22
23
24
25
26
fun.boxplot.breaks <-  function(x,y,Nclass=15,add=TRUE,...){
x1 <- x[!is.na(x) & !is.na(y)] 
y1 <- y[!is.na(x) & !is.na(y)] 
x <- x1
y <- y1
27
28
breakss <- seq(min(x),max(x),length=Nclass+1)
x.cut <- cut(x,breakss)
Georges Kunstler's avatar
Georges Kunstler committed
29
30
mid.points <- breakss[-(Nclass+1)]+abs(breakss[2]-breakss[1])/2
boxplot(y~x.cut,at=mid.points,names=round(mid.points,0),outline=FALSE,add=add,col=grey(1-table(x.cut)/length(x.cut)),...)
31
32
33
34
35
36
37
}

seq.from.to.quantile <- function(x,length.out=20,probs=c(0.001,0.999)){
qq <- quantile(x,probs)
return(seq(from=qq[1],to=qq[2],length.out=length.out))
}

Georges Kunstler's avatar
Georges Kunstler committed
38
fun.load.data.for.partial.residual <- function(trait,set,ecoregion,type.filling){
39
40
df.lmer <- load.and.prepare.data.for.lmer(trait,set,ecoregion,
                                           min.obs=10, sample.size=NA,
41
                                         type.filling,std=FALSE) # return a DF
42
43
44
45
46
47
ERcomp <- readRDS(file.path("output/lmer", set,ecoregion,trait,
                            type.filling,"lmer.LOGLIN.ER.Tf",
                            "results.rds"))
ADcomp <- readRDS(file.path("output/lmer", set,ecoregion,trait,
                            type.filling,"lmer.LOGLIN.AD.Tf",
                            "results.rds"))
Georges Kunstler's avatar
Georges Kunstler committed
48
return(list(df.lmer=df.lmer,
49
50
            ERcomp=ERcomp,ADcomp=ADcomp,trait=trait,
            set=set,ecoregion=ecoregion,
Georges Kunstler's avatar
Georges Kunstler committed
51
52
            type.filling=type.filling,
            set=set,ecoregion=ecoregion))
53
54
}

Georges Kunstler's avatar
Georges Kunstler committed
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
fun.compute.partial.residual.for.one.var <- function(df.lmer,ERcomp,ADcomp,var,var.vec,set,ecoregion){
df.lmer <- df.lmer[!is.na(df.lmer[[var]]),]    
df.lmer <- df.lmer[sort(df.lmer[[var]],index.return=TRUE)$ix,]
df.lmer2 <- df.lmer
df.lmer2[[var]] <- 0
df.lmer3 <- df.lmer
for (i in var.vec[!var.vec %in% var]) df.lmer3[[i]] <- mean(df.lmer[[i]])
df.lmer4 <- df.lmer3
df.lmer4[[var]] <- 0
# compute partial residual
partial.resid <- df.lmer[['logG']] -predict(ERcomp,newdata=df.lmer2,REform=NULL)
predict.line <- predict(ERcomp,newdata=df.lmer3,REform=NA)-predict(ERcomp,newdata=df.lmer4,REform=NA)
par(mar=c(4.5,3.5,0.5,0.5),mgp=c(1.5,0.5,0))
fun.boxplot.breaks(df.lmer[[var]],partial.resid,add=FALSE,xlab=var,
                   ylab="growth residual",
                   main=paste("Effect ",set,ecoregion,sep=""))
lines(df.lmer[[var]],predict.line)
}    
73
74


75
76
77
78
79
80
81

## function to get all set and ecoregion to plot

get.ecoregions.for.set <- function(set, base.dir = "./output/processed/"){
  sub(paste(base.dir,set,"/",sep=""),"",list.dirs(paste(base.dir,set,sep="")))[-1]
}

Georges Kunstler's avatar
Georges Kunstler committed
82
83
84
85
86
87
88
89
fun.plot.partial.resid.var <- function(var.name,set,trait,ecoregions,type.filling,
                                       var.vec=c('logD','sumBn','Tf','sumTfBn','sumTnBn')){
pdf(paste("figs/plot.resid/",
          paste(trait,set,type.filling,"partial.residual.per.ecocode",var.name,"pdf",sep="."),sep=""),width=15,height=8)
par(mfrow=c(2,ceiling(length(ecoregions)/2)))
          for (ecoregion in ecoregions){
try(list.data <- fun.load.data.for.partial.residual(trait,set,ecoregion,type.filling))
try(with(list.data,fun.compute.partial.residual.for.one.var(df.lmer,ERcomp,ADcomp,var=var.name,var.vec=var.vec,set,ecoregion)) )
90
         }
Georges Kunstler's avatar
Georges Kunstler committed
91
dev.off()
92
93
}

Georges Kunstler's avatar
Georges Kunstler committed
94
fun.plot.partial.resid.all <-  function(set,
95
               traits = c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass"),
Georges Kunstler's avatar
Georges Kunstler committed
96
               type.fillings=c("species") ){
97
ecoregions <- get.ecoregions.for.set(set, base.dir = "./output/processed/")    
Georges Kunstler's avatar
Georges Kunstler committed
98
    for(trait in traits){
99
      for (type.filling in type.fillings){
Georges Kunstler's avatar
Georges Kunstler committed
100
101
102
103
104
105
106
107
108
               fun.plot.partial.resid.var(var.name='sumTfBn',set,trait,
                                          ecoregions,type.filling,
                           var.vec=c('logD','sumBn','Tf','sumTfBn','sumTnBn'))

               fun.plot.partial.resid.var(var.name='sumTnBn',set,trait,
                                          ecoregions,type.filling,
                           var.vec=c('logD','sumBn','Tf','sumTfBn','sumTnBn'))
      }   
    }
109
110
111
}


Georges Kunstler's avatar
Georges Kunstler committed
112
113
114
#### DO ALL THE PLOTS
      
fun.plot.partial.resid.all('US')
115

Georges Kunstler's avatar
Georges Kunstler committed
116
117
for (set in c("Sweden","NSW","Mbaiki","Luquillo","Japan")){#"BCI","Canada","France","Fushan","NVS","Paracou","Spain","US","Swiss",
fun.plot.partial.resid.all(set)
118
119
120
}