using different starting points, different results

R/likelihood.R

@@ -120,7 +120,7 @@ wright <- function(data, tree, regimes, alpha=1, sigma=1, Xo=NULL, ...){
 
 
 # OUCH
-ouch <- function(data, tree, regimes, alpha=1, sigma=1, Xo=NULL){
+ouch <- function(data, tree, regimes, alpha=1, sigma=1, Xo=NULL, ...){
 # alpha is fixed at ~zero, sigma is regime dependent, theta is global
 
     # intialize a parameter vector to optimize: 

demos/loop_models_traits_regimes.R

@@ -24,7 +24,7 @@ fit_all <- function(models, traits, regimes, tree){
             out[[i]] <- brown(traits[j], tree)
         } else if (models[[i]]=="ou1"){
           out[[i]] <- hansen(traits[j], tree,  regime=labrid$noregimes,
-                             .01, .01, control=list(maxit=5000))
+                             .01, .01, maxit=5000)
 
 
         ## Simple hansen model
@@ -39,7 +39,7 @@ fit_all <- function(models, traits, regimes, tree){
 #                      "regime=", names(regimes[k])))
           ## use hansen to start with good parameters 
           hansen <- try(fit("hansen", traits[j], regimes[[k]], tree=tree,
-                            control=list(maxit=5000)))
+                            maxit=5000))
           ## hansen will sometimes give very large or negative 
 
           ## Fit one of the generalized models using the initial guess from hansen
@@ -92,7 +92,71 @@ llik_matrix <- function(results){
 }
 
 
+
+sort_lik <- function(likmat, model_names = c("bm", "ou", "theta",
+                      "sigma", "gen", "alpha")){
+## Takes output of llik_matrix and:
+## Reorganize the data to list by ascending score
+## Rename models (with short, parameter-based names) 
+## rescale the data relative to weakest model 
+lliks <- vector("list", length=length(likmat))
+for(i in 1:length(likmat)){
+  tmp <- likmat[[i]]
+  names(tmp) <- model_names
+  tmp <- sort(tmp)
+  shift <- tmp[1] 
+  for(j in 1:length(tmp)){
+    print(tmp[j])
+    tmp[j] <- tmp[j]-shift
+    print(tmp[j])
+  }
+  lliks[[i]] <- tmp
+  print(tmp)
+}
+names(lliks) <- names(likmat)
+lliks
+}
+
+
+alpha_traits <- function(results, model_name="release_constraint"){
+# Takes the output of fit_all and grabs alpha values for specified model
+  k <- 1 # the regime index
+  groups <- levels(results$regimes[[k]])
+  n_groups <- length(groups)
+  i <- match(model_name, results$models)
+  M <- matrix(NA, ncol=n_groups, nrow=length(results$traits))
+  for(j in 1:length(results$traits)){
+    a <- results$fit[[j]][[k]][[i]]
+    M[j,] <- a$alpha
+  }
+  rownames(M) <- names(results$traits)
+  colnames(M) <- groups
+  M
+}
+
+
+conv<- function(results){
+## Takes the output of fit_all and reports which ones didn't converge
+  for(j in 1:length(results$traits)){
+    for(k in 1:length(results$regimes)){
+      for(i in 1:length(results$models)){
+        a <- results$fit[[j]][[k]][[i]]
+        if (is(a, "multiOU"))
+          if (a$convergence != 0)
+            print(paste("trait ", names(results$traits)[j], 
+            ", regime ", names(results$regimes)[k], ", model ", 
+            results$models[i], " didn't converge", sep=""))
+      }
+    }
+  }
+}
+
+
+
+
+
 alpha_matrix <- function(results, trait_id, k = 1){
+## Takes the output of fit_all and computs the matrix of alphas 
   k <- 1
   j <- trait_id 
   models <- results$models
@@ -125,36 +189,5 @@ alpha_matrix <- function(results, trait_id, k = 1){
   M
 }
 
-alpha_traits <- function(results, release=TRUE){
-# 
-  k <- 1 # the regime index
-  groups <- levels(results$regimes[[k]])
-  n_groups <- length(groups)
-  if (release){
-    i <- match("release_constraint", results$models)
-  } else {
-    i <- match("wright", results$models)
-  }
-  M <- matrix(NA, ncol=n_groups, nrow=length(results$traits))
-  for(j in 1:length(results$traits)){
-    a <- results$fit[[j]][[k]][[i]]
-    M[j,] <- a$alpha
-  }
-  rownames(M) <- names(results$traits)
-  colnames(M) <- groups
-  M
-}
 
 
-conv<- function(results){
-for(j in 1:length(results$traits)){
-for(k in 1:length(results$regimes)){
-for(i in 1:length(results$models)){
-  a <- results$fit[[j]][[k]][[i]]
-  if (is(a, "multiOU"))
-    if (a$convergence != 0)
-      print(paste("trait ", names(results$traits)[j], 
-      ", regime ", names(results$regimes)[k], ", model ", 
-      results$models[i], " didn't converge", sep=""))
-}}}
-}

demos/parrotfish.R

@@ -11,30 +11,19 @@ model_list <- list("brown", "hansen", "ouch", "brownie", "wright", "release_cons
 regime_list <-  list(intramandibular=intramandibular)
 
 test <- fit_all(model_list, labrid$data, regime_list, labrid$tree)
+conv(test)
 
-likmat <- llik_matrix(test)
 
 
-## Reorganize the data to list by ascending score
-## Rename models (with short, parameter-based names) 
-## rescale the data relative to weakest model 
-lliks <- vector("list", length=length(likmat))
-for(i in 1:length(likmat)){
-  tmp <- likmat[[i]]
-  names(tmp) <-c("bm", "ou", "theta", "sigma", "gen", "alpha")
-  tmp <- sort(tmp)
-  shift <- tmp[1] 
-  for(j in 1:length(tmp)){
-    print(tmp[j])
-    tmp[j] <- tmp[j]-shift
-    print(tmp[j])
-  }
-  lliks[[i]] <- tmp
-  print(tmp)
-}
-names(lliks) <- names(likmat)
+## Summary matrices
+likmat <- llik_matrix(test)
+lliks <- sort_lik(likmat)
+release_alphas <- alpha_traits(test) 
+wright_alphas <- alpha_traits(test, model_name="wright") 
 
 
+#### Plots
+
 # Body mass, ratios 
 png("parrotfish_ratios.png", width=480*1.5, height=480*1.5)
 par(mfrow=c(2,2))
@@ -63,27 +52,27 @@ for(i in c(1,2,10,11)){
   barplot(lliks[[i]], horiz=T, main=names(lliks)[i])
 }
 dev.off()
-require(socialR)
-flickr(files="parrotfish*.png", tag=tag)
 
 
-release_alphas <- alpha_traits(test) 
-
 png("release_alphas.png", width=480*2)
  barplot(log(release_alphas[,1]/release_alphas[,2]))
 dev.off()
 
 
-wright_alphas <- alpha_traits(test, release=F) 
-
 png("wright_alphas.png", width=480*2)
  barplot(log(wright_alphas[,1]/wright_alphas[,2]))
 dev.off()
 
-flickr(files="*alphas.png", tag=tag)
+#flickr(files="parrotfish*.png", tag=tag)
+#flickr(files="*alphas.png", tag=tag)
+
+
+
+
+
+
 
 
-conv(test)
 ## sometimes hansen returns a silly large value of alpha, which will break the generalized likelihood function (too stiff)
  #h <- fit("hansen", labrid$data[13], intramandibular, labrid$tree, .01, .01)
  # w <- fit("release_constraint", labrid$data[13], intramandibular, labrid$tree, h@sqrt.alpha^2, h@sigma)