diff --git a/NEWS.md b/NEWS.md
index 0159fae3..27e9f120 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,8 @@
+Version 2.0.
+=============
+
+* Added VA implementation of Tweedie
+
 Version 2.0
 =============
 
diff --git a/R/gllvm.R b/R/gllvm.R
index 007a4411..bd841450 100644
--- a/R/gllvm.R
+++ b/R/gllvm.R
@@ -1250,10 +1250,10 @@ gllvm <- function(y = NULL, X = NULL, TR = NULL, data = NULL, formula = NULL, fa
     #   cat("VA method cannot handle", family, " family, so LA method is used instead. \n")
     #   method <- "LA"
     # }
-    if (quadratic != FALSE && family %in% c("tweedie", "beta")){
+    if (quadratic != FALSE && family %in% c("beta")){
       stop("The quadratic model is not implemented for ", family, " family yet. \n")
     }
-    if (method == "VA" && family %in% c("tweedie", "beta")){
+    if (method == "VA" && family %in% c("beta")){
       cat("Note that, the", family, "family is implemented using the extended variational approximation method. \n")
     }
     # if (method == "EVA"){
diff --git a/src/gllvm.cpp b/src/gllvm.cpp
index 5bedfe1a..db572547 100644
--- a/src/gllvm.cpp
+++ b/src/gllvm.cpp
@@ -2704,7 +2704,7 @@ Type objective_function<Type>::operator() ()
           if(!gllvmutils::isNA(y(i,j)))nll -= ( -eta(i,j) - exp(-eta(i,j)+cQ(i,j))*y(i,j) )*iphi(j) + log(y(i,j)*iphi(j))*iphi(j) - log(y(i,j)) -lgamma(iphi(j));
         }
       }
-    } else if(family==5){ // Tweedie EVA
+    } else if((family==5) && (method >1)){ // Tweedie EVA
       //Type ePower = extra(0);
       ePower = invlogit(ePower) + Type(1);
       for (int i=0; i<n; i++) {
@@ -2721,7 +2721,45 @@ Type objective_function<Type>::operator() ()
           }
         }
       }
-    } else if(family==6){ 
+    }else if((family==5) && (method <1)){ // Tweedie VA
+      ePower = invlogit(ePower) + Type(1);
+      for (int i=0; i<n; i++) {
+        for (int j=0; j<p; j++) {
+          if(!gllvmutils::isNA(y(i,j))){
+            if(y(i,j) == 0){
+              nll += pow(exp(eta(i,j)+(2-ePower)*cQ(i,j)), 2-ePower) / (iphi(j)*(2-ePower)); // log(prob(y=0))
+            } else if(y(i,j)>0){
+              CppAD::vector<Type> tx(4);
+              tx[0] = y(i,j);
+              tx[1] = iphi(j);
+              tx[2] = ePower;
+              tx[3] = 0;
+              nll -= atomic::tweedie_logW(tx)[0];
+              nll -= (1/iphi(j)) * (y(i,j)*pow(exp(eta(i,j)+(1-ePower)*cQ(i,j)),1-ePower)/(1-ePower)-pow(exp(eta(i,j)+(2-ePower)*cQ(i,j)),2-ePower)/(2-ePower));
+              nll += log(y(i,j));
+            }
+            // Type p1 = ePower - 1.0, p2 = 2.0 - ePower;
+            // nll += pow(exp(eta(i,j)+p2*cQ(i,j)), p2) / (iphi(j) * p2); // log(prob(y=0))
+            // if (y(i,j) > 0) {
+            //   CppAD::vector<Type> tx(4);
+            //   tx[0] = y(i,j);
+            //   tx[1] = iphi(j);
+            //   tx[2] = ePower;
+            //   tx[3] = 0;
+            //   nll -= atomic::tweedie_logW(tx)[0];
+            //   nll += y(i,j) / (iphi(j) * p1 * pow(exp(eta(i,j)+p1*cQ(i,j)), p1)) + log(y(i,j));
+            //   // if (CppAD::Variable(y)) {
+            //   //   ans += CppAD::CondExpGt(y, Type(0), ans2, Type(0));
+            //   // } else {
+            //   //   ans += ans2;
+            //   // }
+            // }
+              
+            // }
+          }
+        }
+      }
+    }else if(family==6){ 
       iphi = iphi/(1+iphi);
       Type pVA;
       for (int j=0; j<p;j++){