SCM

SCM Repository

[matrix] Diff of /pkg/src/lmeRep.c
ViewVC logotype

Diff of /pkg/src/lmeRep.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 409, Sun Dec 19 02:19:58 2004 UTC revision 410, Tue Dec 28 14:51:42 2004 UTC
# Line 1651  Line 1651 
1651      UNPROTECT(1);      UNPROTECT(1);
1652      return Omg;      return Omg;
1653  }  }
1654    
1655    /**
1656     * Create the crosstabulation.  Should do this early in case we need to
1657     * reorder the factors.
1658     *
1659     * @param flist pointer to the factor list.
1660     * @param nobs number of observations.
1661     * @param nc number of columns in the model matrices.
1662     *
1663     * @return the pairwise crosstabulation in the form of the ZtZ array.
1664     */
1665    static SEXP
1666    lmer_crosstab(SEXP flist, int nobs, const int nc[])
1667    {
1668        int i, nf = length(flist);
1669        int npairs = (nf * (nf + 1))/2;
1670        SEXP val = PROTECT(allocVector(VECSXP, npairs));
1671        SEXP cscbCl = MAKE_CLASS("cscBlocked");
1672        int *Ti = Calloc(nobs, int),
1673            *nlevs = Calloc(nf, int),
1674            **zb = Calloc(nf, int*); /* zero-based indices */
1675        double *ones = Calloc(nobs, double),
1676            *Tx = Calloc(nobs, double);
1677    
1678        for (i = 0; i < nobs; i++) ones[i] = 1.;
1679        for (i = 0; i < nf; i++) {  /* populate the zb vectors */
1680            SEXP fi = VECTOR_ELT(flist, i);
1681            int j;
1682    
1683            zb[i] = Calloc(nobs, int);
1684            nlevs[i] = length(getAttrib(fi, R_LevelsSymbol));
1685            for (j = 0; j < nobs; j++) zb[i][j] = INTEGER(fi)[j] - 1;
1686            for (j = 0; j <= i; j++) {
1687                int *ijp, ind = Lind(i, j), nnz;
1688                SEXP ZZij;
1689    
1690                SET_VECTOR_ELT(val, ind, NEW_OBJECT(cscbCl));
1691                ZZij = VECTOR_ELT(val, ind);
1692                SET_SLOT(ZZij, Matrix_pSym, allocVector(INTSXP, nlevs[j] + 1));
1693                ijp = INTEGER(GET_SLOT(ZZij, Matrix_pSym));
1694                triplet_to_col(nlevs[i], nlevs[j], nobs, zb[i], zb[j], ones,
1695                               ijp, Ti, Tx);
1696                nnz = ijp[nlevs[j]];
1697                SET_SLOT(ZZij, Matrix_iSym, allocVector(INTSXP, nnz));
1698                Memcpy(INTEGER(GET_SLOT(ZZij, Matrix_iSym)), Ti, nnz);
1699                SET_SLOT(ZZij, Matrix_xSym, alloc3Darray(REALSXP, nc[i], nc[j], nnz));
1700                /* The crosstab counts are copied into the first nnz elements */
1701                /* of the x slot.  These aren't the correct array positions */
1702                /* unless nc[i] == nc[j] == 1 but we don't use them. */
1703                Memcpy(REAL(GET_SLOT(ZZij, Matrix_xSym)), Tx, nnz);
1704            }
1705        }
1706    
1707        for (i = 0; i < nf; i++) Free(zb[i]);
1708        Free(zb); Free(nlevs); Free(ones); Free(Ti); Free(Tx);
1709        UNPROTECT(1);
1710        return val;
1711    }
1712    
1713    
1714    /**
1715     * Update the arrays ZtZ, ZtX, and XtX in an lme object
1716     * according to a list of model matrices.
1717     *
1718     * @param x pointer to an lmer object
1719     * @param mmats pointer to a list of model matrices
1720     *
1721     * @return NULL
1722     */
1723    SEXP lmer_update_mm(SEXP x, SEXP mmats)
1724    {
1725        SEXP
1726            ZtZP = GET_SLOT(x, Matrix_ZtZSym),
1727            ZtXP = GET_SLOT(x, Matrix_ZtXSym),
1728            flist = GET_SLOT(x, Matrix_flistSym);
1729        int *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1730            *dims = INTEGER(getAttrib(ZtXP, R_DimSymbol)),
1731            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1732            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
1733            nf = length(flist), nfp1 = nf + 1,
1734            i, ione = 1,
1735            nobs = nc[nfp1],
1736            pp1 = nc[nf];
1737        double
1738            *X,
1739            *XtX = REAL(GET_SLOT(x, Matrix_XtXSym)),
1740            *ZtX = REAL(ZtXP),
1741            one = 1.0, zero = 0.0;
1742    
1743        if (!isNewList(mmats) || length(mmats) != nfp1)
1744            error("mmats must be a list of %d model matrices", nfp1);
1745        for (i = 0; i <= nf; i++) {
1746            SEXP mmat = VECTOR_ELT(mmats, i);
1747            int *mdims = INTEGER(getAttrib(mmat, R_DimSymbol));
1748    
1749            if (!isMatrix(mmat) || !isReal(mmat))
1750                error("element %d of mmats is not a numeric matrix", i + 1);
1751            if (nobs != mdims[0])
1752                error("Expected %d rows in the %d'th model matrix. Got %d",
1753                      nobs, i+1, mdims[0]);
1754            if (nc[i] != mdims[1])
1755                error("Expected %d columns in the %d'th model matrix. Got %d",
1756                      nc[i], i+1, mdims[1]);
1757        }
1758                                    /* Create XtX */
1759        X = REAL(VECTOR_ELT(mmats, nf));
1760        F77_CALL(dsyrk)("U", "T", &pp1, &nobs, &one, X, &nobs, &zero, XtX, nc + nf);
1761                                    /* Zero an accumulator */
1762        memset((void *) ZtX, 0, sizeof(double) * pp1 * Gp[nf]);
1763        for (i = 0; i < nf; i++) {
1764            int *fac = INTEGER(VECTOR_ELT(flist, i)),
1765                j, k, nci = nc[i], ZtXrows = Gp[i+1] - Gp[i];
1766            int ncisqr = nci * nci, nlev = ZtXrows/nci;
1767            double *Z = REAL(VECTOR_ELT(mmats, i)), *ZZx;
1768    
1769            for (k = 0; k < i; k++) {
1770                SEXP ZZxM = VECTOR_ELT(ZtZP, Lind(i, k));
1771                int *rowind = INTEGER(GET_SLOT(ZZxM, Matrix_iSym)),
1772                    *colptr = INTEGER(GET_SLOT(ZZxM, Matrix_pSym));
1773                int *f2 = INTEGER(VECTOR_ELT(flist, k)), nck = nc[k];
1774                double *Zk = REAL(VECTOR_ELT(mmats, k));
1775    
1776                ZZx = REAL(GET_SLOT(ZZxM, Matrix_xSym));
1777                memset(ZZx, 0, sizeof(double) *
1778                       length(GET_SLOT(ZZxM, Matrix_xSym)));
1779                for (j = 0; j < nobs; j++) {
1780                    F77_CALL(dgemm)("T", "N", nc + i, nc + k, &ione, &one,
1781                                    Z + j, &nobs, Zk + j, &nobs, &one,
1782                                    ZZx + Tind(rowind, colptr, fac[j] - 1, f2[j] - 1)
1783                                    * (nci * nck), &nci);
1784                }
1785            }
1786            ZZx = REAL(GET_SLOT(VECTOR_ELT(ZtZP, Lind(i, i)), Matrix_xSym));
1787            memset((void *) ZZx, 0, sizeof(double) * nci * nci * nlev);
1788            if (nci == 1) {         /* single column in Z */
1789                for (j = 0; j < nobs; j++) {
1790                    int fj = fac[j] - 1; /* factor indices are 1-based */
1791                    ZZx[fj] += Z[j] * Z[j];
1792                    F77_CALL(daxpy)(&pp1, Z + j, X + j, &nobs, ZtX + fj, dims);
1793                }
1794            } else {
1795                for (j = 0; j < nobs; j++) {
1796                    int fj = fac[j] - 1; /* factor indices are 1-based */
1797    
1798                    F77_CALL(dsyr)("U", nc + i, &one, Z + j, &nobs,
1799                                   ZZx + fj * ncisqr, nc + i);
1800                    F77_CALL(dgemm)("T", "N", nc + i, &pp1, &ione,
1801                                    &one, Z + j, &nobs,
1802                                    X + j, &nobs, &one,
1803                                    ZtX + fj * nci, dims);
1804                }
1805            }
1806            ZtX += ZtXrows;
1807        }
1808        status[0] = status[1] = 0;
1809        return R_NilValue;
1810    }
1811    
1812    /**
1813     * Create an lmer object from a list grouping factors and a list of model
1814     * matrices.  There is one more model matrix than grouping factor.  The last
1815     * model matrix is the fixed effects and the response.
1816     *
1817     * @param facs pointer to a list of grouping factors
1818     * @param ncv pointer to a list of model matrices
1819     *
1820     * @return pointer to an lmer object
1821     */
1822    SEXP lmer_create(SEXP flist, SEXP mmats)
1823    {
1824        SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("lmer")));
1825        SEXP ZtZ, cnames, fnms, nms;
1826        int *nc, i, nf = length(flist), nobs;
1827    
1828                                    /* Check validity of flist */
1829        if (!(nf > 0 && isNewList(flist)))
1830            error("flist must be a non-empty list");
1831        nobs = length(VECTOR_ELT(flist, 0));
1832        if (nobs < 1) error("flist[[0]] must be a non-null factor");
1833        for (i = 0; i < nf; i++) {
1834            SEXP fi = VECTOR_ELT(flist, i);
1835            if (!(isFactor(fi) && length(fi) == nobs))
1836                error("flist[[%d]] must be a factor of length %d",
1837                      i + 1, nobs);
1838        }
1839        SET_SLOT(val, Matrix_flistSym, flist);
1840        if (!(isNewList(mmats) && length(mmats) == (nf + 1)))
1841            error("mmats must be a list of length %d", nf + 1);
1842        SET_SLOT(val, Matrix_ncSym, allocVector(INTSXP, nf + 2));
1843        nc = INTEGER(GET_SLOT(val, Matrix_ncSym));
1844        nc[nf + 1] = nobs;
1845        for (i = 0; i <= nf; i++) {
1846            SEXP mi = VECTOR_ELT(mmats, i);
1847            int *dims;
1848    
1849            if (!(isMatrix(mi) && isReal(mi)))
1850                error("mmats[[%d]] must be a numeric matrix", i + 1);
1851            dims = INTEGER(getAttrib(mi, R_DimSymbol));
1852            if (dims[0] != nobs)
1853                error("mmats[[%d]] must have %d rows", i + 1, nobs);
1854            if (dims[1] < 1)
1855                error("mmats[[%d]] must have at least 1 column", i + 1);
1856            nc[i] = dims[1];
1857        }   /* Arguments have now been checked for type, dimension, etc. */
1858                                    /* Create pairwise crosstabulation in ZtZ */
1859        SET_SLOT(val, Matrix_ZtZSym, lmer_crosstab(flist, nobs, nc));
1860        lmer_populate(val);
1861        ZtZ = GET_SLOT(val, Matrix_ZtZSym);
1862        /* FIXME: Check for possible reordering of the factors to maximize the
1863         * number of levels in the initial sequence of nested factors. */
1864        fnms = getAttrib(flist, R_NamesSymbol);
1865                                    /* Allocate and populate nc and cnames */
1866        SET_SLOT(val, Matrix_cnamesSym, allocVector(VECSXP, nf + 1));
1867        cnames = GET_SLOT(val, Matrix_cnamesSym);
1868        setAttrib(cnames, R_NamesSymbol, allocVector(STRSXP, nf + 1));
1869        nms = getAttrib(cnames, R_NamesSymbol);
1870        for (i = 0; i <= nf; i++) {
1871            SEXP mi = VECTOR_ELT(mmats, i);
1872            SET_VECTOR_ELT(cnames, i,
1873                           duplicate(VECTOR_ELT(getAttrib(mi, R_DimNamesSymbol),
1874                                                1)));
1875            if (i < nf)
1876                SET_STRING_ELT(nms, i, duplicate(STRING_ELT(fnms, i)));
1877            else
1878                SET_STRING_ELT(nms, nf, mkChar(".fixed"));
1879        }
1880        lmer_update_mm(val, mmats);
1881        UNPROTECT(1);
1882        return val;
1883    }
1884    
1885    /**
1886     * Create and insert initial values for Omega.
1887     *
1888     * @param x pointer to an lmer object
1889     *
1890     * @return NULL
1891     */
1892    SEXP lmer_initial(SEXP x)
1893    {
1894        SEXP Omg = GET_SLOT(x, Matrix_OmegaSym);
1895        int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)), i, nf = length(Omg);
1896    
1897        for (i = 0; i < nf; i++) {
1898            SEXP ZZxP = GET_SLOT(VECTOR_ELT(GET_SLOT(x, Matrix_ZtZSym), Lind(i, i)),
1899                                 Matrix_xSym);
1900            int *dims = INTEGER(getAttrib(ZZxP, R_DimSymbol));
1901            int j, k, nzc = dims[0], nlev = dims[2];
1902            int nzcsqr = nzc * nzc, nzcp1 = nzc + 1;
1903            double *Omega = REAL(VECTOR_ELT(Omg, i)),
1904                mi = 0.375 / ((double) nlev);
1905    
1906            memset((void *) Omega, 0, sizeof(double) * nzc * nzc);
1907            for (j = 0; j < nlev; j ++) {
1908                for (k = 0; k < nzc; k++) {
1909                    Omega[k * nzcp1] += REAL(ZZxP)[k * nzcp1 + j * nzcsqr] * mi;
1910                }
1911            }
1912        }
1913        status[0] = status[1] = 0;
1914        return R_NilValue;
1915    }
1916    
1917    /**
1918     * Copy the diagonal blocks from ZtZ to ZZpO and inflate by Omega.
1919     * Update devComp[1].
1920     *
1921     * @param x pointer to an lmer object
1922     */
1923    SEXP
1924    lmer_inflate(SEXP x)
1925    {
1926        SEXP Omg = GET_SLOT(x, Matrix_OmegaSym),
1927            ZZpO = GET_SLOT(x, Matrix_ZZpOSym),
1928            ZtZ = GET_SLOT(x, Matrix_ZtZSym);
1929        int *dcmp = INTEGER(GET_SLOT(x, Matrix_devCompSym)),
1930            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1931            i, nf = length(Omg);
1932    
1933        dcmp[0] = dcmp[1] = dcmp[2] = dcmp[3] = 0.;
1934        for (i = 0; i < nf; i++) {
1935            int ind = Lind(i, i), j, nci = nc[i], ncisqr = nci * nci;
1936            SEXP ZZOel = VECTOR_ELT(ZZpO, ind);
1937            SEXP ZZm = GET_SLOT(VECTOR_ELT(ZtZ, ind), Matrix_xSym),
1938                ZZOm = GET_SLOT(ZZOel, Matrix_xSym);
1939            int nlev = INTEGER(getAttrib(ZZm, R_DimSymbol))[2],
1940                *ZZOd = INTEGER(getAttrib(ZZOm, R_DimSymbol));
1941            int *colptr = INTEGER(GET_SLOT(ZZOel, Matrix_pSym)),
1942                *rowind = INTEGER(GET_SLOT(ZZOel, Matrix_iSym));
1943            double *ZZ = REAL(ZZm), *ZZO = REAL(ZZOm),
1944                *Omega = REAL(VECTOR_ELT(Omg, i));
1945    
1946            /* zero the whole of the ZZpO component */
1947            memset(ZZO, 0, sizeof(double) * ZZOd[0] * ZZOd[1] * ZZOd[2]);
1948            if (nci == 1) {
1949                dcmp[1] += nlev * log(Omega[0]);
1950                for (j = 0; j < nlev; j++) { /* inflate */
1951                    ZZO[Tind(rowind, colptr, j, j)] = ZZ[j] + Omega[0];
1952                }
1953            } else {
1954                double *tmp = Memcpy(Calloc(ncisqr, double), Omega, ncisqr);
1955    
1956                F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
1957                if (j)
1958                    error("Leading minor of size %d of Omega[[%d]] is not positive definite",
1959                          j, i + 1);
1960                for (j = 0; j < nci; j++) { /* nlev * logDet(Omega_i) */
1961                    dcmp[1] += nlev * 2. * log(tmp[j * (nci + 1)]);
1962                }
1963                Free(tmp);
1964    
1965                for (j = 0; j < nlev; j++) {
1966                    int ii, jj;
1967                    double
1968                        *ZZOj = ZZO + Tind(rowind, colptr, j, j) * ncisqr,
1969                        *ZZj = ZZ + j * ncisqr;
1970    
1971                    for (jj = 0; jj < nci; jj++) { /* Copy ZZi to ZZOj and inflate */
1972                        for (ii = jj; ii < nci; ii++) { /* upper triangle only */
1973                            int ind = ii + jj * nci;
1974    
1975                            ZZOj[ind] = ZZj[ind] + Omega[ind];
1976                        }
1977                    }
1978                }
1979            }
1980        }
1981        return R_NilValue;
1982    }
1983    
1984    /**
1985     * If status[["factored"]] is FALSE, create and factor Z'Z+Omega, then
1986     * create RZX and RXX, the deviance components, and the value of the
1987     * deviance for both ML and REML.
1988     *
1989     * @param x pointer to an lmeRep object
1990     *
1991     * @return NULL
1992     */
1993    SEXP lmer_factor(SEXP x)
1994    {
1995        int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
1996    
1997        if (!status[0]) {
1998            SEXP DP = GET_SLOT(x, Matrix_DSym),
1999                LP = GET_SLOT(x, Matrix_LSym),
2000                Linv = GET_SLOT(x, Matrix_LinvSym),
2001                RZXsl = GET_SLOT(x, Matrix_RZXSym),
2002                ZZOP = GET_SLOT(x, Matrix_ZZpOSym);
2003            int *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
2004                *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
2005                *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
2006                i, j, nf = length(DP);
2007            int nml = nc[nf + 1], nreml = nml + 1 - nc[nf];
2008            double
2009                *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
2010                *RZX = REAL(RZXsl),
2011                *dcmp = REAL(GET_SLOT(x, Matrix_devCompSym)),
2012                *deviance = REAL(GET_SLOT(x, Matrix_devianceSym)),
2013                minus1 = -1., one = 1.;
2014    
2015    
2016            lmer_inflate(x);
2017            for (i = 0; i < nf; i++) {
2018                SEXP ZZOiP = VECTOR_ELT(ZZOP, Lind(i, i));
2019                SEXP DiP = VECTOR_ELT(DP, i);
2020                int nlev = INTEGER(getAttrib(DiP, R_DimSymbol))[2],
2021                    *colptr = INTEGER(GET_SLOT(ZZOiP, Matrix_pSym)),
2022                    *rowind = INTEGER(GET_SLOT(ZZOiP, Matrix_iSym));
2023                int jj, nci = nc[i], ncisqr = nci * nci;
2024                double *D = REAL(DiP), *ZZOi = REAL(ZZOiP);
2025    
2026                if (nci == 1) {
2027                    for (j = 0; j < nlev; j++) {
2028                        double ZZO = ZZOi[Tind(rowind, colptr, j, j)];
2029                        dcmp[0] += log(ZZO);
2030                        D[j] = sqrt(ZZO);
2031                    }
2032                } else {
2033                    for (j = 0; j < nlev; j++) {
2034                        double *Dj = Memcpy(D + j * ncisqr,
2035                                            ZZOi + Tind(rowind, colptr, j, j) * ncisqr,
2036                                            ncisqr);
2037    
2038                        F77_CALL(dpotrf)("U", &nci, Dj, &nci, &jj);
2039                        if (jj)
2040                            error("D[ , , %d], level %d, is not positive definite",
2041                                  j + 1, i + 1);
2042                        for (jj = 0; jj < nci; jj++) /* accumulate determinant */
2043                            dcmp[0] += 2. * log(Dj[jj * (nci + 1)]);
2044                    }
2045                }
2046            }
2047                                    /* solve for RZX */
2048            Memcpy(RZX, REAL(GET_SLOT(x, Matrix_ZtXSym)), dims[0] * dims[1]);
2049            L_inverse(nf, LP, Linv);
2050            lmeRep_L_sm(nf, LP, Linv, RZXsl);
2051            for (i = 0; i < nf; i++) {
2052                SEXP DiP = VECTOR_ELT(DP, i);
2053                int nci = nc[i], ncisqr = nci * nci,
2054                    nlev = INTEGER(getAttrib(DiP, R_DimSymbol))[2];
2055    
2056                for (j = 0; j < nlev; j++) {
2057                    F77_CALL(dtrsm)("L", "U", "T", "N", &nci, dims + 1,
2058                                    &one, REAL(DiP) + j * ncisqr, &nci,
2059                                    RZX + Gp[i] + j * nci, dims);
2060                }
2061            }
2062                                    /* downdate and factor XtX */
2063            Memcpy(RXX, REAL(GET_SLOT(x, Matrix_XtXSym)), dims[1] * dims[1]);
2064            F77_CALL(dsyrk)("U", "T", dims + 1, dims,
2065                            &minus1, REAL(RZXsl), dims,
2066                            &one, RXX, dims + 1);
2067            F77_CALL(dpotrf)("U", dims + 1, RXX, dims + 1, &j);
2068            if (j) {
2069                warning("Leading minor of size %d of downdated X'X is indefinite",
2070                        j);
2071                dcmp[2] = dcmp[3] = deviance[0] = deviance[1] = NA_REAL;
2072            } else {
2073                for (j = 0; j < (dims[1] - 1); j++) /* 2 logDet(RXX) */
2074                    dcmp[2] += 2 * log(RXX[j * (dims[1] + 1)]);
2075                dcmp[3] = 2. * log(RXX[dims[1] * dims[1] - 1]); /* 2 log(ryy) */
2076                deviance[0] =       /* ML criterion */
2077                    dcmp[0] - dcmp[1] + nml*(1.+dcmp[3]+log(2.*PI/nml));
2078                deviance[1] = dcmp[0] - dcmp[1] + /* REML */
2079                    dcmp[2] + nreml*(1.+dcmp[3]+log(2.*PI/nreml));
2080            }
2081            status[0] = 1; status[1] = 0; /* factored but not inverted */
2082        }
2083        return R_NilValue;
2084    }

Legend:
Removed from v.409  
changed lines
  Added in v.410

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge