lmer.run.nolog.all.R 7.82 KB
Newer Older
Georges Kunstler's avatar
Georges Kunstler committed
1
2
3
4
5
6
7
8
9
10
11
12
13
###########################
###########################
### FUNCTION TO RUN LMER ESTIMATION WITH NO logBA for all in one big model
library(lme4)



run.models.for.set.all.traits  <- function(model.file,fun.model,  traits =
         c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass"),type.filling, ...){
  for(trait in traits)
   run.multiple.model.for.set.one.trait(model.file,fun.model, trait,  type.filling=type.filling,...)
}

Georges Kunstler's avatar
Georges Kunstler committed
14
run.multiple.model.for.set.one.trait <- function(model.files,fun.model, trait, type.filling, ...){
Georges Kunstler's avatar
Georges Kunstler committed
15
  for (m in model.files)
Georges Kunstler's avatar
Georges Kunstler committed
16
    (run.model.for.set.one.trait (m, fun.model,trait, type.filling=type.filling, ...))
Georges Kunstler's avatar
Georges Kunstler committed
17
18
19
}


Georges Kunstler's avatar
Georges Kunstler committed
20
run.model.for.set.one.trait <- function(model.file,fun.model, trait, set,type.filling, ...){
Georges Kunstler's avatar
Georges Kunstler committed
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
    fun.model <- match.fun(fun.model)
    try(fun.model(model.file, trait,  type.filling=type.filling,...))
}


#=====================================================================
# Run lmer() (in package lme4) for one ecoregion, one trait, one model
#=====================================================================

model.files.lmer.Tf.1 <- c("R/analysis/model.lmer.all/model.lmer.LOGLIN.E.Tf.R",
                 "R/analysis/model.lmer.all/model.lmer.LOGLIN.R.Tf.R",
                 "R/analysis/model.lmer.all/model.lmer.LOGLIN.ER.Tf.R")
model.files.lmer.Tf.2 <- c("R/analysis/model.lmer.all/model.lmer.LOGLIN.nocomp.Tf.R",
                 "R/analysis/model.lmer.all/model.lmer.LOGLIN.AD.Tf.R",
                 "R/analysis/model.lmer.all/model.lmer.LOGLIN.HD.Tf.R",
                 "R/analysis/model.lmer.all/model.lmer.LOGLIN.simplecomp.Tf.R")



fun.call.lmer.and.save <- function(formula,df.lmer,path.out){
41
    lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE,
Georges Kunstler's avatar
Georges Kunstler committed
42
                       control=lmerControl(optCtrl=list(maxfun=40000) ) )
Georges Kunstler's avatar
Georges Kunstler committed
43
44
45
46
   print(summary(lmer.output))
  saveRDS(lmer.output,file=file.path(path.out, "results.nolog.all.rds"))
}

Georges Kunstler's avatar
Georges Kunstler committed
47
run.lmer <- function (model.file, trait,
Georges Kunstler's avatar
Georges Kunstler committed
48
49
50
51
52
53
54
55
56
57
58
59
                      min.obs=10, sample.size=NA,
                      type.filling) {
    require(lme4)
    source(model.file, local = TRUE)
    model <- load.model()
    #= Path for output
    path.out <- output.dir.lmer(model$name, trait, 'all',
                                type.filling=type.filling)
    print(path.out)
    dir.create(path.out, recursive=TRUE, showWarnings=FALSE)
    cat("run lmer for model",model.file," for trait",
         trait,"\n")
Georges Kunstler's avatar
Georges Kunstler committed
60
      df.lmer <- load.and.prepare.data.for.lmer(trait,
Georges Kunstler's avatar
Georges Kunstler committed
61
62
63
64
65
                                           min.obs, sample.size,
                                         type.filling=type.filling) # return a DF
     cat("Ok data with Nobs",nrow(df.lmer),
         "\n")
        #= Run model
Georges Kunstler's avatar
Georges Kunstler committed
66
        fun.call.lmer.and.save(formula=model$lmer.formula.tree.id,df.lmer,path.out)
Georges Kunstler's avatar
Georges Kunstler committed
67
68
69
70
71
72
73
74
75
76
77
78
79
}


#========================================================================

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


#============================================================
# Function to prepare data for lmer
#============================================================
Georges Kunstler's avatar
Georges Kunstler committed
80
81
load.and.prepare.data.for.lmer <- function(trait,
                                  min.obs, sample.size, type.filling,
Georges Kunstler's avatar
Georges Kunstler committed
82
83
                                  base.dir = "output/processed/"){
    ### load data
Georges Kunstler's avatar
Georges Kunstler committed
84
##     library(data.table)
Georges Kunstler's avatar
Georges Kunstler committed
85

Georges Kunstler's avatar
Georges Kunstler committed
86
87
88
89
90
91
92
93
94
95
96
97
## data.tree.tot <- fread(file.path(base.dir,"data.all.csv"),
##                               stringsAsFactors = FALSE)
data.all.sample <- read.csv(file=file.path(base.dir, "data.all.csv"),
                       stringsAsFactors=FALSE,nrows=1000)
classes <- sapply(data.all.sample,class)
classes[classes=='integer'] <- "numeric"
nrows <- as.numeric(system('wc -l < output/processed/data.all.csv',intern=TRUE))

data.tree.tot <- read.csv(file=file.path(base.dir, "data.all.csv"),
                     stringsAsFactors=FALSE,
                     nrows=nrows,
                     colClasses=classes)
Georges Kunstler's avatar
Georges Kunstler committed
98
99
100
101
102

    fun.data.for.lmer(data.tree.tot,trait,type.filling=type.filling)
}

fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.neigh,BATOT,min.obs,
Georges Kunstler's avatar
Georges Kunstler committed
103
                                         min.perc.neigh=0.90,min.BA.G=-60,max.BA.G=500){
Georges Kunstler's avatar
Georges Kunstler committed
104
105
106
107
## remove tree with NA
data.tree <- subset(data.tree,subset=(!is.na(data.tree[["BA.G"]])) &
                                     (!is.na(data.tree[["D"]])) )
## remove tree with zero
108
109
data.tree <- subset(data.tree,subset=data.tree[["BA.G"]]>min.BA.G &
                    data.tree[["BA.G"]]<max.BA.G &
Georges Kunstler's avatar
Georges Kunstler committed
110
                    data.tree[["D"]]>9.9)
Georges Kunstler's avatar
Georges Kunstler committed
111
112
113
114
## select species with no missing traits
select.temp <- (!is.na(data.tree[[abs.CWM.tntf]]) &
                    !is.na(data.tree[[BATOT]]))
data.tree <- data.tree[select.temp,]
Georges Kunstler's avatar
Georges Kunstler committed
115
# select  obs abov min perc.neigh
116
117
data.tree <- subset(data.tree,subset=data.tree[[perc.neigh]] > min.perc.neigh &
                    !is.na(data.tree[[perc.neigh]]))
Georges Kunstler's avatar
Georges Kunstler committed
118
# select  species with minimum obs
Georges Kunstler's avatar
Georges Kunstler committed
119
data.tree <- subset(data.tree,subset=data.tree[["sp"]] %in%
120
121
                    names(table(factor(data.tree[["sp"]])))[
                                   table(factor(data.tree[["sp"]]))>min.obs])
Georges Kunstler's avatar
Georges Kunstler committed
122
123
124
125
126
127
128
129
return(data.tree)
}

fun.center.and.standardized.var <- function(x){
return((x-mean(x))/sd(x))
}

### get variables for lmer
Georges Kunstler's avatar
Georges Kunstler committed
130
131
fun.get.the.variables.for.lmer.tree.id <- function(data.tree,BATOT, CWM.tn,
                                                   abs.CWM.tntf, tf,
Georges Kunstler's avatar
Georges Kunstler committed
132
                                                   min.BA.G=60){
Georges Kunstler's avatar
Georges Kunstler committed
133
logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]]+min.BA.G))
134
135
136
137
logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
species.id <- factor(data.tree[["sp"]])
tree.id <- factor(data.tree[["tree.id"]])
plot.id <- factor(data.tree[["plot"]])
Georges Kunstler's avatar
Georges Kunstler committed
138

139
140
141
142
143
144
145
146
set.id <- factor(data.tree[["set"]])
ecocode.id <- factor(data.tree[['ecocode']])

 #= multiply CWMs by BATOT
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*
                   data.tree[[BATOT]]/max(data.tree[[BATOT]])
sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]/max(data.tree[[BATOT]])
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]/max(data.tree[[BATOT]])
Georges Kunstler's avatar
Georges Kunstler committed
147
sumTnTfBn.diff <- sumTnBn-sumTfBn
148
sumBn <- data.tree[[BATOT]]/max(data.tree[[BATOT]])
Georges Kunstler's avatar
Georges Kunstler committed
149
150
151
152
153
154
return(data.frame(logG= logG,
            logD= logD,
            species.id= species.id,
            tree.id= tree.id,
            set.id= set.id,
            ecocode.id= ecocode.id,
155
            plot.id= plot.id,
Georges Kunstler's avatar
Georges Kunstler committed
156
157
158
159
160
161
            sumTnTfBn.diff= sumTnTfBn.diff,
            sumTnTfBn.abs= sumTnTfBn.abs,
            Tf= data.tree[[tf]],
            sumTnBn= sumTnBn,
            sumTfBn= sumTfBn,
            sumBn= sumBn))
Georges Kunstler's avatar
Georges Kunstler committed
162
163
}

Georges Kunstler's avatar
Georges Kunstler committed
164

Georges Kunstler's avatar
Georges Kunstler committed
165
166
167
168
#============================================================
# Function to prepare data for lmer with new CWM data
# that will be used in a simple linear model
#============================================================
Georges Kunstler's avatar
Georges Kunstler committed
169
170
171
172
173
174
fun.data.for.lmer <-  function(data.tree, trait, min.obs=10,
                               type.filling='species') {
if(! trait %in%  c("SLA", "Leaf.N","Seed.mass",
                   "SLA","Wood.density","Max.height"))
    stop("need trait to be in SLA Leaf.N Seed.mass
         SLA Wood.density or  Max.height ")
Georges Kunstler's avatar
Georges Kunstler committed
175
176
177
178
179
180
# get vars names
CWM.tn <- paste(trait,"CWM",'fill',sep=".")
abs.CWM.tntf <- paste(trait,"abs.CWM",'fill',sep=".")
tf <- paste(trait,"focal",sep=".")
BATOT <- "BATOT"
perc.neigh <- paste(trait,"perc",type.filling,sep=".")
Georges Kunstler's avatar
Georges Kunstler committed
181
182
183
184
185
186
187
188
189
190
191
data.tree <- fun.select.data.for.analysis(data.tree,
                                          abs.CWM.tntf,
                                          perc.neigh,
                                          BATOT,
                                          min.obs)
#= DATA LIST FOR LMER
lmer.data <- fun.get.the.variables.for.lmer.tree.id(data.tree,
                                                    BATOT,
                                                    CWM.tn,
                                                    abs.CWM.tntf,
                                                    tf)
Georges Kunstler's avatar
Georges Kunstler committed
192
193
194
195
196
    return(lmer.data)
}