lmer.output.all-fun.R 3.71 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#### function to analyse lmer output for GLOBAL ANALYSIS


library(lme4)


read.lmer.output <- function(file.name){
  tryCatch(readRDS(file.name), error=function(cond)return(NULL))   # Choose a return value in case of error
}



summarise.lmer.output <- function(x){
 list( nobs = nobs(x),
       R2m =Rsquared.glmm.lmer(x)$Marginal,
       R2c =Rsquared.glmm.lmer(x)$Conditional,
       AIC = AIC(x),
       deviance = deviance(x),
       conv=x@optinfo$conv,
       effect.response.var=variance.fixed.glmm.lmer.effect.and.response(x),
       fixed.coeff.E=fixef(x),
       fixed.coeff.Std.Error=sqrt(diag(vcov(x))),
       fixed.var=variance.fixed.glmm.lmer(x),
       ecocode.BLUP=coef(x)$ecocode.i)
}



summarise.lmer.output.all.list <- function(f ){
    output.lmer <- read.lmer.output(f)
    if(!is.null(output.lmer)){
	res <- list(files.details=files.details.all(f),
		   lmer.summary=summarise.lmer.output( output.lmer))
    }else{
        res <- NULL
    }
    return(res)
}




files.details.all <- function(x){
	s <- data.frame(t(strsplit(x, "/", fixed=TRUE)[[1]]), stringsAsFactors= FALSE)
	names(s)  <- c("d1", "d2", "set", "trait", "filling", "model", "file" )
	s[-(1:2)]
}


#### R squared functions

# Function rsquared.glmm requires models to be input as a list (can include fixed-
# effects only models,but not a good idea to mix models of class "mer" with models
# of class "lme") FROM http://jslefche.wordpress.com/2013/03/13/r2-for-linear-mixed-effects-models/

Rsquared.glmm.lmer <- function(i){
# Get variance of fixed effects by multiplying coefficients by design matrix
      VarF <- var(as.vector(fixef(i) %*% t(i@pp$X)))
      # Get variance of random effects by extracting variance components
      VarRand <- colSums(do.call(rbind,lapply(VarCorr(i),function(j) j[1])))
      # Get residual variance
      VarResid <- attr(VarCorr(i),"sc")^2
      # Calculate marginal R-squared (fixed effects/total variance)
      Rm <- VarF/(VarF+VarRand+VarResid)
      # Calculate conditional R-squared (fixed effects+random effects/total variance)
      Rc <- (VarF+VarRand)/(VarF+VarRand+VarResid)
      # Bind R^2s into a matrix and return with AIC values
      Rsquared.mat <- data.frame(Class=class(i),Family="Gaussian",Marginal=Rm,
                              Conditional=Rc,AIC=AIC(i))
      return(Rsquared.mat)
}


variance.fixed.glmm.lmer <- function(i){

# Get variance of for each fixed effects by multiplying coefficients by design matrix
var.vec <- apply(fixef(i) * t(i@pp$X),MARGIN=1,var)
# Get variance of fixed effects by multiplying coefficients by design matrix
VarF <- var(as.vector(fixef(i) %*% t(i@pp$X)))
# Get variance of random effects by extracting variance components
VarRand <- colSums(do.call(rbind,lapply(VarCorr(i),function(j) j[1])))
# Get residual variance
VarResid <- attr(VarCorr(i),"sc")^2
var.vec <- var.vec/(VarF+VarRand+VarResid)
names(var.vec) <- paste(names(var.vec),"VAR",sep=".")
return(var.vec)
}

variance.fixed.glmm.lmer.effect.and.response <- function(i){
if(sum(c("sumTfBn","sumTnBn") %in% names(fixef(i)))==2){
# Get variance of for each fixed effects by multiplying coefficients by design matrix
var.vec <- var(as.vector(fixef(i)[c("sumTfBn","sumTnBn")] %*% t(i@pp$X[,c("sumTfBn","sumTnBn")])))
# Get variance of fixed effects by multiplying coefficients by design matrix
VarF <- var(as.vector(fixef(i) %*% t(i@pp$X)))
# Get variance of random effects by extracting variance components
VarRand <- colSums(do.call(rbind,lapply(VarCorr(i),function(j) j[1])))
# Get residual variance
VarResid <- attr(VarCorr(i),"sc")^2
var.vec <- var.vec/(VarF+VarRand+VarResid)
}else{
var.vec <- NA
}    
names(var.vec) <- paste("effect.response","VAR",sep=".")
return(var.vec)
}





coef(.)$ecocode.i
ranef(.,whichel='ecocode.id',condVar=TRUE)