title: Leave-one-out cross validation

text: I am working about elaborate a simple linear regression so I need to evaluate many models. I am asking how I can use leave-one-out cross validation to validate a simple linear regression.


title: Spatially explicit model to test effects of multiple variables

text: What is the appropriate model to test the effects of location, species, and size on tree growth? So I have both categorical and continuous variables I want in the model. The following graphic just shows a single plot measured in two consecutive time periods, each point scaled by its size and colored by species. I want to model the change tree (point) size as a function of neighbor size, distance, species, and initial size. ![enter image description here](http://i.stack.imgur.com/h0vhj.png) Thanks!


title: How to interpret Cochran-Mantel-Haenszel test?

text: I'm testing the independence of two variables, A and B, stratified by C. A and B are binary variables and C is categorical (5 values). Running Fisher's exact test for A and B (all strata combined), I get:               ##          (B)     ##      (A) FALSE TRUE     ##    FALSE  1841   85     ##    TRUE    915   74          OR: 1.75 (1.25 --  2.44), p = 0.0007 *      where OR is the odds ratio (estimate and 95% confidence interval), and `*` means that p < 0.05. Running the same test for each stratum (C), I get:               C=1, OR: 2.31 (0.78 --  6.13), p = 0.0815     C=2, OR: 2.75 (1.21 --  6.15), p = 0.0088 *     C=3, OR: 0.94 (0.50 --  1.74), p = 0.8839     C=4, OR: 1.48 (0.77 --  2.89), p = 0.2196     C=5, OR: 3.38 (0.62 -- 34.11), p = 0.1731      Finally, running Cochran-Mantel-Haenszel (CMH) test, using A, B, and C, I get:               OR: 1.56 (1.12 --  2.18), p = 0.0089 *      The result from the CMH test suggests that A and B are **not** independent at each stratum (p < 0.05); however, most of the within stratum tests were non- significant, which would suggest that we don't have enough evidence to discard that A and B are independent at each stratum. So, what conclusion is right? How to report the conclusion given those results? Can C be considered a confounding variable or not? EDIT: I performed the Breslow-Day test for the null hypothesis that the odds ratio is the same across strata, and the p-value was 0.1424.


title: What is usefulness of time series decomposition if building a multivariate model?

text: I am at the early stages of building a panel regression model of sales data. I know my final model dataset will consist of log sales, control variables and log media variables.   I am planning to use stl() to decompose the log sales time series. My question: If I am going to create a multivariate time series model anyway, does it make sense to decompose the dependent variable during the Exploratory Data Analysis (EDA) phase? Does the seasonal/ trend component in the decomp tell me something about my timeseries different than I would expect controlling for these factors on the RHS of the model equation?


title: How to do and interpret a generalized estimating equation

text: My study is looking at the effects of enclosure type (1 IV (penned v not penned)) on the stereotypical behaviours and interactions (2 DV's (counts)) of elephants, however using rain (yes/no) and high and low temperatures as covariates. The same 9 elephants were studied last year while in pens and again this year when not in pens. I have been told I need to do a GEE but I have no idea how to go about doing this, let alone understand the output. Nor do I know how this should be set out in SPSS.


title: Suitable method for modelling (underdispersed?) count data with lots of zeros and long tail

text: I have a small data set of counts of bees. I tried a simple Poisson model without random effects but it was very overdispersed (3.95). When I fit a GLMM with random effects (using glmer in lme4) it was then underdispersed (0.19). I need to include the random effects because of my set-up. The data is distributed like this, with lots of zeroes and two large values: ![histogram of count data](http://i.stack.imgur.com/ecklZm.jpg) In addition, most of the zeros occur in one level of the main predictor variable, like this:                   with(SB, table(solbee, treatment) )      ![tabulated counts](http://i.stack.imgur.com/qZ3Go.png) I tried a negative binomial mixed model (glmer.nb in lme4) which was a better fit but still not right, also ZIP and ZINB (using package glmmADMB) but none are a good fit. I tried a hurdle model in glmmADMB but got the following error message: Model:               hurdle1 <- glmmadmb(solbee ~    treatment + landuse + snh +                              (1|site/dayfac),                             data=subset(SB,solbee>0),                             family="truncnbinom1")      Error message:               Error in model.frame.default(formula = solbee ~ treatment + landuse +  :      variable lengths differ (found for 'treatment')      I think this is because most of the zero counts are in one level of the treatment factor. Any advice as to what to try next would be very welcome, as would any info on why my data is overdispersed _without_ random effects, and underdispersed _with_ random effects.


title: Word Sense Disambiguation in Practice

text: I have a question that might seem very obvious but I don't really have a good answer for it. There are many algorithms out there that deal with word sense disambiguation but all of the ones that I have seen assume that the ambiguous words are known apriori. i.e. usually they are applied to a dataset where there are known ambiguous words (eg. plant, jet, bank, etc). My question is, what if you are given some completely new text dataset that might have ambiguous words in it or might not. Obviously before you can apply a word sense disambiguation algorithm to them you must know which words to apply it to. So my question is, given an unknown piece of text, how do we know which words are ambiguous (words we should apply a WSD algorithm to) and which are not. I can think of two possible methods 1) compile a large list of ambiguous words and check if any of them are in the text 2) do it in a more machine learning oriented like judging how ambiguous a word is based on its collocations somehow.


title: Suppose X and Y are each continuous. Are they jointly continuous?

text: My question is in the title. I do not think that X and Y are jointly continuous but I cannot think of a counterexample where x and y would not be jointly continuous.


title: Proper similarity measure for clustering

text: I have problems in finding a proper similarity measure for clustering. I have around 3000 arrays of sets, where each set contains features of certain domain (e.g., number, color, days, alphabets, etc). I'll explain my problem with an example. Lets assume i have only 2 arrays(a1 & a2) and I want to find the similarity between them. each array contains 4 sets (in my actual problem there are 250 sets (domains) per array) and a set can be empty.               a1: {a,b}, {1,4,6}, {mon, tue, wed}, {red, blue,green}     a2: {b,c}, {2,4,6}, {}, {blue, black}      I have come with a similarity measure using Jaccard index (denoted as J):               sim(a1,a2) = [J(a1[0], a2[0]) + J(a1[1], a2[1]) + ... + J(a1[3], a2[3])]/4      note:I divide by total number of sets (in the above example 4) to keep the similarity between 0 and 1. `Is this a proper similarity measure and are there any flaws in this approach`. I am applying Jaccard index for each set separately because I want compare the similarity between related domains(i.e. color with color, etc...) I am not aware of any other proper similarity measure for my problem. Further, `can I use this similarity measure for clustering purpose?`


title: Are PCA solutions unique?

text: When I run PCA on a certain data set, is the solution given to me unique? I.e., I obtain a set of 2d coordinates, based on interpoint distances. Is it possible to find at least one more arrangement of the points that would meet these constraints? If the answer is yes, how can I find such different solution?


title: Which logit/probit model do I use for multiple reponse/dependent variables?

text: I have $300$ time series objects that constitute the $300$ columns of matrix $X$. This matrix has $5$ rows and represents $5$ days of time series information for each $300$ columns. I set up a $300\textrm{x}5$ matrix of binary values. So the first row might look something like $(1 1 0 0 1)$, which would mean that column $1$ of $X$ had negative elements in the 1st, 2nd and 5th rows (modeled by category "$1$"), and positive elements in the 3rd and 4th rows (modeled by category "$0$"). I have $8$ predictor variables that explains/relates well to whether there are more negative elements early (i.e. if we observe something like $(1 1 0 0 0)$) versus whether there are more negative elements later (i.e. we observe something like $(0 0 0 1 1)$). How do I model this? My knowledge is not expansive enough to know if there's any statistical framework where I can use this matrix as my response: $ Y = \left( \begin{array}{ccccc} 1 & 1 & 0 & 0 & 1 \\\ 1 & 0 & 1 & 1 & 0 \\\ 0 & 0 & 0 & 1 & 1 \\\ . & . & . & . & . \\\ \end{array} \right) $ Using this I want to do something like: $\mathbf{Y} = a + \mathbf{B}\mathbf{Z}$ (I understand that this is not a correct logit/probit specification, but I think it's fine for communicating what I want). Then the ideal interpretation of some coefficient $b_i \in \mathbf{B}$ would be; if it's negative, then we have increased odds of seeing something like $(11000)$ instead of something like $(00011)$ if the associated predictor $x_i \in \mathbf{Z}$ is larger.


title: How robust are multivariate methods to violations of normality?

text: In many cases, multivariate methods are used without normality tests. How are following methods robust if data are not normal?   * Principal components analysis   * Canonical correlations analysis   * Factor analysis   * Multivariate regression analysis   * Multivariate analysis of variance   * Discriminant analysis


title: Similarity measure between two colors?

text: If one feature of my objects is color, how to compute similarity between two objects? It is correct to transform it to RGB and use Euclidean distance in 3D? Or is there any acceptable one dimensional ordering of colors (as in the case of rainbow)?


title: Large scale ridge regression

text: I'm trying to solve a problem of the form $\min_x \frac{1}{2}||Ax-b||^2_2 + \frac{\rho}{2}||x-z||^2_F$ where both $x$ and $b$ are high dimensional, and $b$ is much higher dimensional than $x$. The solution is given by $x^* = (A^T A+\rho I)^{-1}(A^T b + z)$, but the problem is so large that even inverting $A^T A + \rho I$ is infeasible. However, due to structure in the problem we can efficiently multiply by $A$ and $A^T$. Basically this is large scale linear ridge regression. What would be the ideal algorithm for efficiently implementing this minimization? Would something like biconjugate gradient work?


title: predicting manufacturing survival time

text: My task is to improve widget quality from a high volume manufacturing process. The survival data is 99% right censored, since most products do not fail, and rapid analysis long before failure is necessary to get the business value. The goals are to:   1. change the manufacturing process to stop making high risk widgets. This needs "effect size" in the language of Survival Analysis   2. select widgets at high risk of return for separate handling. This needs widget level predictions. The right censoring makes Survival Analysis appropriate (which I haven't done before) relative to ordinary prediction (which I've done lots of). I have thousands of covariates about the widgets, making closed form methods desirable. Accelerated Failure Time models seem popular in manufacturing. The log-logistic distribution is a commonly used Accelerated Failure Time (AFT) model with a closed form. My questions are:   1. Can I get effect sizes from AFT?   2. Can I get observation level predictions from AFT?   3. Can I scale AFT to tens of thousands of covariates?   4. Can you propose a Python implementation of an AFT?   5. Can you propose an R implementation of an AFT?   6. Do you have an alternate to AFT to propose? thanks


title: Partitioning error in mixed factorial anova in R

text: I need help understanding a 3-way anova with 2 fixed factors and 1 random factor. I am analyzing an experiment looking at the effect of diet on feeding rate in snails. The two fixed factors are diet (2 levels) and herbivore (2 levels snails present, snails absent). The experiment was repeated three times, and in each trial the fixed factors were fully crossed and replicated 3 times. I would like to analyze "trial" (3 levels) as a random factor. Here's a my data: > nochoiceb3                  Trial Diet Container Herbivore  lmg.day     1      1   22         1    snails 6.513823     2      1   22         2    snails 6.282641     3      1   22         3    snails 6.087683     4      1   22         4   control 5.132263     5      1   22         5   control 3.407842     6      1   22         6   control 4.987708     7      1   10         7    snails 3.632309     8      1   10         8    snails 4.301359     9      1   10         9    snails 4.364372     10     1   10        10   control 4.295924     11     1   10        11   control 2.433613     12     1   10        12   control 3.015535     13     2   22         2    snails 4.646312     14     2   22        10    snails 5.370638     15     2   22        17    snails 4.888844     16     2   22         1   control 5.129899     17     2   22         6   control 4.284965     18     2   22        13   control 4.237001     19     2   10         3    snails 4.239887     20     2   10        11    snails 4.133565     21     2   10        19    snails 4.130355     22     2   10         4   control 2.415914     23     2   10         7   control 2.219203     24     2   10        14   control 2.844909     25     3   22         2    snails 5.906451     26     3   22        11    snails 5.558371     27     3   22        18    snails 5.015954     28     3   22         5   control 4.284965     29     3   22        12   control 4.678421     30     3   22        13   control 4.371976     31     3   10         3    snails 3.775057     32     3   10        15    snails 4.028917     33     3   10        23    snails 4.433195     34     3   10         1   control 1.974081     35     3   10         4   control 2.797281     36     3   10        16   control 3.427515      As I understand, including a random factor changes the partitioning of the error in an anova, and the F-ratio should then be calculated based on the the correct error term. In this case the F-ratio for the fixed factors should be the ratio of the mean squares of the fixed factor to the mean squares of the fixed factor:random factor interaction. Ex. for diet F= Mean squares diet/Mean squares diet:trial I ran the following model: > model2.1<-aov(lmg.day~Diet*Herbivore+Error(Trial), data=nochoiceb3) And get these results:               Error: Trial               Df Sum Sq Mean Sq F value Pr(>F)     Residuals  2  1.543  0.7717                         Error: Within                    Df Sum Sq Mean Sq F value   Pr(>F)         Diet            1 22.283  22.283  75.268 1.12e-09 ***     Herbivore       1 12.686  12.686  42.853 3.06e-07 ***     Diet:Herbivore  1  0.096   0.096   0.324    0.573         Residuals      30  8.881   0.296       I wanted to check the results by calculating the F-ratios myself, so I ran a 3-way anova with trial as a fixed factor to generate sums of squares for the interactions:               model2.2<-aov(lmg.day~Diet*Herbivore*Trial, data=nochoiceb3)                      Df Sum Sq Mean Sq F value   Pr(>F)         Diet                  1 22.283  22.283  78.350 5.03e-09 ***     Herbivore             1 12.686  12.686  44.607 6.58e-07 ***     Trial                 2  1.543   0.772   2.713   0.0866 .       Diet:Herbivore        1  0.096   0.096   0.338   0.5666         Diet:Trial            2  0.135   0.067   0.237   0.7909         Herbivore:Trial       2  0.111   0.056   0.196   0.8235         Diet:Herbivore:Trial  2  1.810   0.905   3.181   0.0595 .       Residuals            24  6.826   0.284       The f-ratio given by the mixed model is not the same as you get by calculating the f-ratio based on dividing the mean squares by the correct error term (ex. for diet F= Mean squares diet/Mean squares diet:trial) but instead using the residuals as the error term. I don't know where to go from here. So: Why does this difference happen? Is what r is doing correct, or does this mean that my mixed model is wrong? How can I ensure that I am using the right syntax for my model in R? Thanks for any help!


title: Calculating random effects from a glmerMod object (r package lme4)

text: Using the `lme4` package, how does `ranef()` calculate (or extract) estimates of random effects from a `glmerMod` object? The random effects seem to be related to `model@u`, but what is the relationship? On a related note, does anyone know of any sources that explain in detail what each of the components of a `glmerMod` object actually are? I can figure out some of them (e.g. `model@beta`), but most of them remain a mystery.


title: Fisher Scoring v/s Coordinate Descent for MLE in R

text: R base function glm() uses Fishers Scoring for MLE, while the glmnet uses the coordinate descent method to solve the same equation ? Coordinate descent is more time efficient than Fisher Scoring as fisher scoring calculates the second order derivative matrix and some other matrix operation which makes it space and time expensive, while coordinate descent can do the same task in O(np) time. Why R base function uses Fisher Scoring or this method has advantage over other optimization methods? What will be comparison between coordinate descent and Fisher Scoring ? I am relatively new to do this field so any help or resource will be helpful


title: Tail probabilities and the GHK simulator

text: I am trying to use the GHK simulator to estimate the probabilities $F(\mathbf{x} > k\mathbf{a})$ that the values of a high dimensional ($n>1000$), correlated random vector $\mathbf{x}$ will exceed some threshold vector $\mathbf{a}$, linearly scaled by $k$ with $k \rightarrow \infty$. The GHK approach seems simple enough (e.g. here), however the results of my MATLAB implementation are pretty disappointing. For my scenario I know that for large values of $k$, $F$ will be fairly small (e.g. $<10^{-6}$). (I can verify this using a brute force Monte-Carlo approach.) However, after some good initial estimates for small $k$, the probabilities returned by GHK simulator begin fluctuating around some value (e.g. $ \sim 10^{-3}$) and eventually drop to zero (as the calculated CDFs become ones). My questions are:   * Is the GHK simulator fundamentally unsuitable for use with high dimensional data, or is there some workaround which avoids numerical problems?    * How many sample vectors should I be generating? Currently I use 1000 samples for each value of $k$. Using a larger number improves the results slightly, however then, in terms of computation time, the advantage of the GHK approach compared to a brute-force MC approach is negated.    * Is there a more suitable approach for this type of problem?


title: Find cumulative distribution function

text: ![enter image description here](http://i.stack.imgur.com/NNEU4.png) It is the second part that I am having troubles with. I've tried Y from x to x+u and X from 0 to 1: didn't work ! So what are the correct limits in this case(that give you the required answer) ?


title: Measuring the "distance" between two multivariate distributions

text: I'm looking for some good terminology to describe what I'm trying to do, to make it easier to look for resources. So, say I have two clusters of points A and B, each associated to two values, X and Y, and I want to measure the "distance" between A and B - i.e. how likely is it that they were sampled from the same distribution (I can assume that the distributions are normal). For example, if X and Y are correlated in A but not in B, the distributions are different. Intuitively, I would get the covariance matrix of A, and then look at how likely each point in B is to fit in there, and vice-versa (probably using someting like Mahalanobis distance). But that is a bit "ad-hoc", and there is probably a more rigorous way of describing this (of course, in practice I have more than two datasets with more than two variables - I'm trying to identify which of my datasets are outliers). Thanks!


title: Optimal lag length in VECM using vars R package

text: I have some series that are cointegrated, so I know that I should fit a VECM model. Nevertheless I found no guidance in finding the optimal lag length, say lagLength. I am using vars R package. In order to check cointegration I used ca.jo(..,K=cointegrationLength) then I used cajorls to fit the VECM cajoorls(...,K=lagLength). I do not understand: 1\. the right interaction between K in both functions. 2\. A "optimal" criteria to choose the lag lengths


title: How do I determine how well a discrete variable is predicted by a continuous response variable?

text: I have a response variable which is the number of days it takes for an insect to emerge from the pupal stage. I measured this daily, and most of the individuals emerged after about 20-21 days, with a tail of some individuals trailing as far off as 33 days. For ease of analysis, I have a variable called, "ClaustPer", which is this number of days minus 20. ClaustPer is not continuous (it can only have a narrow range of integer values), but it also has quite a number of "categories", if they can be called that. I want to figure out whether MLH, which is the average heterozygosity of these individual insects (across a number of genes), is a significant predictor of the number of days it take for an insect to emerge from the pupal stage. I'm currently using generalized linear models (GLMs), but I'm not sure whether that's appropriate and/or how to designate the link function, so I thought I would make this more a question of how best to figure out which type of analysis is best suited to a not-quite-categorical response variable. I provide the data set below in R (note that there are missing data designated as NA and some individuals have been omitted prior to analysis):               data<-structure(list(ID = c(1, 2, 3, 4, 7, 9, 10, 12, 13, 14, 15, 16,      17, 18, 20, 21, 22, 23, 24, 25, 27, 28, 29, 31, 34, 37, 38, 39,      40, 41, 43, 44, 45, 46, 47, 48, 49, 52, 55, 56, 59, 60, 61, 62,      63, 65, 66, 67, 68, 69, 71, 72, 73, 74, 75, 76, 77, 79, 80, 81,      82, 83, 84, 85, 86, 87, 88, 89, 90, 92, 93, 94, 95, 98, 99, 100,      102, 103, 104, 105, 107, 108, 109, 111, 113, 114, 115, 116, 117,      118, 119, 121, 123, 124, 125, 126, 127, 128, 129, 130, 131, 133,      135, 136, 137, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148,      149, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162,      163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175,      177, 178, 179, 182, 183, 184, 185, 186, 187, 188, 189, 191, 192,      194, 195, 197, 198, 199, 202, 204, 205, 206, 208, 209, 210, 211,      212, 213, 214, 215, 216, 217, 218, 219, 222, 224, 226, 228, 229,      232, 233, 234, 235, 236, 237, 238, 239, 241, 242, 244, 245, 247,      248, 251, 254, 255, 259, 260, 261, 262, 266, 267, 268, 269, 270,      271, 273, 275, 276, 278, 279, 280, 281, 282, 285, 286, 287, 289,      291, 292, 295, 296, 297, 298, 299, 302, 303, 305, 306, 307, 309,      310, 312, 313, 314, 315, 316, 318, 319, 320, 321, 322, 324, 325,      327, 328, 329, 330, 333, 334, 335, 336, 337, 338, 339, 340, 341,      343, 344, 345, 346, 347, 349, 350, 351, 352, 353, 354, 355, 356,      357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 369, 371,      372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384,      385, 387, 388, 389, 391, 392, 393, 394, 395, 396, 397, 401, 402,      403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415,      416, 417, 418, 419, 420, 421, 422, 423, 426, 427, 430, 434, 437,      438, 439, 440, 442, 443, 445, 447, 448, 449, 450, 451, 452, 453,      454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 467, 468,      469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481,      482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494,      495, 496, 497, 499, 500, 501, 502, 503, 505, 507, 508, 509, 510,      511, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 522, 523,      524, 525, 526, 527, 528, 529, 530, 531, 532, 533, 534, 535, 536,      537, 538, 539, 540, 541, 542, 543, 544, 545, 546, 547, 548, 549,      550, 551, 552, 555, 556, 557, 558, 559, 560, 561, 562, 563, 564,      565, 566, 567, 568, 569, 570, 571, 572, 573, 574, 575, 576, 577,      578, 579, 580, 581, 582, 583, 584, 585, 586, 587, 588, 589, 590,      591, 592, 593, 594, 595, 596, 597, 598, 599, 600, 601, 603, 604,      605, 606, 607, 608, 609, 610, 611, 612, 613, 614, 615, 616, 617,      618, 619, 620, 621, 622, 623, 624, 625, 626, 627, 628, 630, 631,      632, 633, 634, 635, 636, 637, 638, 639, 641, 642, 643, 644, 645,      646, 647, 648, 649, 650, 651, 652), MLH = c(0.534090909090909,      0.5, NA, 0.40506329113924, 0.298507462686567, 0.410958904109589,      0.293103448275862, 0.442105263157895, 0.48, 0.554347826086957,      0.453488372093023, 0.535353535353535, 0.443298969072165, 0.304878048780488,      0.457446808510638, 0.455555555555556, 0.397849462365591, 0.494252873563218,      0.48314606741573, 0.377777777777778, 0.457446808510638, 0.445652173913043,      0.3, 0.412371134020619, 0.354838709677419, 0.464646464646465,      0.474226804123711, 0.43010752688172, 0.46078431372549, 0.541666666666667,      0.494623655913978, 0.446808510638298, 0.413793103448276, 0.462962962962963,      0.426966292134831, 0.494736842105263, 0.522222222222222, 0.460674157303371,      0.524390243902439, 0.506493506493506, 0.505747126436782, 0.395348837209302,      0.419354838709677, 0.456521739130435, 0.469387755102041, 0.375,      0.4125, 0.295454545454545, 0.436363636363636, 0.486486486486487,      0.411111111111111, 0.338028169014085, 0.319587628865979, 0.368421052631579,      0.422222222222222, 0.381443298969072, 0.479591836734694, 0.483870967741935,      0.469387755102041, 0.397849462365591, 0.427083333333333, 0.340659340659341,      0.44, 0.26027397260274, 0.408602150537634, 0.463917525773196,      0.36734693877551, 0.5, 0.402061855670103, 0.454545454545455,      0.382022471910112, 0.494736842105263, 0.458333333333333, 0.434343434343434,      0.505154639175258, 0.397959183673469, 0.494949494949495, 0.484848484848485,      0.460674157303371, 0.389473684210526, 0.571428571428571, 0.417582417582418,      0.533980582524272, 0.538461538461538, 0.5, 0.52, 0.421686746987952,      0.451923076923077, 0.51, 0.456310679611651, 0.422680412371134,      0.427184466019417, 0.526315789473684, 0.480392156862745, 0.376344086021505,      0.463414634146341, 0.455445544554455, 0.505154639175258, 0.520408163265306,      0.442307692307692, 0.431372549019608, 0.46, 0.39, 0.464646464646465,      0.427184466019417, 0.354430379746835, 0.505050505050505, 0.49,      0.451923076923077, 0.416666666666667, 0.52, 0.389473684210526,      0.355769230769231, 0.47, 0.407766990291262, 0.443298969072165,      0.297872340425532, 0.386138613861386, 0.43298969072165, 0.44,      0.391304347826087, 0.436170212765957, 0.46078431372549, 0.446808510638298,      0.461538461538462, 0.42156862745098, 0.423076923076923, 0.529411764705882,      0.445544554455446, 0.436893203883495, 0.504854368932039, 0.436893203883495,      0.563106796116505, 0.38, 0.436170212765957, 0.336538461538462,      0.490384615384615, 0.442307692307692, 0.46078431372549, 0.461538461538462,      0.346153846153846, 0.447916666666667, 0.427083333333333, 0.343434343434343,      0.387755102040816, 0.307692307692308, 0.442105263157895, 0.436170212765957,      0.445544554455446, 0.480769230769231, 0.495145631067961, 0.42156862745098,      0.402061855670103, 0.408163265306122, 0.382352941176471, 0.45,      0.524752475247525, 0.423076923076923, 0.444444444444444, 0.456521739130435,      0.422680412371134, 0.524271844660194, 0.42156862745098, 0.480769230769231,      0.574257425742574, 0.445544554455446, 0.490384615384615, 0.407766990291262,      0.509615384615385, 0.451923076923077, 0.344086021505376, 0.5,      0.519230769230769, 0.519607843137255, 0.470588235294118, 0.407766990291262,      0.466019417475728, 0.310679611650485, 0.231578947368421, 0.423076923076923,      0.368932038834951, 0.489795918367347, 0.5, 0.455445544554455,      0.475728155339806, 0.470588235294118, 0.378640776699029, 0.423529411764706,      0.485148514851485, 0.529411764705882, 0.465346534653465, 0.474226804123711,      0.47, 0.535353535353535, 0.446601941747573, 0.388349514563107,      0.444444444444444, 0.509803921568627, 0.405940594059406, 0.459770114942529,      0.42156862745098, 0.4375, 0.509803921568627, 0.438775510204082,      0.450980392156863, 0.484848484848485, 0.391752577319588, 0.452631578947368,      0.468085106382979, 0.406976744186047, 0.427184466019417, 0.538461538461538,      0.461538461538462, 0.425742574257426, 0.505376344086022, 0.495049504950495,      0.38, 0.488636363636364, 0.434782608695652, 0.531914893617021,      0.419354838709677, 0.515463917525773, 0.410526315789474, 0.446808510638298,      0.494623655913978, 0.430232558139535, 0.438775510204082, 0.394230769230769,      0.403846153846154, 0.415841584158416, 0.455445544554455, 0.455445544554455,      0.484848484848485, 0.379310344827586, 0.393939393939394, 0.504854368932039,      0.378640776699029, 0.441176470588235, 0.408163265306122, 0.439024390243902,      0.495145631067961, 0.432692307692308, 0.378640776699029, 0.468085106382979,      0.5, 0.473684210526316, 0.479591836734694, 0.480392156862745,      0.494949494949495, 0.365853658536585, 0.475728155339806, 0.333333333333333,      0.333333333333333, 0.514563106796116, 0.464646464646465, 0.461538461538462,      0.442307692307692, 0.494252873563218, 0.456310679611651, 0.54,      0.441176470588235, 0.495145631067961, 0.355555555555556, 0.461538461538462,      0.43, 0.407766990291262, 0.431372549019608, 0.427184466019417,      0.466019417475728, 0.519607843137255, 0.366336633663366, 0.403846153846154,      0.407766990291262, 0.405940594059406, 0.489583333333333, 0.567307692307692,      0.398058252427184, 0.423076923076923, 0.480769230769231, 0.425742574257426,      0.466019417475728, 0.528846153846154, 0.524271844660194, 0.414634146341463,      0.384615384615385, 0.469879518072289, 0.465116279069767, 0.517647058823529,      0.387096774193548, 0.448275862068966, 0.443181818181818, 0.41304347826087,      0.449438202247191, 0.402439024390244, 0.463157894736842, 0.395061728395062,      0.376344086021505, 0.479591836734694, 0.519607843137255, 0.529411764705882,      0.538461538461538, 0.446601941747573, 0.398058252427184, 0.470588235294118,      0.424242424242424, 0.461538461538462, 0.398058252427184, 0.520408163265306,      0.408163265306122, 0.445544554455446, 0.474226804123711, 0.414141414141414,      0.5, 0.505050505050505, 0.422680412371134, 0.425287356321839,      0.404040404040404, 0.428571428571429, 0.514563106796116, 0.388349514563107,      0.411764705882353, 0.396039603960396, 0.432692307692308, 0.383838383838384,      0.5, 0.46875, 0.435643564356436, 0.478723404255319, 0.480769230769231,      0.432692307692308, 0.509803921568627, 0.485436893203884, 0.461538461538462,      0.451612903225806, 0.5, 0.47, 0.455445544554455, 0.397959183673469,      0.49, 0.395833333333333, 0.4375, 0.473684210526316, 0.436170212765957,      0.407766990291262, 0.440860215053763, 0.39, 0.384615384615385,      0.574257425742574, 0.44, 0.397849462365591, 0.416666666666667,      0.376470588235294, 0.471910112359551, 0.478723404255319, 0.406593406593407,      0.425531914893617, 0.41, 0.484848484848485, 0.469387755102041,      0.366666666666667, 0.457446808510638, 0.465346534653465, 0.408163265306122,      0.454545454545455, 0.398058252427184, 0.466019417475728, 0.470588235294118,      0.444444444444444, 0.469387755102041, 0.373737373737374, 0.402061855670103,      0.45, 0.435643564356436, 0.431578947368421, 0.377551020408163,      0.415841584158416, 0.411764705882353, 0.450549450549451, 0.441176470588235,      0.450980392156863, 0.363636363636364, 0.414141414141414, 0.401960784313726,      0.392156862745098, 0.490196078431373, 0.41747572815534, 0.435643564356436,      0.372549019607843, 0.475728155339806, 0.46, 0.450980392156863,      0.411764705882353, 0.509803921568627, 0.309278350515464, 0.564356435643564,      0.474226804123711, 0.504854368932039, 0.466019417475728, 0.441176470588235,      0.455445544554455, 0.419354838709677, 0.401960784313726, 0.466019417475728,      0.46078431372549, 0.510204081632653, 0.494949494949495, 0.490196078431373,      0.438775510204082, 0.431372549019608, 0.445544554455446, 0.484848484848485,      0.524271844660194, 0.427184466019417, 0.441176470588235, 0.474747474747475,      0.42156862745098, 0.46, 0.411764705882353, 0.392156862745098,      0.495145631067961, 0.49, 0.407766990291262, 0.520408163265306,      0.42, 0.445544554455446, 0.442105263157895, 0.463917525773196,      0.386138613861386, 0.47, 0.485148514851485, 0.5, 0.529411764705882,      0.333333333333333, 0.4, 0.514563106796116, 0.45360824742268,      0.46, 0.490196078431373, 0.553398058252427, 0.544554455445545,      0.465346534653465, 0.525252525252525, 0.404040404040404, 0.418367346938776,      0.461538461538462, 0.368932038834951, 0.53125, 0.347826086956522,      0.485148514851485, 0.475247524752475, 0.444444444444444, 0.434210526315789,      0.424242424242424, 0.375, 0.447761194029851, 0.493975903614458,      0.357894736842105, 0.385416666666667, 0.565217391304348, 0.43956043956044,      0.348314606741573, 0.450549450549451, 0.358695652173913, 0.479166666666667,      0.518072289156627, 0.408163265306122, 0.458823529411765, 0.407894736842105,      0.520408163265306, 0.434782608695652, 0.323529411764706, 0.442105263157895,      0.44047619047619, 0.373626373626374, 0.456521739130435, 0.411764705882353,      0.443298969072165, 0.352941176470588, 0.383720930232558, 0.438356164383562,      0.529411764705882, 0.387755102040816, 0.428571428571429, 0.443298969072165,      0.383838383838384, 0.414141414141414, 0.393939393939394, 0.414141414141414,      0.484848484848485, 0.438775510204082, 0.412371134020619, 0.306122448979592,      0.494949494949495, 0.371134020618557, 0.438775510204082, 0.473684210526316,      0.444444444444444, 0.418367346938776, 0.459183673469388, 0.510416666666667,      0.285714285714286, 0.3125, 0.369565217391304, 0.375, 0.326923076923077,      0.38961038961039, 0.210526315789474, 0.324675324675325, 0.344086021505376,      0.360824742268041, 0.414141414141414, 0.5, 0.454545454545455,      0.46875, 0.414141414141414, 0.387755102040816, 0.505154639175258,      0.474226804123711, 0.363636363636364, 0.414141414141414, 0.454545454545455,      0.418367346938776, 0.285714285714286, 0.438775510204082, 0.387755102040816,      0.40625, 0.425531914893617, 0.515151515151515, 0.43010752688172,      0.402173913043478, 0.336734693877551, 0.459183673469388, 0.479166666666667,      0.296296296296296, 0.424242424242424, 0.505050505050505, 0.540816326530612,      0.458333333333333), ClaustPer = c(NA, 1L, 2L, NA, 3L, 0L, 2L,      0L, 1L, 0L, 0L, 0L, 1L, NA, 0L, 7L, 1L, 0L, 1L, 0L, 1L, 2L, 2L,      NA, 2L, 3L, 2L, 2L, NA, 0L, 1L, NA, NA, 0L, 0L, 0L, 0L, 3L, 3L,      3L, 1L, 0L, 2L, NA, 1L, 0L, 1L, 1L, 3L, 1L, 2L, 0L, 2L, 1L, 0L,      6L, 0L, 0L, NA, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, NA, 0L, 0L, NA,      NA, 0L, 0L, 0L, 0L, 0L, 1L, 0L, NA, 1L, 0L, 2L, 1L, 2L, 1L, 1L,      2L, 1L, 1L, NA, 1L, NA, NA, 2L, 1L, 1L, NA, 0L, 2L, 3L, 1L, NA,      1L, NA, 4L, NA, 0L, 1L, 1L, NA, 1L, 0L, 0L, 3L, 3L, 2L, NA, 1L,      0L, 0L, 1L, 1L, 0L, 0L, 1L, 2L, 4L, 1L, 1L, 0L, 0L, NA, 0L, 2L,      1L, 2L, 2L, 0L, 0L, 1L, 0L, 2L, 2L, 0L, 0L, NA, 2L, 2L, 4L, NA,      0L, 1L, 3L, 0L, 2L, NA, 2L, 3L, 2L, 2L, 3L, 2L, NA, NA, NA, 3L,      NA, 1L, 0L, 8L, 2L, 3L, 1L, NA, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 5L,      0L, 0L, 0L, 0L, 0L, 1L, 1L, NA, 1L, 3L, 1L, 0L, 1L, 1L, 0L, 0L,      NA, 2L, 2L, 2L, 3L, NA, 1L, NA, 2L, 1L, 13L, 2L, 1L, NA, 2L,      1L, NA, 2L, 2L, 1L, NA, 4L, 0L, 1L, 1L, 1L, 3L, 2L, 1L, 4L, 4L,      NA, 6L, NA, NA, 2L, 1L, 1L, NA, 2L, 0L, 2L, 1L, 6L, 2L, NA, 4L,      NA, 1L, NA, NA, 0L, 5L, 3L, 2L, NA, 2L, NA, NA, 4L, NA, 3L, 1L,      6L, 5L, NA, 0L, 1L, 0L, NA, 0L, 2L, 0L, 0L, 1L, 0L, NA, 2L, 1L,      NA, 1L, 0L, 3L, 2L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA)), .Names = c("ID", "MLH", "ClaustPer"), class = "data.frame", row.names = c(1L,      2L, 3L, 4L, 7L, 9L, 10L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 20L,      21L, 22L, 23L, 24L, 25L, 27L, 28L, 29L, 31L, 34L, 37L, 38L, 39L,      40L, 41L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 52L, 55L, 56L, 59L,      60L, 61L, 62L, 63L, 65L, 66L, 67L, 68L, 69L, 71L, 72L, 73L, 74L,      75L, 76L, 77L, 79L, 80L, 81L, 82L, 83L, 84L, 85L, 86L, 87L, 88L,      89L, 90L, 92L, 93L, 94L, 95L, 98L, 99L, 100L, 102L, 103L, 104L,      105L, 107L, 108L, 109L, 111L, 113L, 114L, 115L, 116L, 117L, 118L,      119L, 121L, 123L, 124L, 125L, 126L, 127L, 128L, 129L, 130L, 131L,      133L, 135L, 136L, 137L, 139L, 140L, 141L, 142L, 143L, 144L, 145L,      146L, 147L, 148L, 149L, 151L, 152L, 153L, 154L, 155L, 156L, 157L,      158L, 159L, 160L, 161L, 162L, 163L, 164L, 165L, 166L, 167L, 168L,      169L, 170L, 171L, 172L, 173L, 174L, 175L, 177L, 178L, 179L, 182L,      183L, 184L, 185L, 186L, 187L, 188L, 189L, 191L, 192L, 194L, 195L,      197L, 198L, 199L, 202L, 204L, 205L, 206L, 208L, 209L, 210L, 211L,      212L, 213L, 214L, 215L, 216L, 217L, 218L, 219L, 222L, 224L, 226L,      228L, 229L, 234L, 235L, 236L, 237L, 238L, 239L, 240L, 241L, 243L,      244L, 246L, 247L, 249L, 250L, 253L, 256L, 257L, 261L, 262L, 263L,      264L, 268L, 269L, 270L, 271L, 272L, 273L, 275L, 277L, 278L, 280L,      281L, 282L, 283L, 284L, 287L, 288L, 289L, 291L, 293L, 294L, 297L,      298L, 299L, 300L, 301L, 304L, 305L, 307L, 308L, 309L, 311L, 312L,      314L, 315L, 316L, 317L, 318L, 320L, 321L, 322L, 323L, 324L, 326L,      327L, 329L, 330L, 331L, 332L, 335L, 336L, 337L, 338L, 339L, 340L,      341L, 342L, 343L, 345L, 346L, 347L, 348L, 349L, 351L, 352L, 353L,      354L, 355L, 356L, 357L, 358L, 359L, 360L, 361L, 362L, 363L, 364L,      365L, 366L, 367L, 368L, 369L, 371L, 373L, 374L, 375L, 376L, 377L,      378L, 379L, 380L, 381L, 382L, 383L, 384L, 385L, 386L, 387L, 389L,      390L, 391L, 393L, 394L, 395L, 396L, 397L, 398L, 399L, 403L, 404L,      405L, 406L, 407L, 408L, 409L, 410L, 411L, 412L, 413L, 414L, 415L,      416L, 417L, 418L, 419L, 420L, 421L, 422L, 423L, 424L, 425L, 428L,      429L, 432L, 436L, 439L, 440L, 441L, 442L, 444L, 445L, 447L, 449L,      450L, 451L, 452L, 453L, 454L, 455L, 456L, 457L, 458L, 459L, 460L,      461L, 462L, 463L, 464L, 465L, 466L, 469L, 470L, 471L, 472L, 473L,      474L, 475L, 476L, 477L, 478L, 479L, 480L, 481L, 482L, 483L, 484L,      485L, 486L, 487L, 488L, 489L, 490L, 491L, 492L, 493L, 494L, 495L,      496L, 497L, 498L, 499L, 501L, 502L, 503L, 504L, 505L, 507L, 509L,      510L, 511L, 512L, 513L, 514L, 515L, 516L, 517L, 518L, 519L, 520L,      521L, 522L, 523L, 524L, 525L, 526L, 527L, 528L, 529L, 530L, 531L,      532L, 533L, 534L, 535L, 536L, 537L, 538L, 539L, 540L, 541L, 542L,      543L, 544L, 545L, 546L, 547L, 548L, 549L, 550L, 551L, 552L, 553L,      554L, 557L, 558L, 559L, 560L, 561L, 562L, 563L, 564L, 565L, 566L,      567L, 568L, 569L, 570L, 571L, 572L, 573L, 574L, 575L, 576L, 577L,      578L, 579L, 580L, 581L, 582L, 583L, 584L, 585L, 586L, 587L, 588L,      589L, 590L, 591L, 592L, 593L, 594L, 595L, 596L, 597L, 598L, 599L,      600L, 601L, 602L, 603L, 605L, 606L, 607L, 608L, 609L, 610L, 611L,      612L, 613L, 614L, 615L, 616L, 617L, 618L, 619L, 620L, 621L, 622L,      623L, 624L, 625L, 626L, 627L, 628L, 629L, 630L, 632L, 633L, 634L,      635L, 636L, 637L, 638L, 639L, 640L, 641L, 643L, 644L, 645L, 646L,      647L, 648L, 649L, 650L, 651L, 652L, 653L, 654L))          summary(data$ClaustPer)         hist(data$ClaustPer)          summary(data$MLH)         hist(data$MLH)


title: How to calculate group statistics in R

text: > **Possible Duplicate:**   >  R: compute correlation by group I have a data set with a 6 or so data points for groups of data. e.g.:               email_id   date_sent  number_sent  number_of_views  number_of_responses     1           5/4         600            25                6     1           5/5         500            22                8     1           5/6         450            23                4     1           5/7         700            34               12     2           5/5         900            30               10     2           5/6         750            28               11     ...      (this is made up data that illustrates the point) Assuming I have this in a data frame in R, I'd like to write something which will give me stats by group. I'm most interested in the correlation coefficient between some of the columns. I know how to do this with a data frame that contains only one group: cor(col1, col2) but I'd like to learn a technique that will allow me to extract data that looks something like this:               email_id    cor(col3, col4)     1             .73     2             .85     3             .98      and so on. Thanks,


title: How do we know that a classifier is linear or not?

text: I searched on the net, but I could not find any really good answer to my question which is very easy. A linear classifier is a classifier which can classify the samples correctly. But Given a function, like for example an f that uses sigmoid, and classifies to girl/boy for example if the result of sigmoid is less or more than 0.5, how do we know if f is actually a linear function? More generally, given a function (and not the data), how we are able to know if that function is a linear classifier or not? Should we prove something? Are there any references/bibliography that can help me understand better? Any other example that can help me better?


title: Averaging time series with different sampling interval

text: I have a few time series that were (for technical reasons) acquired with slightly different time intervals, ranging between 19 and 21 seconds. Now, I would like to average the values of these different time series over time, so I thought that I could do some sort of interpolation of the values at some regular interval (e.g. each 20 seconds). Could someone point out a good way to do this? Thank you nico


title: SAS: How to create variables from the values given in a specific column of a table?

text: ![http://i.stack.imgur.com/hWSze.png](http://i.stack.imgur.com/hWSze.png) I'm using SAS and have come across a problem. I have a dataset with some variables / columns, one of which has values of only the following: `blue`, `red`, and `gray`. How can I convert these values into respective variables and assign the `count` as values, as shown in the second table?


title: Evaluating an integral using MCMC algorithm

text: Consider the problem of evaluating the following integral $$ \int_T f(\theta_1, \theta_2, \theta_3| D) d\Theta$$ where $f$ is a posterior density and the values of $D = (d_1, ..., d_n)$, given $\Theta$ are $i.i.d$ with a given distribution and $T$ has the following form: $$T = \\{(\theta_1, \theta_2, \theta_3) \in M | f(\theta_1, \theta_2, \theta_3|D) > k \\} $$ for $k$ a known constant bigger than zero and $M=(0,\infty)$x$(1,\infty)$x$(0,\infty)$ I don't know how to apply $MC$ integration method to solve this integral, since I can't see how is this region $T$ or even how to write its volume. I thought we could use the distribution of $d|\Theta$, but I don't know how to. I would be glad if I could have some help understanding better how to calculate this integral Thanks! (P.S: I will use the metropolis hastings to sample from $f$, thats why I wrote MCMC)


title: What is the best technique for panel regressions (ols, fixed, between, random effects)?

text: I have panel data which includes American states (1-48) and years (1900-1917). All the variables are time-varying with one exception. This exception is time invariant and a three level categorical variable measuring regional designations for the states tested using two dummy variables. I also want look at interactions between one of the dummies and several of the time-varying variables. OLS estimates of this model indicate heteroscedasticity. So given the organization of the data my question is which of the panel techniques is best and why? I am using Stata. Thanks, RB


title: Identifying outliers based on standard error of residuals vs sample standard deviation

text: If this is too elementary a question for CrossValidated and should go elsewhere please let me know. I've been given a linear model in R similar to:               model = lm(y ~ x1 + x2 + x3, data=sampleset)      ...that when fed a set of data and plotted produces the graph below: ![Predicted vs observed](http://i.stack.imgur.com/1Oya8.png) The light blue center line represents the linear regression prediction. The blue area represents 2 standard deviations from expected in either direction, using "standard error of residuals":               # se_: number of standard errors away from the predicted value     se_thresh = 2          # standard error of residuals     res_se = sd(model$residuals)          # for each member in the sample     if(model$fitted.values[j] < (i - se_thresh*res_se)){           dataset$mark_good[j] = TRUE     }      I've been told that we can infer which members of the sample are unusually "good" or "bad" on the basis of whether the actual scores are 2 SEs above or below the predicted value. I'm wondering:   * why comparing members on the basis of the SE of residuals is better/more valid than simply comparing n with the standard deviation, which is simpler and would not require the model?   * what standard error of residuals represents intuitively in the model


title: lsqnonneg – variation between results in R and MATLAB

text: We are trying to do nonlinear least squares fitting in R. We have reference code in MATLAB as below.                A =  [ -4.09549023 -1.5924967 -5.2775267 0.7195365 -1.4681932;      -0.09302291  0.2538085  0.6080219 0.1413484 -0.4614447;         1.00000000  1.0000000  1.0000000 1.0000000  1.0000000;]           b = [-0.52701539;0.08974224;1.00000000]          [w,res] = lsqnonneg(A, b)          w =                   0              0         0.1377         0.6700         0.1922               res =             1.2519e-32      If we try the same thing in R, the results are entirely different.               A <- rbind(cbind(4.09549023, -1.5924967, -5.2775267, 0.7195365, -1.4681932),     cbind(-0.09302291,  0.2538085,  0.6080219, 0.1413484, -0.4614447),      cbind(1.00000000,  1.0000000,  1.0000000, 1.0000000,  1.0000000))          b <- rbind(-0.52701539,0.08974224,1.00000000)          sol <- lsqnonneg(A,as.vector(b))          $x     [1] 0.2212015 0.1986838 0.1890215 0.2080815 0.1830117          $resnorm     [1] 4.889664e-10          $residual                   [,1]     [1,]  2.871547e-11     [2,]  4.866847e-10     [3,] -3.743827e-11      The problem is, though the input values of `A` and `b` are the same, and the methods invoked both in R and MATLAB have the same purpose (viz., non-negative least squares fitting) the residuals are different in both cases. `resnorm` = 4.889664e-10 in R and `res` = 1.2519e-32 in MATLAB. If anyone can point out some obvious mistake in the results, or suggest some alternative method in R to match the method in MATLAB, it would be really appreciated.


title: Including nested random effect in a GAM

text: We have 6 transects at 6 different sites (6 transects per site). Each site subdivided into two zones therefore in our model we have: Site (fixed categorical factor) and Zone (random categorical factor nested within Site). During each transect we measured occurrence of diseases (our continuous response variable) and benthic cover (our continuous covariate). To test the potential relationship between disease prevalence and benthic coverage we used a generalized additive model (non parametric distribution of our response variable and no a priori assumptions on the shape of the relationship). I have tried many models on R but to be honest, I am a bit lost as it seems to me that there are infinite possibilities. First I ran the following model: model_1<-gam(Occurrence ~ s(Cover) + Site + Zone, data=df) but I am not sure whether I should also use smoothing splines on the categorical factors (like in model_2). I have seen it in other posts but I am wondering if it makes any sense. model_2<\- gam(Occurrence ~ s(Cover) + s(Site) + s(Zone), data=df) Neither model_1 nor model_2 take into account the fact that Zone is nested within Site and is a random factor though. I have tried this: model_3<\- gam(Occurrence ~ s(Cover) + s(Site) + s(Zone,bs=”re”), data=df) model_4<-gam(Occurrence ~ s(Cover)+ Site + s(Site, Zone, bs="re"),data=df) I was advised to use a generalized additive mixed model since I had both fixed and random factors in my design, this did not make things easier. Would you have any advice on what model would suit best and how to include properly Zone random and nested within Site? Of course I am aware that there are many other parameters to consider (family distribution, link function, knots etc.) but for now I would like to limit the question. These issues are probably confusing for a lot of people new to GAM and hopefully your help will be useful to numerous users.


title: Data fitting with repeated measurements

text: we have an experiment to excite a system with some energy, then measure the decay as a function of time. We are measuring data at 4 times $t$ to fit to an exponential model: $y = a \exp(-t/p) + c$, where we fit the parameters $a, p$, and $c$. I should probably mention that there is not one excitation and 4 measurements, but rather we excite 4 times, and adjust our time-from- excitation to measure -- so each measurement is taken independently. We have pretty good signal-to-noise with respect to the measurement device, however the object under consideration may or may not move during the experiment resulting in some unpredictable error in one or more data points. Parameter $p$ is the most important, and we know a physical range a-priori for this parameter. We must make this measurement in a fixed amount of time and we have time to collect a maximum of 8 measurements. Right now, we are collecting the data twice at each of the 4 time points, so we have basically two distinct experiments. I have a lot of questions, actually :), but one I find interesting is whether it is better to collect two measurements at 4 time points, or if it is better to collect 8 measurements at different time points (supposedly to better characterize the decay curve). We've done the experiment at up to 16 time points and found good agreement with the 4 time point experiment. Other questions I have pertain to how to detect bad measurements in the 2X 4 time point data sets. Any ideas would be greatly appreciated! I am not too well-learned in statistics, so if you have some suggestions on some background reading or material pertinent to this problem, that would be awesome! Thanks!!


title: Statistical differences between two percentages

text: I am doing an experiment on mice. I found 11 recombinants within 142 mice. I calculated the recombination frequency and found it to be 7.7%. Other studies (many studies, not one or two) had 1000 mice and 58 recombinants with recombination frequency 5.8%. My question is whether 7.7% is significantly different from 5.8%.


title: How to test if two coins have different biases?

text: probably a noob question because I'm a noob. I have two coins, A and B. For each coin I have a sample of the results I get by tossing it. My null hypothesis is that A has equal or higher probability to land on HEAD than B, but that probability is unknown (as well as B's probability). However, the data shows me that B seems to be more probable to land on HEAD. How can I test the null hypothesis? It seems to me that a simple t-test won't work - A coin toss is Bernoulli distributed, and not normally distributed. I thought about dividing each sample to sub-samples, consider the ratio of HEADs in each subsample as a random variable that is distributed approximately like a normal variable, and then use t-test. However, this seems like a complicated (and inaccurate) solution for a basic (the most basic?) question. bonus points: what's the effect size?


title: What statistics should I use to compare percentages?

text: I am trying to see if there is a statistical significant difference between the amount of bycatch in 5 different types of gear used by fishermen. For each gear I have varying number of replicates as there are higher number of fishermen using certains types of gears compared to others. I have converted my raw data into 'percentage of time occurred' in order to compare them. What stats test should i use? I have access to minitab, but I am a novice with it.


title: How do I interpret the results from a basic interaction plot from the R effects package?

text: Using the method in this post, I have made a plot to visualize the interaction between two predictor variables using the effects package in r, but I'm not really sure what I am looking at. Tide heights and rain averages are continuous. 8 bins were the maximum the function would allow me to use. The following is the call to `effect` producing this plot:                R > plot(effect(term="rain.avg2:tide.avg",mod=bkrain9.lm,default.levels=8),               main="", xlab="Precipitation - 24hr average (cm)",               ylab=expression("TCB Concentration - CFU*100m"*L^-1),multiline=TRUE)      **Plot updated.** Could somebody explain the purpose of this plotting feature in context to predictor interactions? ![P](http://i.stack.imgur.com/neZql.jpg) Background: For a class project, I have created a linear regression model to evaluate the effects of the interaction between two predictor variables (tide height and precipitation) on bacteria concentrations. Thermo-tolerant coliform bacteria concentrations were sampled at 5 sites in a day, where a sample time was recorded at sampling completion. I took an average of these for roughly 20 days, and calculated an associated 24hr average of precipitation before sampling completion, and a 50 minute (the sampling duration) average of tide height before sampling completion.


title: Multiple imputation and model selection

text: Multiple imputation is fairly straightforward when you have an _a priori_ linear model that you want to estimate. However, things seem to be a bit trickier when you actually want to do some model selection (e.g. find the "best" set of predictor variables from a larger set of candidate variables - I am thinking specifically of LASSO and fractional polynomials using R). One idea would be to fit the model in the original data with missing values, and then re-estimate this model in MI datasets and combine estimates as you normally would. However, this seems problematic since you are expecting bias (or else why do the MI in the first place?), which could lead to selecting a "wrong" model from the start. Another idea would be to go through whatever model selection process you are using in each MI dataset - but how would you then combine results if they include different sets of variables? One thought I had was to stack a set of MI datasets and analyze them as one large dataset that you would then use to fit a single, "best" model, and include a random effect to account for the fact you are using repeated measures for each observation. Does this sound reasonable? Or perhaps incredibly naive? Any pointers on this issue (model selection with multiple imputation) would be greatly appreciated.


title: is there any package / simulator for discrete time markov chain simulation

text: is there any package or simulator for discrete time markov chain simulation. Matlab/ Python solutions ?


title: confusion related to higher order markov chain

text: I was reading this book related to machine learning. It's given that for Mth order markov chain, the number of parameters = $K^{M-1}(K-1)$ where M is the order. I am not sure how this is derived. For eg, lets say I have three states {Sunny, Cloudy, Rainy}. So at each time stamp I will look at the previous two states. So if I use the above formula, I will have 3*(2-1)*(3-1) = 6 which is less I guess. I should have the following parameters               P(Sunny|Cloudy,Rainy)     P(Sunny|Cloudy,Sunny)     P(Sunny|Cloudy,Cloudy)     P(Sunny|Rainy,Rainy)     P(Sunny|Rainy,Sunny)     P(Sunny|Rainy,Cloudy)     P(Sunny|Sunny,Rainy)     P(Sunny|Sunny,Sunny)     P(Sunny|Sunny,Cloudy)      and so on. It has lots of parameters isn't it higher than 6. I am bit confused. IT should be $K^{M}(K-1)$ I guess


title: How to analyze categorical variables with multiple levels

text: I have 14 different habitat classes and two activity states (so two variables – habitat and activity). For activity state A, I have a data count of over 300, but with B, I only have around 30 (so unequal sample sizes). I really want to see if one or more habitat/s is used significantly more than others during each activity state, but then I also want to compare habitat count between activity states, to see if they are significantly different in this regard. I have two other variables, site location and gender, and I would also like to see if the habitat count varies significantly between these two. Obviously it would also be nice to combine some of the variables, i.e., does habitat count vary between gender between sites, if this makes sense. I am currently using JMP8 as I have very limited experience with R. I have been using chi squared and contingency analysis (although a number of my cells have quite low counts, less than five, which chi squared struggles with). What would be the best way to approach this statistically? How would I need to format the data (it is currently in text rather than count form, as there are so many classes in the habitat variable that it seems otherwise difficult to manage)? Any help would be hugely appreciated!!


title: Do I need P value Adjustment for 100 pairwise comparisons?

text: I made 100 comparisons by applying chi-square test (tajima relative rate test) and got 100 p values. Do I need P-value adjustment? If yes then which one correction will be suitable? I have explored internet and also on this forum but didnt get answer of this particular question.


title: Poisson vs Binomial for rare events

text: From Poisson's postulates, we know Poisson works for rare events. However, we also know binomial is an approximation of Poisson when the probability of an event is small. So can we use binomial and Poisson interchangeably for rare event? What is the benefit of using one rather than the other?


title: Variable selection for time covariate

text: I'm fitting a linear model where the response is a function both of time and of static covariates (i.e. ones that are independent of time). The ultimate goal is to identify significant effects of the static covariates. Is this the best general strategy for selecting variables (in R, using the `nlme` package)? Anything I can do better?   1. Break the data up by groups and plot it against time. For continuous covariates, bin it and plot the data in each bin against time. Use the group-specific trends to make an initial guess at what time terms to include-- time, time^n, sin(2*pi*time)+cos(2*pi*time), log(time), exp(time), etc.   2. Add one term at a time, comparing each model to its predecessor, never adding a higher order in the absence of lower order terms. Sin and cos are never added separately. **_Is it acceptable to pass over a term that significantly improves the fit of the model if there is no physical interpretation of that term?_**.   3. With the full dataset, use forward selection to add static variables to the model and then relevant interaction terms with each other and with the time terms. **_I've seen some strong criticism of stepwise regression, but doesn't forward selection ignore significant higher order terms if the lower order terms they depend on are not significant? And I've noticed that it's hard to pick a starting model for backward elimination that isn't saturated, or singular, or fails to converge. How do you decide between variable selection algorithms?_**   4. Add random effects to the model. **_Is this as simple as doing the variable selection using`lm()` and then putting the final formula into `lme()` and specifying the random effects? Or should I include random effects from the very start?_**. Compare the fits of models using a random intercept only, a random interaction with the linear time term, and random interaction with each successive time term.    5. Plot a semivariogram to see if an autoregressive error term is needed. **_What should a semivariogram look like if the answer is 'no'? A horizontal line? How straight, how horizontal? Does including autoregression in the model again require checking potential variables and interactions to make sure they're still relevant?_**   6. Plot the residuals to see if the variance changes as a function of fitted value, time, or any of the other terms. If it does, weigh the variances appropriately (for `lme()`, use the `weights` argument to specify a `varFunc()`) and compare to the unweighted model to see if this improves the fit. **_Is this the right sequence in which to do this step, or should it be done before autocorrelation?_**.   7. Do `summary()` of the fitted model to identify significant coefficients for numeric covariates. Do `Anova()` of the fitted model to identify significant effects for qualitative covariates.


title: Detecting noisy patterns in document images

text: I am looking for features to extract to distinguish between text objects and arbitrary noisy patterns in degraded document images. ![enter image description here](http://i.stack.imgur.com/zTKbx.png) This is an example of a document with some parts of noise, I have extracted manually 4 windows, three of them correspond to noise, and the last one to text. ![noise1](http://i.stack.imgur.com/cPjop.png) ![noise2](http://i.stack.imgur.com/L9jHM.png) ![noise3](http://i.stack.imgur.com/SJRHc.png) ![text](http://i.stack.imgur.com/zJBvf.png) Any suggestion ?


title: handling multicollinearity by backwards regression and omitted variable bias

text: Suppose I try to estimate a production function as follows: $logY=b_0+b1*logX_1+b_2*logX_2+b_{11}*(logX_1)^2+b_{22}*(logX_2)^2+b_{12}*(logX_1)*(logX_2)+u$, where $Y$ is the output, $X_1$, $X_2$ are inputs, and for simplicity, $E(u|X_1, X_2)$ is assumed. $b_{ij}$ are parameters to be estimated. Since the log-inputs and their interactions are usually highly correlated, the multicollinearity is likely to be a serious problem in the sense that lots of coefficients could be insignificant under conventional alphas. To fix this problem, I follow the routine to eliminate the least significant variables backwards. However, I'm not sure whether this practice is fully legitimate, especially when the coef. of remaining variables change sharply. This is because (though without convincing arguments) I tend to associate the "large" changes with the Omitted Variable Bias (OVB) in econometrics literature (roughly speaking, by omitting a variable belonging to the underlying population model, the new error term could be correlated with covariates, thus estimates could be biased). So my questions about the backwards regression and the OVB are: (1) Is there any theory to identify "large" changes in estimates? Whether increase from 0.5 to 0.6 (20%) can be viewed as a large change? What about a decrease from 0.5 to 0.1? (2) To what extent should I rely on the imprecise estimates to calculate such changes? (3) Whether I should associate "large" changes in estimates using sample data with dropping variables belonging to underlying population model? Thank you !


title: Best Regression with Binary Features

text: I'm seeking to do a linear regression for an evaluation function in a board game. My features are all (signed) binary 1 0 0 1 -1 1 0 0 0. Mostly zeros. Around 200 to an observation. I have 10 million observations with target values [-1000, 1000]. I'm wondering how this should affect my choice of regression model. I've tried with least squares and it gave very bad and unstable results. Can you recommend me some good, linear, robust regression models, which perhaps take advantage of the simple format of my inputs?


title: Bayesian Survival Analysis: please, write me a prior for Kaplan Meier!

text: Consider right-censored observations, with events at times $t_1, t_2, \dots$. The number of susceptible individuals at time $i$ is $n_i$, and the number of events at time $i$ is $d_i$. The Kaplan-Meier or product estimator arises naturally as a MLE when the survival function is a step function $S(t) = \prod_{i : t_i < t} \alpha_i$. The likelihood is then $$ L(\alpha) = \prod_i (1-\alpha_i)^{d_i} \alpha_i^{n_i-d_i} $$ and the MLE is $\widehat\alpha_i = 1 - {d_i\over n_i}$. OK, now assume that I want to go Bayesian. I need some kind of ``natural'' prior with which I will multiply $L(\alpha)$, right? Googling the obvious keywords I found that the Dirichlet process is a good prior. But as far as I understand, it is also a prior on the discontinuity points $t_i$? This is surely very interesting and I am eager to learn about it, however I would settle for something simpler. I begin to suspect that’s not so easy as I first thought, and it’s time to ask for your advice... Many thanks in advance! **PS: A few precision on what I am hoping** I am interested in (as simple as possible) explanations about the way to handle the Dirichlet process prior, however I think it should be possible to use simply a prior on the $\alpha_i$ — that is a prior on the step functions with discontinuities in $t_i$. I think that the "global shape" of the step functions sampled in the prior should not depends on the $t_i$'s -- there should be an underlying family of continuous functions which are approximated by these step functions. I don't know if the $\alpha_i$ should be independent (I doubt it). If they are, I think this implies that the prior $\alpha_i$ depends on$\Delta t_i = t_i - t_{i-1}$, and if we denote its distribution by $A(\Delta t)$ then the product of a $A(\Delta_1)$ variable by an independent $A(\Delta_2)$ variable is a $A(\Delta_1+\Delta_2)$ variable. It seems here that log-$\Gamma$ variables can be useful. But here basically I am stuck. I did not type this at first because I did not want to direct all answers in this direction. I would particularly appreciate answers with bibliographic references to help me justifying my final choice.


title: Transform LN variates to match reference mean and SD

text: I have some measurements of a biological analyte that exhibits a log-normal distribution. There was both a level shift and some long-term drift in the values that I believe is due to the measurement process rather than to biology. I would like to use recent values as the reference distribution and transform within months or quarters so that the mean(log(Y)) and sd(log(Y)) will match. Experimentation makes me think I can do the transformation with:               exp( sd(log(X_ref))/sd(log(Y) *( log(Y) -mean(log(X_ref) ) ) -> Y_transformed      Can this be supported in theory as opposed to experiment?


title: Logistic regression with LBFGS solver

text: Is there any open source library or code which implements Logistic Regression using L-BFGS solver? I would prefer Python, but other languages are welcome, too.


title: A mathematical formula for K-fold cross-validation prediction error?

text: Could anyone elaborate on a theoretical perspective of K-fold cross-validation and especially a mathematical formula for k-fold CV prediction error? **Updated** :Could any one of you help me to understand the formula written on latest edition (Feb, 2011) of the book at page 242. ![enter image description here](http://i.stack.imgur.com/a9rE0.jpg)


title: Ignorability in Rubin's theory of missing data mechanisms

text: I am trying to understand Rubin's theory of bayesian inference with missing data, specifically how the missing data mechanism affects the inference on a superpopulation parameter. The theory is exposed for example in the chapter 7 of [1]. The complete $N$-dimensional data vector $y$ is splitted into an observed and unobserved components: $y = (y_{\mathrm{obs}}, y_{\mathrm{nob}})$ and $I$ is the $N$-vector of indicators whose components $I_{i}$ equal 0 when $y_{i}$ is unobserved and 1 when it is observed. The sample space is thus the product of the $y$ space and the $I$ space. For simplicity, I consider a missing data mechanism $p(I ~|~ y_{\mathrm{obs}})$ which does not depend on a fully observed covariate $x$ or an unknown parameter $\phi$. The key result in Rubin is that when $p(I ~|~ y) = p(I ~|~ y_{\mathrm{obs}})$, the missing data mechanism is ignorable and the analysis can proceed as usual. This can be seen in the following development: $$ \begin{align*} p(\theta ~|~ y_{\mathrm{obs}}, I) &~=~ \frac{p(\theta) \int p(y~|~\theta)\, p(I ~|~ y)~d y_{\mathrm{nob}}~}{\int\int p(y~|~\theta)\,p(\theta) ~d\theta p(I ~|~ y)~d y_{\mathrm{nob}}~} \\\ \end{align*} $$ If $p(I ~|~ y) = p(I ~|~ y_{\mathrm{obs}})$, it can be taken out of the integrals and cancelled out such that the posterior distribution $p(\theta ~|~ y_{\mathrm{obs}}, I) = p(\theta ~|~ y_{\mathrm{obs}})$. So far so good. But what I don't understand is the following derivation that contradicts the preceding result: $$ \begin{align*} p(\theta ~|~ y_{\mathrm{obs}}, I) &~=~ \frac{p(y_{\mathrm{obs}}, I, \theta)}{p(y_{\mathrm{obs}}, I)} \\\ &~=~ \frac{p(y_{\mathrm{obs}} ~|~ \theta) p(\theta) p(I ~|~ y_{\mathrm{obs}})}{\int p(y_{\mathrm{obs}} ~|~\theta)p(\theta) ~d\theta~ p(I ~|~ y_{\mathrm{obs}})} \\\ &~=~ \frac{p(y_{\mathrm{obs}} ~|~ \theta)p(\theta)}{\int p(y_{\mathrm{obs}} ~|~\theta)p(\theta) ~d\theta} \\\ &~=~ p(\theta ~|~ y_{\mathrm{obs}}). \end{align*} $$ Here the sampling mechanism cancels out no matter what. It does not matter if $p(I ~|~ y) = p(I ~|~ y_{\mathrm{obs}})$. I must be doing something wrong, but I can't see what. [1] Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B. (2004). Bayesian Data Analysis (2nd ed.). Boca Raton: Chapman & Hall/CRC.


title: What is the distribution of norm induced by an inverse Wishart?

text: Suppose $S$ is distributed as a Wishart matrix with $n$ degrees of freedom and scale matrix $\Sigma$, and let $\vec{a}$ be a fixed vector. It is well known that $\vec{a}^{\top}S\vec{a}$ is equal to $\vec{a}^{\top}\Sigma\vec{a}$ times a Chi-square random variable with $n$ degrees of freedom. I am curious about the distribution of $\vec{a}^{\top}S^{-1}\vec{a}$. My guess was that it would be disributed as an inverse gamma, with scale $\vec{a}^{\top}\Sigma^{-1}\vec{a} / 2$, and shape $n/2$. My simulations (which are somewhat questionable) indicate this is _not_ the case. Does this follow a well-known distribution?


title: How to know which variables are more important in a process?

text: I have a process with 15 effective variables. I could record 9 variables to study its effect on process.   I am looking for an appropriate factor to estimate the value of effectiveness of each factor. I believe p-values, standardized coefficients, partial correlation and partial $R^2$ are the potential solutions in this regard. Which one is more useful in informing the degree of effectiveness of each variable on scale of 0 to 100?


title: Recalibrating weights from complex survey with external data

text: Hello CrossValidated Users, I am using HCUP NIS data to determine various effects on health outcomes using SAS software. I recognize that this is a complex survey and that the design variables need to be incorporated using PROC SURVEYlogistic (etc.). Usually we are interested in subpopulations so the DOMAIN statement is necessary. I also read that HCUP has a way to make dummy values for hospitals to achieve the proper standard errors. Instead, can I calculate different weights based on existing distributions from external national data on my particular disease cohort?


title: Can I use Kolmogorov-Smirnov to compare two empirical distributions?

text: Is it OK to use the Kolmogorov-Smirnov goodness-of-fit test to compare two empirical distributions to determine whether they appear to have come from the same underlying distribution, rather than to compare one empirical distribution to a pre-specified reference distribution? Let me try asking this another way. I collect N samples from some distribution at one location. I collect M samples at another location. The data is continuous (each sample is a real number between 0 and 10, say) but not normally distributed. I want to test whether these N+M samples all come from the same underlying distribution. Is it reasonable to use the Kolmogorov- Smirnov test for this purpose? In particular, I could compute the empirical distribution $F_0$ from the $N$ samples, and the empirical distribution $F_1$ from the $M$ samples. Then, I could compute the Kolmogorov-Smirnov test statistic to measure the distance between $F_0$ and $F_1$: i.e., compute $D = \sup_x |F_0(x) - F_1(x)|$, and use $D$ as my test statistic as in the Kolmogorov-Smirnov test for goodness of fit. Is this a reasonable approach? (I read elsewhere that the Kolmogorov-Smirnov test for goodness of fit is not valid for discrete distributions, but I admit I don't understand what this means or why it might be true. Does that means my proposed approach is a bad one?) Or, do you recommend something else instead?


title: How do I visualize changes in proportions compared to another period?

text: I have a hypothetical set of diet composition, showing how much of each fruit does a particular group of people eat something like: ![table percentages](http://i.stack.imgur.com/jAyxP.png)   * How can I visualize in one page, all the compositions for groups A to G and more importantly, their change over the past year?    * I know having more than 4 colors in design is a no-no, but what can I do here?


title: Meaning of the "existence" proof

text: After doing some reading on stochastic processes for work I've found that a proof that the specific process _exists_ is often one of the first things presented. Could someone please explain, in layman's terms, the purpose/necessity of this proof?


title: Bayesian estimation of a mean with tighter constraints than the observations

text: The basic problem I have is very similar to a classic Bayesian estimation problem: there is a real valued parameters $p\in[0,1]$, that I would like to estimate via observations $o_i$. It is also the case that $$ p=\lim_{N\to\infty}{1\over N}\sum_{i=1}^N o_i, $$ so this problem is very closely analogous to estimating a probability. Here is the catch: $o_i\in[-1,1]$ Due to the similarity with estimating a probability, I would like to do something like a Bayesian update about my knowledge of $p$, where the prior would be a Beta distributed RV. However, since $p\in[0,1]$ it would seem like the prior would only have support in $[0,1]$, but that would mean that any $o_i\in[-1,0]$ would have not lead to any useful information about $p$, which is simply incorrect. On the other hand, if I have the prior have support over the whole $[-1,1]$ interval and perform the usual Bayesian update, the posterior distribution will indicate $\Pr(p\in[-1,0])\not=0$, which is also incorrect. What is the appropriate prior and update rule in this case? Has anyone seen problems similar to this one?


title: What's a good introduction to simulated method of moments and the extended path technique?

text: I'm reading a paper by Stephane Adjémian on DSGE modeling with a zero lower bound for the nominal interest rate, and he's using what he describes as the simulated method of moments / extended path. Has anyone worked with these techniques? What would you say is the next step toward gaining familiarity with them for someone who has a bit of a background in GMM estimation. I know that the main paper is McFadden (1989), but does anyone know of a textbook treatment of this material? I'd like to avoid having to worry about probability theory if possible.


title: How to show that $E\left[ y E(y|x) \right] = E\left[ E(y|x)E(y|x) \right]$ - linearity of conditional expectations?

text: I'm trying to finish a proof for a review exercise and I'm asked to show that $E\left[(y-E(y|x))(E(y|x)-f(x))\right]=0$ where $y$ is the dependent variable and $f(x)$ is a linear predictor of $y$. I'm almost finished, but I just want to check whether or not $E\left[ E \left[ y E(y|x)\,|\,x \right] \right]=E\left[ E(y|x)E \left[ E(y|x)\,|\,x \right] \right]$ Basically, can I take $y$ out of the conditional expectation as $E(y|x)$, or worded differently, are conditional expectations multiplacative in this way? Ordinarily you can back a function out of the expectation if it is a function of the variable being conditioned on - i.e., $E\left[ f(x)y|x \right]=f(x)E(y|x)$, but I think what I'm trying to do above is different. Also, a second, related question – if $E\left[E(y|x)\right]=E(y)$ by the Law of Total Expectations, then does $E\left[E(y|x)E(y|x)\right]=E(y^2)$ by the same token?


title: Heckman selection model with difference-in-differences specification

text: Following my question on Tobit with DiD specification I am wondering if it is possible to estimate a heckman sample selection model with a Difference in Differences specification? For example in STATA, using the heckman command how do I implement a DiD specification in the selection equation and how do I interpret my results?


title: Find a general normal stationary process

text: I am wondering how to find the general normal stationary process satisfying $X_{n+2} + X_{n} = 0$. Any help would be much appreciated, although I am relatively new to this space so some details how to do this in general would be very nice. Alternatively please feel free to point out relevant literature online.


title: Calculate prediction interval for ridge regression?

text: Does anyone know how to compute the prediction interval for ridge regression? and what's its relation with prediction interval of OLS regression?


title: Simulate dependent variable using a linear regression and a correlation coefficient

text: I create a vector $X$ with $n$ values and I would like to generate another vector $Y$ such as $Y = aX+b$ and add some noise to $Y$ through a given correlation coefficient (e.g. $r=0.8$). It is essentially a simple linear regression. My question is can I find the standard deviation of the noise from a correlation coefficient without using a brute force approach like I did below in `R`?               r <- 0.8     slope <- 0.003     intercept <- 0     n <- 100     X <- runif(n, 50, 100)          rnormcor <- function(x,sigma) rnorm(1,x,sigma)          f <- function(x, indpt,r,b1,b0){       y <- indpt*b1+b0       y <- sapply(y, rnormcor, x)       abs(r-cor(indpt,y))     }     sigmas = seq(0.0000001,1,by=0.001)     sigmas.cor = sapply(sigmas, f, indpt=X, r=rho, b1=slope,b0=intercept)     plot(sigmas , sigmas.cor)     sigmas[which.min(sigmas.cor)]


title: r, How to create two correlated variables that are distributed jointly normal (mean 0, var 1)

text: I want to create two random variables $X \sim N(0,1)$ and $Y \sim N(0,1)$ that satisfy $E(X,Y)=0.5 $ That is I want to create $Z=(X,Y)^\top $ with a joint bivariate normal distribution $ Z \sim N\left( \left(\begin{array}{c} 0\\\ 0 \end{array}\right) , \left(\begin{array}{cc} 1 & 0.5\\\ 0.5 & 1 \end{array}\right) \right) $. How do I code this in R? Thanks.


title: How to generate two negatively correlated random series of numbers?

text: > **Possible Duplicate:**   >  How to generate correlated random numbers (given means, variances and > degree of correlation)? I have been trying to generate two random variables with correlation of -0.9 but have not been successful. I have tried:               x1 <- rnorm(200, mean=0, sd=1)     x2 <- rnorm(200, mean=0, sd=1)      z <- -0.9*x1+sqrt((1-(-0.9)))*x2     cor(x1,z)      But it did not work. I tried also the ?mvrnorm function, but it did not work either. How I can get this done??


title: Comparing two means: Same group, different variance (Welch's t-test?)

text: **Intro:** I have recently made a simulation program which simulates patients with Type 1 diabetes. In that context I am creating artifical patients. Let's denote three of these as $\text{p}_1$, $\text{p}_2$, and $\text{p}_3$. Assume that $\text{p}_1$, $\text{p}_2$ and $\text{p}_3$ are treated with treatment $A$. Treatment $A$ treats $\text{p}_1$, $\text{p}_2$ and $\text{p}_3$'s diabetes to a certain degree. We evaluate the treatment with some measures conducted in a virtual trial; let's denote one of these measures $\text{M}1_A$. $\text{M}1_A$ has an average value and a standard deviation for Treatment A for the three patients. Now let's introduce treatment B for the same patient group, $\text{p}_1$, $\text{p}_2$ and $\text{p}_3$. Again we put the treatment through a trial test and receive a list of measures, this time denoted $\text{M}1_B$. I am interested in testing whether the difference in $\text{M}1_A$ and $\text{M}1_B$ is significant. I know: $N=3, \text{std(M}1_A), \text{std(M}1_B), \text{mean(M}1_A), \text{mean(M}1_B)$ **Question:** I want to compare the two measures tested on the same group of 'people' and use a statistical tool to figure out if the difference in the two measures are significant. I am considering using Welch's t-test – would this be correct?


title: Generalized linear mixed model: what is base dependent variable?

text: I have run Generalized linear mixed model with glmer in lme4. I use R version 3.0.1. My dependent variable is binary (correct or wrong). And this is my results:               > glmer16 <- glmer(result ~ (1|item) + (1|speaker) + vowel + sex + cat + dog + exposure + frequency + v00004 + v00024 + v00034 + v00044, data=data1.frame, family=binomial)     > summary(glmer16)     Generalized linear mixed model fit by maximum likelihood ['glmerMod']      Family: binomial ( logit )     Formula: realisation ~ (1 | item) + (1 | speaker) + vowel + sex + cat +      dog + exposure + frequency + v00004 + v00024 + v00034 +      v00044         Data: data1.frame                 AIC       BIC    logLik  deviance       881.7026  958.6402 -426.8513  853.7026           Random effects:      Groups  Name        Variance Std.Dev.      speaker (Intercept) 7.0291   2.651         item    (Intercept) 0.5084   0.713        Number of obs: 1800, groups: speaker, 50; item, 12          Fixed effects:                 Estimate Std. Error z value Pr(>|z|)         (Intercept) 15.52018    5.33634   2.908  0.00363 **      vowelhigh    0.16750    0.55907   0.300  0.76449         vowellow     0.70981    0.63194   1.123  0.26134         sexmale      1.37080    1.03228   1.328  0.18420         cat         -0.11460    0.09537  -1.202  0.22953         dog         -0.05460    0.03633  -1.503  0.13286         exposure    -0.00404    0.01564  -0.258  0.79613         frequency   -0.01709    0.15594  -0.110  0.91272         v00004      -2.83445    0.66039  -4.292 1.77e-05 ***     v00024       0.29687    0.55868   0.531  0.59515         v00034       0.43899    0.58656   0.748  0.45421         v00044       0.36663    0.65130   0.563  0.57349         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      My questions are: 1) Does v00004 decrease the result of 'wrong' or 'correct'? and 2) what does it mean by significant at the intercept?


title: Confidence interval of two variances

text: I have two normally distributed samples $[x_1, x_2, \ldots, x_n]$ with sample variance $s^2_1$ and another sample $[y_1, y_2,\ldots, y_n]$ with sample variance $s^2_2$. I know how to calculate the confidence interval for each of the variances. However i do not know how to calculate the confidence interval of the sum of the two variances. In other words, can someone help me to calculate the confidence interval for $\sigma^2_1+\sigma^2_2$.


title: What test of significance is best suited to compare concentrations, with only one sample of each population?

text: I have two large reservoirs of water with unknown concentrations of bacteria. I have only one sample of each reservoir: Sample of reservoir A: volume of the sample $=V_1$, number of bacteria in the sample $=N_1$ Sample of reservoir B: volume of the sample $=V_2$, number of bacteria in the sample $=N_2$ $N_1$ or $N_2$ are not necessarily large. What test of significance should I use to estimate the probability that the two reservoirs have different concentrations of bacteria?


title: Empirical justification for the one standard error rule when using cross-validation

text: Are there any empirical studies justifying the use of the one standard error rule in favour of parsimony? Obviously it depends on the data-generation process of the data, but anything which analyses a large corpus of datasets would be a very interesting read.


title: Distribution over a matrix of probabilities

text: I am interested in a problem where I need to infer the distribution of a matrix of probabilities, $\matrix{\Theta}$, where the rows and columns must sum to 1 and every entry lie in the range 0 and 1. Essentially I want to analyze the performance of a set of classifiers where $\theta_{ij}$ is the probability that a pattern belonging to class $i$ will be classified as belonging to class $j$ (for simplicity, I am assuming that $\matrix{\Theta}$ is the same for all classifiers). The problem is that I don't know the true class of any of the observations, so I want to infer $\matrix{\Theta}$ from the labels assigned by each classifier (most of the observations will have been classified by two or more of the classifiers) in a Bayesian manner. Hence I need to be able to define a distribution over $\matrix{\Theta}$. I was thinking about some sort of product of Dirichlet distributions, representing each row and each column. I'd be grateful for any better suggestions.


title: Treating data that is not normally distributed

text: Please allow me to start off by saying that I'm not particularly familiar with statistics and data analysis. I have a data set which has a population size of 266, a mean of 24.8 and a standard deviation of 6.86 on MS Excel 2007. A Normal Quantile plot I drew for this data having undergone a Kolmogorov- Smirnov test indicated that it was not normally distributed (at alpha levels of 0.01, 0.05 and 0.1). How do I find the Confidence Intervals for such a set of data? The graph showing the probability density function for the data is below: ![enter image description here](http://i.stack.imgur.com/qBnIv.png) **Edit 1** A bit more detail as to the nature of the work: I am doing an exercise in analysing a fictional supply chain so as to become familiar with the concepts of statistical analysis and modelling. I want to understand the risk associated with minimum re-order levels decreasing from 2 standard deviations to 1 (suppose the standard deviation for a particular supplier/customer is one week, then the minimum re-order level would be two weeks previously). So I theorize (and please feel free to correct any fallacious logic here) that the risk associated would be the difference in percentages given by the 3-sigma rule (i.e.5% to 32% from 2-sigma to 1). The second part of this exercise requires me to find the optimal minimum re- order time for each product given a set of data for each supplier/customer that is not normally distributed (given some standard financial assumptions). I assume that the logic remains as before, but I am uncertain as to how to treat the data so as to make the 3-sigma rule still apply to it. **Edit 2** Plotting the raw data yields the following graph: ![enter image description here](http://i.stack.imgur.com/Sasne.png)


title: Quantile regression cross validation

text: In a previous question i asked how to calculate the quantiles of my data so that i can do cross validation on a hold out set of data for my quantile regression model. But i think i understood what i need to do incorrectly. Suppose i do three quantile models, with quantiles 0.25, 0.5, 0.75 and got there coefficients and test the coefficients of all three on a hold out data set (data that is not used to calculate coefficients). Do i then just choose the model that give me the best accuracy (actuals vs predicted values) of the three models as my best? I will be using Root mean squared error as a test of accuracy. In my simulations i did, i found that the accuracy is much better for smaller actuals with the 0.25 model and worse for the 0.75 model. And also for the higher actuals the predicted is much closer to the actuals for the 0.75 model than than the 0.25 model. So i thought of combining the models and use the coefficients of 0.25 on smaller values and the coefficients of 0.75 on the larger ones, but i am not sure if this is the correct way.


title: Caret retraining of neural network after finding optimal parameters?

text: I am applying a neural network and logistic regression to a classification problem. In order to evaluate the performance of the two classifiers I'm using 5-fold cross-validation (roughly 800 samples in the test data for each fold). Using caret I am finding the best parameters for my ANN in an inner loop in the cross-validation (roughly 600 samples in test data for each fold). However, I would really like to get a better understanding of how caret handles this. Seeing as neural networks can benefit from retraining in order to avoid local minima I would like to do some retraining. I am not sure whether caret does this after having selected the optimal parameters or simply trains the model once on the larger sample of test data (800 samples) with the found parameters. Can anyone clarify this for me? EDIT: And if it doesn't do any retraining, how can I do this in the best possible way? **Code:**                             library(caret)                   library(nnet)                   tunGrid <-  expand.grid(size  = c(1, 2, 3, 4, 5, 6, 7),                           decay = c(0, 10^(-2), 10^(-3), 10^(-4)))                        fitControl <- trainControl(method = "repeatedcv", number = 5, repeats=3)                   set.seed(825)                   fmla <- as.formula("y_BBR_train ~ .")                     annFit1 <- train(fmla, data = X_BBR_train,                   method = "nnet",                   trControl = fitControl,                   tuneGrid  = tunGrid)      Thank you!


title: Correlation Between Random Variables

text: I recently received the following claim from a work colleague with regard to a loose definition of the correlation between two random variables: “A 60% correlation between two variables approximately implies that in 60% of cases the movement away from the variables' mean position is aligned.” I have not seen correlations defined in such a way, but I was unable to provide an answer as to a) whether this is correct, b) the alternative official definition.


title: Arithmetic mean and Gini index in a grouped distribution

text: In the following table is represented the relative frequency distribution of 100 families related to the class of annual income (in thousands of dollars):               Income:                         10-25   25-40   40-60   60-100  100-200     Relative frequencies:           0,48    0,25    0,15    0,10    0,02      Calculate arithmetic mean and Gini index. I started calculating the middle value of the classes and the cumulative frequency:               Income (middle value):          17,5    32,5    50      80      150     Cumulative frequency:           0,48    0,73    0,88    0,90    1      Problem is, how do I calculate the arithmetic mean? I though it was the summation of the middle values * their relative frequency, which outputs 35,025 (8,4 + 8,125 + 7,5 + 8 + 3). The problem is that the sum of the residuals from the estimated mean isn't zero!               17,5-35,025 = -17,525   +     32,5-35,025 = -2,525    +     50-35,025   = 14,975    +     80-35,025   = 44,975    +     150-35,025  = 114,975   =          154,875 ≠ 0      Is 35,025 still the correct value for arithmetic mean even if the sum of residuals isn't 0? The other question: Gini index. It's a concentration estimator. It's calculated in this way: ![http://upload.wikimedia.org/math/7/5/f/75f99e8670d20b9a7b95bf19dbd31802.png](http://i.stack.imgur.com/9XDPe.png) Where Qi are cumulative frequency and Pi are cumulative frequency in case of equidistribution. Can you help me? :)


title: RandomForestClassifier Parameter Optimization

text: I'm a ML **novice** and I'm wondering if someone can critique what i'm doing (this is a bit open-ended).   * I have a very small corpus of text documents ( **n = 122** ).    * There is a **binary decision** associated with each document.    * I have created a "bag of words" representation of each document and I'm using python's `RandomForestClassifier` to make models to classify the data. I'm tinkering with the parameters in the `RandomForestClassifier` in the following way:   * Run the `RandomForestClassifer` on 200 random subsets of the data (n = 112 for each of the 200 runs) (these #s were chosen arbitrarily).   * Rank the **importance** of each word in my bag of words matrix based on the average importance of each word in the 200 runs. Now I want to see if there is an "optimal" # of feature/words for my data set using the RandomForestClassifier. This is done as follows:   * Generate 500 random forests ( **600 trees per forest** ). Each of the 500 forests uses 112 randomly chosen documents as a training set and the remaining 10 docs as a test set.   * Measure the **average accuracy** of these 500 forests **as a function of # of words/features** used to generate the models. Here is what I see. The optimal average accuracy is around n=80 words/features. ![enter image description here](http://i.stack.imgur.com/2DE6x.png) Questions:   * I'm sure my approach is unorthodox. Is there a **better way** to optimize the RandomForest parameters?   * Is there any "intuitive" **explanation** for why my average accuracy seems to be optimal at around 80 words and then tails off? Is it simply that when n-features gets too large, my forests don't incorporate enough of the good features and so accuracy suffers?   * Any **other parameters** that are worth modifying here?   * Any other classification models worth looking at? Thank you for any thoughts.


title: Exploratory factor analysis and eigenvalues

text: I have reviewed a number of questions regarding this issue and just want to throw this one in the mix - so I apologize for its redundancy. So, I ran an EFA on 60 items. Analysis resulted in 19 components with an eigenvalue of a score greater than 1. The only factors that theoretically make sense and that include more then 3 items have eigenvalues greater than 3 - can I use these first three components in my analysis or do I have to rerun the analysis? And please feel free to let me know how stupid this question is!


title: Highlighting significant results from non-parametric multiple comparisons on boxplots

text: I have box plots of 13 groups that I show in one plot. The groups have unbalanced populations and are not normally distributed. I want to show which pairs are statistically similar (i.e., have kruskal.test p-value < 0.05) by putting a,b,c, etc. on top of boxes that match. Here is a pseudo code to show what I have:               A = c(1, 5, 8, 17, 16, 3, 24, 19, 6)      B = c(2, 16, 5, 7, 4, 7, 3)      C = c(1, 1, 3, 7, 9, 6, 10, 13)      D = c(2, 15, 2, 9, 7)      junk = list(g1=A, g2=B, g3=C, g4=D)      boxplot(junk)       Here is an plot I found that does what I want (except I have 13 groups in one row): ![](https://www.agronomy.org/images/publications/jeq/32/3/1036f3.jpeg)


title: Testing two-tailed p-values using Stouffer's approach

text: About combining multiple test statistics, Wikipedia says > This Z-score (for the overall meta-analysis) is appropriate for one-sided > right-tailed p-values; minor modifications can be made if two-sided or left- > tailed p-values are being analyzed. What are those minor modifications?


title: Random walk on $\mathbb{Z}$

text: Let $X_n=\sum\limits_{k=1}^n Y_k = X_{n-1} + Y_n$ and $X_0=0$, $(Y_n)$ i.i.d with $P(Y_n=1)=p=1-P(Y_n=-1)=1-q, p \in (0,1)$, $(X_n)$ is a random walk on $\mathbb{Z}$. Why is this Markov chain irreducible? How do we know that all states (integers) communicate with each other? How do we show this intuitively? Thank you for your help!