title: What is ga() in the gamlss package doing?
body: <p>I have been looking into the gamlss package for fitting semiparametric models and came across something strange in the ga() function. Even if the model is specified as having a gamma distribution, fitted using REML, the output for the model is Gaussian, fitted using GCV.</p>
  <h2>Example::</h2>
  <pre><code>library(mgcv)
  library(gamlss)
  library(gamlss.add)
  data(rent)
  ga3 &lt;- gam(R~s(Fl)+s(A), method=&quot;REML&quot;, data=rent, family=Gamma(log))
  gn3 &lt;- gamlss(R~ga(~s(Fl)+s(A), method=&quot;REML&quot;), data=rent, family=GA)
  </code></pre>
  <h2>Model summary for the GAM::</h2>
  <pre><code>summary(ga3)
  </code></pre>
  <blockquote>
  <pre><code>Family: Gamma 
  Link function: log 
  
  Formula:
  R ~ s(Fl) + s(A)
  
  Parametric coefficients:
              Estimate Std. Error t value Pr(&gt;|t|)    
  (Intercept) 6.667996   0.008646   771.2   &lt;2e-16 ***
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  Approximate significance of smooth terms:
          edf Ref.df      F p-value    
  s(Fl) 1.263  1.482 442.53  &lt;2e-16 ***
  s(A)  4.051  4.814  36.34  &lt;2e-16 ***
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  R-sq.(adj) =  0.302   Deviance explained = 28.8%
  -REML =  13979  Scale est. = 0.1472    n = 1969
  </code></pre>
  </blockquote>
  <h2>Model summary for the GAMLSS::</h2>
  <pre><code>summary(getSmo(gn3))
  </code></pre>
  <blockquote>
  <pre><code>Family: gaussian 
  Link function: identity 
  
  Formula:
  Y.var ~ s(Fl) + s(A)
  
  Parametric coefficients:
               Estimate Std. Error t value Pr(&gt;|t|)
  (Intercept) 6.306e-13  8.646e-03       0        1
  
  Approximate significance of smooth terms:
          edf Ref.df      F p-value    
  s(Fl) 1.269  1.492 440.14  &lt;2e-16 ***
  s(A)  3.747  4.469  38.83  &lt;2e-16 ***
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  R-sq.(adj) =  0.294   Deviance explained = 29.6%
  GCV = 0.97441  Scale est. = 0.97144   n = 1969
  </code></pre>
  </blockquote>
  <h2>Question::</h2>
  <p>Why is the model output giving the incorrect distribution and fitting method? Is there something that I am missing here and this is correct?</p>