<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
  <head>
    <meta content="text/html; charset=UTF-8" http-equiv="Content-Type">
    <title></title>
  </head>
  <body bgcolor="#ffffff" text="#000000">
    On 12/06/2010 10:00 PM, Peter Tittmann wrote:
    <blockquote cite="mid:80C1F8072EB74C4F903BC96C09FB1D6A@gmail.com"
      type="cite">
      <p> </p>
      Gentlemen,
      <div><br>
      </div>
      <div><span id="goog_749721654"></span>I've decided to switch to
        the OLS method, thought I did get the NNLS method that Skipper
        proposed working. I was not prepared to spend more time trying
        to make sense of the resulting array for ancova, etc. (also
        could not figure out how to interpret the resulting coefficient
        array as I was expecting a 2d array representing the a and b
        coefficient values but it returned a 1d array). I have hopefully
        simple follow up questions:<span id="goog_749721664"></span><span
          id="goog_749721665"></span></div>
      <div><span id="goog_749721666"></span><span id="goog_749721667"></span><br>
      </div>
      <div>1. Is there a method to define explicitly the function used
        in OLS? I know numpy.linalg.lstsq is the way OLS works but is
        there another function where I can define the form? <span
          id="goog_749721672"></span><span id="goog_749721673"></span><span
          id="goog_749721669"></span></div>
      <div><span id="goog_749721674"></span><span id="goog_749721675"></span><br>
      </div>
      <div>2. I'm still interested in interpreting the results of the
        NNLS method, so if either of you can suggest what the resulting
        arrays mean id be grateful. I've attached the output of NNLS </div>
      <div><span id="goog_749721670"></span><span id="goog_749721671"></span><br>
      </div>
      <div>warm regards,<span id="goog_749721678"></span><span
          id="goog_749721679"></span></div>
      <div><span id="goog_749721680"></span><span id="goog_749721681"></span><br>
      </div>
      <div>Peter<span id="goog_749721682"></span><span
          id="goog_749721683"></span></div>
      <div><span id="goog_749721684"></span><span id="goog_749721685"></span><br>
      </div>
      <div><br>
      </div>
      <div><span id="goog_749721662"></span>Here is the working version
        of NNLS:<span id="goog_749721676"></span><span
          id="goog_749721677"></span></div>
      <div><br>
      </div>
      <div>
        <div>def getDiam2(ht,*b):</div>
        <div>    return b[0] * ht[:,1]**b[1] + np.sum(b[2:]*ht[:,2:],
          axis=1)</div>
        <div><br>
        </div>
        <div>dt = np.genfromtxt('/home/peter/Desktop/db_out.csv',
          delimiter=",", names=True, dtype=None)</div>
        <div> </div>
        <div>indHtPlot = adt['height']</div>
        <div>depDbh = adt['dbh']</div>
        <div>plot_dummies, col_map = sm.tools.categorical(dt['plot],
          drop=True, dictnames=True)</div>
        <div><br>
        </div>
        <div><br>
        </div>
        <div>def nnlsDummies():</div>
        <div>    '''this function returns coefficients and covariance
          arrays'''</div>
        <div>    plot_dummies, col_map = sm.tools.categorical(indPlot,
          drop=True, dictnames=True)</div>
        <div>    X = np.column_stack((indHt, plot_dummies))</div>
        <div>    Y = depDbh</div>
        <div>    coefs, cov = curve_fit(getDiam2, X, Y, p0=
          [0.]*X.shape[1])</div>
        <div>    return coefs, cov</div>
      </div>
      <div><span id="goog_749721658"></span><span id="goog_749721659"></span><span
          id="goog_749721660"></span><span id="goog_749721661"></span><br>
        <br>
        <div id="BBCDC7FCEA99477283C223B44F299898">-- <br>
          Peter Tittmann
          <div><br>
          </div>
        </div>
        <br>
        <p>On Monday, December 6, 2010 at 4:55 PM, <a class="moz-txt-link-abbreviated" href="mailto:josef.pktd@gmail.com">josef.pktd@gmail.com</a>
          wrote:</p>
        <blockquote type="cite" style="border-left-style: solid;
          border-width: 1px; margin-left: 0px; padding-left: 10px;">
          <div>
            <div>
              <div>On Mon, Dec 6, 2010 at 7:41 PM, Skipper Seabold &lt;<a
                  moz-do-not-send="true"
                  href="mailto:jsseabold@gmail.com">jsseabold@gmail.com</a>&gt;
                wrote:<br>
                <blockquote type="cite">
                  <div> On Mon, Dec 6, 2010 at 7:31 PM, Peter Tittmann <ptittmann@gmail.com>
                      wrote:<br>
                      &gt; thanks both of you,<br>
                      &gt; Josef, the data that I sent is only the first
                      100 rows of about 1500, there<br>
                      &gt; should be sufficient sampling in each plot.<br>
                      &gt; Skipper, I have attempted to deploy your
                      suggestion for not linearizing the<br>
                      &gt; data. It seems to work. I'm a little confused
                      at your modification if the<br>
                      &gt; getDiam function and I wonder if you could
                      help me understand. The form of<br>
                      &gt; the equation that is being fit is:<br>
                      &gt; Y= a*X^b<br>
                      &gt; your version of the detDaim function:<br>
                      &gt;<br>
                      &gt; def getDiam(ht, *b):<br>
                      &gt;    return ht[:,0]**b[0] +
                      np.sum(b[1:]*ht[:,1:], axis=1)<br>
                      &gt;<br>
                      &gt; Im sorry if this is an obvious question but I
                      don't understand how this<br>
                      &gt; works as it seems that the "a" coefficient is
                      missing.<br>
                      &gt; Thanks again!<br>
                      <br>
                      Right.  I took out the 'a', because as I read it
                      when I linearized (I<br>
                      might be misunderstanding ancova, I never recall
                      the details), if you<br>
                      include 'a' and also all of the dummy variables
                      for the plot, then you<br>
                      will have a the problem of multicollinearity.  You
                      could also include<br>
                      'a' and drop one of the plot dummies, but then 'a'
                      is just your<br>
                      reference category that you dropped.  So now b[0]
                      is the nonlinear<br>
                      effect of your main variable and b[1:] contains
                      linear shift effects<br>
                      of all the plots.  Hmm, thinking about it some
                      more, though I think<br>
                      you could include 'a' in the non-linear version
                      above (call it b[0]<br>
                      and shift everything else over by one), because
                      now 'a' would be the<br>
                      effect when the current b[0] is zero.  I was just
                      unsure how you meant<br>
                      'a' when you had a*ht**b and were trying to
                      include in ht the plot<br>
                      variable dummies.<br>
                    </ptittmann@gmail.com></div>
                </blockquote>
                <br>
                As I understand it, the intention is to estimate
                equality of the slope<br>
                coefficients, so the continuous variable is multiplied
                with the dummy<br>
                variables. In this case, the constant should still be
                added. The<br>
                normalization question is whether to include all
                dummy-cont.variable<br>
                products and drop the continuous variable, or include
                the continuous<br>
                variable and drop one of the dummy-cont levels.<br>
                <br>
                Unless there is a strong reason to avoid log-normality
                of errors, I<br>
                would work (first) with the linear version.<br>
                <br>
                Josef<br>
                <br>
                <br>
                <blockquote type="cite">
                  <div><br>
                    Skipper<br>
                    _______________________________________________<br>
                    SciPy-User mailing list<br>
                    <a moz-do-not-send="true"
                      href="mailto:SciPy-User@scipy.org">SciPy-User@scipy.org</a><br>
                    <a moz-do-not-send="true"
                      href="http://mail.scipy.org/mailman/listinfo/scipy-user">http://mail.scipy.org/mailman/listinfo/scipy-user</a><br>
                    <br>
                  </div>
                </blockquote>
                _______________________________________________<br>
                SciPy-User mailing list<br>
                <a moz-do-not-send="true"
                  href="mailto:SciPy-User@scipy.org">SciPy-User@scipy.org</a><br>
                <a class="moz-txt-link-freetext" href="http://mail.scipy.org/mailman/listinfo/scipy-user">http://mail.scipy.org/mailman/listinfo/scipy-user</a><br>
              </div>
            </div>
          </div>
        </blockquote>
        <br>
      </div>
      <pre wrap="">
<fieldset class="mimeAttachmentHeader"></fieldset>
_______________________________________________
SciPy-User mailing list
<a class="moz-txt-link-abbreviated" href="mailto:SciPy-User@scipy.org">SciPy-User@scipy.org</a>
<a class="moz-txt-link-freetext" href="http://mail.scipy.org/mailman/listinfo/scipy-user">http://mail.scipy.org/mailman/listinfo/scipy-user</a>
</pre>
    </blockquote>
    I do think this is starting to be an off-list discussion because
    this is really about statistics and not about numpy/scipy (you can
    contact me off-list if you want).<br>
    <br>
    I am not sure what all the variables are so please excuse me but I
    presume you want to model dbh as a function of height, plot and
    species. Following usual biostatistics interpretation, 'plot' is
    probably treated as random effect but you probably have to use R/SAS
    etc for that for both linear and nonlinear models or some spatial
    models. <br>
    <br>
    Really you need to determine whether or not a nonlinear model is
    required. With the few data points you provided, I only see a linear
    relationship between dbh and height with some outliers and perhaps
    some heterogeneity of variance. Often doing a simple
    polynomial/spline can help to see if there is any evidence for a
    nonlinear relationship in the full data - a linear model or
    polynomial with the data provided does not suggest a nonlinear
    model. Obviously a linear model is easier to fit and interpret
    especially if you create the design matrix as estimable functions
    (which is rather trivial once you understand using dummy variables).<br>
    <br>
    The most general nonlinear/multilevel model proposed is of the form:<br>
    dbh= C + A*height^B<br>
    Obviously if B=1 then it is a linear model and the parameters A, B
    and C can be modeled with a linear function of intercept, plot and
    species. Although, if 'plot' is what I think it is then you probably
    would not model the parameters A and B with it. <br>
    <br>
    Without C you are forcing the curve through zero which is biological
    feasible if you expect dbh=0 when height is zero. However, dbh can
    be zero if height is not zero just due to the model itself or what
    dbh actually is (it may take a minimum height before dbh is greater
    than zero). With the data you provided, there are noticeable
    differences between species for dbh and height so you probably need
    C in your model.<br>
    <br>
    For this general model you probably should just fit the curve for
    each species alone but I would use a general stats package to do
    this. This will give you a good starting point to know how well the
    curve fits each species as well as the similarity of parameters and
    residual variation. Getting convergence with a model that has B
    varying across species may be rather hard so I would suggest
    modeling A and C first.  <br>
    <br>
    Bruce<br>
    <br>
    <br>
    <br>
    <br>
    <br>
    <br>
    <br>
  </body>
</html>