title: Error term interpretation in the Cox PH model

text: I am preparing a presentation on Survival Analysis models, with specific focus on the Cox model. I have been using such model (including its competing risks generalization) for some time now, but I still cannot find an easy way to think about the error term of the model, $\varepsilon$. Suppose I am in the simple bivariate survival analysis case (with $x_1$ and $x_2$). The Cox model functional form for the hazard rate is: $$ h(t| x_1, x_2) = h_0(t) \exp(x_1 \beta_1 + x_2 \beta_2). $$ My question is: if you had to explain to graduate students with a relatively good understanding of Econometrics/Statistics how the error term enters in the functional form, how would you do? I always thought at the baseline hazard $h_0(t)$ as to "containing" the error component of the model, but assuming this is correct, how to formalize this? Can you please suggest a good textbook as reference? (either from a Biostatistics or from an Econometrics perspective). Thank you vey much.


title: unique spline for different groups in a linear model (and no spline at all for one group)

text: I have a problem that has puzzled me for a long time. It involves linear models and spline functions. I need to model "time since diagnosis" when some individuals never had a diagnosis. I use Poisson model where I have split time on age and time since diagnosis. The data set may look like this (with millions upon millions of rows,and with 95 % having no diagnosis ) : SUBJECT AGE DIAGNOSIS TIME_SINCE_DIAGNOSIS EVENT/OUTCOME 1 30 A 0 0 1 31 A 1 0 1 32 A 2 1 2 45 B 0 0 2 46 B 1 0 103 22 NON "0 or -" 0 103 23 NON "0 or -" 0 . . . You get the idea .. Now I want to have a model that looks something like this: log ( OUTCOME ) = Spline( AGE ) + DIAGNOSIS * Spline( TIME_SINCE_DIAGNOSIS ) + offset ( logrisktime ) and then plotting the predicted risk for the different diagnoses of A and B as a function of " time since diagnosis " given that the diagnosis was at age = 30, and in the same plot for diagnosis = NON but then instead as a function of age ( with year since 30 on x-axis). I like to have both those with a diagnosis and those with diagnosis=NON in the same model to be able to do contrasts/linear combinations of effects. Model and plot is done, but my problem is that the model above also creates a spline over the " time since diagnosis " for those who have no diagnosis (i.e. DIAGNOSIS = NON ). So, I would like to have a model of this type which I created indicator variables for the different diagnoses : R: factor(A) factor(B) fator(NON) SAS: class A B NON; log ( OUTCOME ) = Spline ( AGE ) + NON + A*Spline( TIME_SINCE_DIAGNOSIS ) + B*Spline ( TIME_SINCE_DIAGNOSIS ) + offset ( logrisktime ) Are you with me ? So I want to create a unique spline of the variable "time since diagnosis" for each level of the variable "diagnosis" (which is done with diagnosis*spline(time_since_diagnosis) ), but with the exception that one level (those who had no diagnosis) does not to have any spline assigned. Idea ? Thanks, Peter


title: A random censored regression problem?

text: Assume I have random variables $X$ and $Y$. $X$ is observable; $Y$ is not. We also know that, if $X \leq Y$, $Z = f(X,Y)$ where $Z$ is a third random variable and $f(\cdot)$ is a known function; otherwise $Z$ is missing. $Z$ is also observable (missing or a real number). Assume I know the distribution of $X$ and I can estimate the distribution of $Z$ from data. I wonder if I can estimate the distribution of $Y$. My first question is: what is this problem? It is similar to but different from a textbook censored regression problem. Is it a random censored regression problem? My next question is, is the problem identifiable? I realize that, for any given set of random samples $X_1, X_2, X_3, ...$, there may exist multiple sequences $Y_1, Y_2, Y_3, ...$ that yield the same output $Z_1, Z_2, Z_3,...$. So I guess it will be difficult to find _the best_ solution unless some constraints are imposed. If my guess is right, what are typical constraints people apply? If we can appropriate apply some constraints to the problem so that we can try to find optimal parameters that minimizes some pre-defined error metric. Hopefully we can "estimate" the distribution of $Y$. Your suggestions, books, and papers are more than welcome. Any idea?


title: Efficient online linear regression

text: I'm analysing some data where I would like to perform ordinary linear regression, however this is not possible as I am dealing with an on-line setting with a continuous stream of input data (which will quickly get too large for memory) and need to update parameter estimates while this is being consumed. i.e. I cannot just load it all into memory and perform linear regression on the entire data set. I'm assuming a simple linear multivariate regression model, i.e. y = Ax + b + e What's the best algorithm for creating a continuously updating estimate of the linear regression parameters A and b? Ideally:   * I'd like an algorithm that is most O(N*M) space and time complexity per update, where N is the dimensionality of the independent variable (x) and M is the dimensionality of the dependent variable (y).   * I'd like to be able to specify some parameter to determine how much the parameters are updated by each new sample, e.g. 0.000001 would mean that the next sample would provide one millionth of the parameter estimate. This would give some kind of exponential decay for the effect of samples in the distant past.


title: Can you still use the least significant difference test even if the F test is not significant?

text: Can you still used the least significant difference test even if the F test is not significant?


title: Sample size for several different populations

text: I would like to find out how quickly etymological (history of words) knowledge becomes obsolete. I thought I could compare several editions of one dictionary published over a long period of time. My problem is that the size of this dictionary has grown during this period from four to twelve thousand words, and I do not know how to choose the sample size in this situation.


title: Support Vector Machine for Longitude-Latitude data

text: I am wondering if someone could point me to the direction of support vector machines being used for longitude latitude. It seems logical that the possible complexity in SVM would be great for modelling variables dependent of location. In this case, I am trying to model depression dependent on US zip codes. My idea is to convert these zip codes to longitude/latitude, and then use an SVM. Has this approach been used before, would someone recommend a different approach?


title: Proper name for this distribution type

text: I need help in defining the type of distribution used in software so that I could use some standard distribution library for this purpose. I apologize for not using proper terms. It takes "center" value and an alpha parameter (0 to 1) and produces N values within a given range, so that values are "clustered" around this "center" value. The less alpha the less they are "centered". For example, if I feed a range of 1-5 with center at 4, alpha=0.5 and ask for 20 values the result will look something like this:               1:1 times     2:3 times     3:4 times     4:7 times     5:5 times      You get the idea. Thanks! P.S. In case it might help, I provided some extracts from code The code for probability density function:               private function assessProbabilityDensityFunction(){             $probabilityDensityFunction=array();             for ($i = 0; $i < $this->valueRange->getNumberOfValues(); $i++)                 $probabilityDensityFunction[$i] = $this->probabilityFunction() * $this->probabilityCenteredCoefficient($i);             return $probabilityDensityFunction;         }          private function probabilityFunction(){             return (1 - $this->settings->getAlpha()) / (1 - pow($this->settings->getAlpha(), $this->valueRange->getNumberOfValues() - $this->CentralValueSerialNumber() + 1) + $this->settings->getAlpha() - pow($this->settings->getAlpha(), $this->CentralValueSerialNumber()));         }     private double centerBasedCoefficient(int valueIndex) {                 return Math.pow(parameters.alpha, distanceFromCenter(valueIndex)+1);             }          private function distanceFromCentralValue($serialNumber){             return abs($this->CentralValueSerialNumber() - $serialNumber - 1);         }      And here's the cumulative function code               private function assessCumulativeDistributionFunction(){         $cumulativeDistributionFunction = array();         $cumulativeDistributionFunction[0] = $this->probabilityFunction->getValue(0);         for ($i = 1; $i < $this->getValueRange()->getNumberOfValues(); $i++)             $cumulativeDistributionFunction[$i] = $cumulativeDistributionFunction[$i - 1] + $this->probabilityFunction->getValue($i);         return $cumulativeDistributionFunction;     }      And here's how we get distributed values               public function getValue(){        return  $this->multinomialDistribution()+$this->cumulativeFunction->getValueRange()->getFromValue();     }          private function multinomialDistribution(){        $rnd= lcg_value(); // random (0,1)        for ($i=0; $i<$this->cumulativeFunction->getValueRange()->getNumberOfValues(); $i++)          if ($this->cumulativeFunction->getValue($i) > $rnd) return $i;        throw new Exception("Should always return a value");     }


title: Significance of Difference between Two Extremes (Maxima/Minima)

text: I am comparing two different populations (with different sizes) of models, for each individual I have a score. I'd like to compare not only the difference of the two means but the maximum value. I've learnt in the past to perform significance test when comparing averages between groups but this is the first time that I'm trying to do the same with an extreme value (maximum in this case). To summarise, my question is: what is the right way to know whether the difference $max(P_1) - max(P_2)$ is significant?


title: Why is compound symmetry sufficient to establish sphericity?

text: At this website it is stated that "If compound symmetry is met then sphericity is also met." Is it impossible to construct a data set that has perfect compound symmetry but not perfect sphericity? Why/why not? It's not intuitively obvious to me because at face value compound symmetry (variances across conditions approximately equal, and covariances between pairs of conditions also approximately equal) and sphericity (approximately equal variances of the differences between data from the same participant at different times) seem like quite different things. I'm thinking about the context of repeated-measures ANOVA.


title: DeLong's Test to compare ROC curves

text: I am trying to compare the performance of a short and long form of a questionnaire. The long form has 25 items and the short form has 8 of the 25 items. Can DeLong's Test be used to compare the ROC curves? If not which test is most suitable for this analysis?


title: What is the best way to reduce false negative percent in the model?

text: My logistic regression model seems to identify successes very well (about 85% - 94%), but fails to identify the failures (only identifying 18% - 32% correctly). I have thought of weighting the success and failures differently, but I am not sure of how to do this R. I have tried:               Fr <- ifelse (y==0, 2, 1)     model <- glm(y~factor(x1) + factor(x2), weights=Fr,                                              family = "binomial",                                              data = data)      but did not get any improvement this way.


title: Confidence interval for hidden markov model (MATLAB preferred)

text: I'm trying to uncover the transition parameters of data of a hidden Markov Model using MATLAB. Using the built in hmmtrain function, I can estimate the parameters quite well (I already know what they are supposed to be as it's simulated data for now), but they are always off by a little. So what I want to do is construct a confidence interval for these parameters. Now, I've searched the web for this quite a bit and it seems as if this is not easy to do. The 'easiest' solution appears to be bootstrapping, but I personally don't know how to go about and implement something like that in MATLAB, for a model as complex as a HMM. Does anyone have any references I could look at that could help me write a script for this, or does anyone know of a package out there that already has this built in? The main issue I have is that I wouldn't know how a bootstrap method would tackle the problem of HMM. Isn't the whole idea of HMM that you get a specific sequence from which you learn about the parameters? Wouldn't cutting it up for the purpose of bootstrapping automatically create some errors?


title: Help me come up with a witty title for my paper

text: My honours thesis is due in soon and I still haven't thought of a great title. I want something that's clever (perhaps using a pun or two), but also gives insight into what I'm doing. Preferably in the form of "Thesis: A Brilliant Thesis" Everyone loves colons. So what is my thesis about? It is a statistical dissertation that is exploring the relationships between a number of economic variables and the price of Gold. What I have found is that the last 12 years has destroyed all the relationships between Gold and other factors. It is also worth noting that Gold has been booming in the last ten years. I was thinking something along the lines of: The Gold Rush: Modelling the Price of Gold but it's not great.


title: Statistical Significance Test on Multiple Algorithms and Single Dataset

text: I have N binary classifier algoritmhs trained on same single train set and tested on single test set. 1) Which statistical test is appropriate for calculating their statistical difference? 2) Using McNemmar's test with pairwise manner is a valid option?   * I know ANOVA, which works on K train/test sets whereas i have single train/test set.


title: Do all interactions terms need their individual terms in regression model?

text: I am actually reviewing a manuscript where the authors compare 5-6 logit regression models with AIC. However, some of the models have interaction terms without including the individual covariate terms. Does it ever make sense to do this? For example (not specific to logit models):               M1: Y = X1 + X2 + X1*X2     M2: Y = X1 + X2     M3: Y = X1 + X1*X2 (missing X2)     M4: Y = X2 + X1*X2 (missing X1)     M5: Y = X1*X2 (missing X1 & X2)      I've always been under the impression that if you have the interaction term X1*X2 you also need X1 + X2. Therefore, models 1 and 2 would be fine but models 3-5 would be problematic (even if AIC is lower). Is this correct? Is it a rule or more of a guideline? Does anyone have a good reference that explains the reasoning behind this? I just want to make sure I don't miscommunicate anything important in the review. Thanks for any thoughts, Dan


title: Can we skip the lower order terms in interactions?

text: This question is about three-way interaction and the possibility of applying without second lower terms with keeping the main variables in the equation not like the other questions. In fact the other answers suggest there is possibility of applying . I am not here to find the best solution because I know it and I already included in my question, but to know whether is it possible regardless if it is preferable or not. thank you and please open my question for discussion The widely known regression equation for assessing the three-way interaction is $$ Y= B_1 X+B_2 Z+B_3 W +B_4XZ+B_5XW+B_6ZW+B_7XZW+B_0 $$ All lower order terms is included in the regression equation for the B7 coefficient to represent the effect of the three-way interaction on Y. Is there possible way to skip the lower order terms and include only the higher term? as in: $$ Y= B_1 X+B_2 Z+B_3 W +B_4XZW+B_0 $$ And how many observations do I need to perform such equation if X & Z are continuous variables and W is dummy variable ? I will be thankful if anyone can provide me with any suggestions


title: Significance test for multiclass classifier

text: In a multiclass classification problem, I want to measure the significance of my classifier against the null hypothesis (in this case, chance level). In this paper, in section 3.4, for a binary classifier, it is proposed as a binomial test. The resultant p-value is a measure. For a multiclass version, how can I make such a measure?   * Is a multinomial test appropriate?    * Or should I apply a binomial test in one-vs.-all fashion and evaluate p-values for each class? Here a chi-squared test is applied on a confusion matrix. Can this evaluation method be a measure for me? Thanks in advance.


title: ar.ols() or arima() for modeling time series

text: When trying to fit a AR(p) model to a time series in R, it seems that both ar.ols() and arima() will work. Is there some consideration of when to use which then? ar.ols() seems to use least square estimation, while arima() uses maximum likelihood or minimize conditional sum-of-squares (default method). Is minimizing conditional sum-of-squares equivalent to conditional MLE for Gaussian stationary process? So is arima() essentially MLE or conditional MLE? Thanks!


title: Is McNemar's Test appropriate for evaluating pre-post differences in intervention studies with a control arm?

text: I recently came across a paper that uses the McNemar test to evaluate the effectiveness of an intervention to improve adherence to treatment. The study used a pre/post test design whereby the adherence behaviour for both the control and intervention group was measured using the same instrument and expressed as continuous data. It then compares the relative improvement (post / pre) between the two groups. I've never seen the McNemar test used for this before, but am not sufficiently familiar with statistics. Is it really an appropriate method? Wouldn't a (paired) t-test or something like that make more sense?


title: Calculating AUPR in R

text: It is easy to find a package calculating area under ROC, but is there a package that calculates the area under precision-recall curve?


title: Visualizing a spline basis

text: Textbooks typically have nice example plots of of the basis for uniform splines when they're explaining the topic. Something like a row of little triangles for a linear spline, or a row of little humps for a cubic spline. This is a typical example: http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_introcom_a0000000525.htm I'm wondering if there is an easy way to generate a plot of the spline basis using standard R functions (like bs or ns). I guess there's some simple piece of matrix arithmetic combined with a trivial R program which will spit out a pretty plots of a spline basis in an elegant way. I just can't think of it!


title: Bivariate Poisson Regression Given Scoring Rates and Probabilities of each Score and not the outcomes of matches

text: I'm trying to come up with a bivariate poisson model for predicting the outcomes of soccer games. I have goal scoring rates for the home and away teams, and the probabilities of each score (i.e. 1-0, 2-1, 3-3) being scored (given by bookmakers) for thousands of games. I am assuming that the probabilities that the bookmakers come up with for each possible score and match my scoring rates for each team are correct. I have read several articles on this topic, the most helpful was Bivariate Poisson and Diagonal Inflated Bivariate Poisson Regression Models in R. They all use the **results** of past soccer games in order to come up with scoring rates for each team, and then perform a regression to fit the rates to the outcomes. I want to use the scoring rates I already have, and the probability distributions given by bookmakers, in order to construct a bivariate poisson model. Everyone seems to be using the results of past games to create a model because they don't have the scoring rates of each team. I do have them, so I believe it is more efficient to come up with a model this way. The end goal is to have a model that I can input the scoring rates for team A and B and into the model and output a probability distribution of each score occurring (i.e. 0-0, 1-3, 2-1, ...). How should I go about doing this? The bivpois package the creators of the article made use of is incredibly useful if you're using the results of past games. Can you suggest another package for R that I can use? Thanks for your help.


title: How to examine predictive value of variables for an outcome that can only be observed as a group?

text: Suppose I have a bunch of people who did a few tests (let's call it tests A). Now I'll do a second test (test B) in order to categorize these people into two categories (that have a pre-defined "true" answer). _Note that score from test A can only predict the signal-to-noise ratio of test B. I.e you can't use score from test A to do the categorization itself._ My aim is to select a group of people that are as easily categorized into correct categories as possible (I.e actual score of test B does not matter), using the score of test A. Now, I want to examine how test A (which consists of multiple variables A1, A2...) predicts the final outcome. How could I do that? I thought about creating bunch of subgroups of a fixed sample size randomly and measure the mean of test score A and the percentage I got the categorization right for these subgroups, and then apply regression analysis on that. Would this be a valid approach? If not, what is the correct way of doing this?


title: Understanding the Meta-analysis Model output in layman terms

text: I am conducting a **meta-analysis** from a large number of studies. In each study are **compared weights of two groups** (fishes with and without internal parasite). I am interested if the **weight can explain the presence/absence of a parasite**. From forest plot it seems to be clear that in each continent (and for world as a whole) there is a **clear preference towards bigger fishes**. I have studies from all over the world, that means **the term "bigger" fish is a relative**. The bigger fish in Africa could be very small in comparison with smallest fish in Australia. **Because of this fact I have used random-effect model**. It should be a good choice according to some textbooks. **Is random-effect model appropriate for my situation?** Could you please help me with **interpretation of model output**? There are some things I do not understand: P-value from model results is highly significant. So it is highly probable that the weight has something to do with parasite infestation. However, the Test for Heterogeneity is also significant. Does it means that **my model do not meet the assumption for normal distribution of residuals?** Similarly like a in the case of some simple regression? **And what about tau^2, tau, I^2 and H^2? Is any of them similar to R-squared?** Please give me some guidance (in layman terms if possible).               library(metafor)     mod_weight <- rma(yi, vi, data = dat_weight); summary(mod_weight)          Random-Effects Model (k = 62; tau^2 estimator: REML)            logLik  deviance       AIC       BIC      AICc       -65.0051  130.0103  134.0103  138.2320  134.2172            tau^2 (estimated amount of total heterogeneity): 0.3618 (SE = 0.0808)     tau (square root of estimated tau^2 value):      0.6015     I^2 (total heterogeneity / total variability):   86.67%     H^2 (total variability / sampling variability):  7.50          Test for Heterogeneity:      Q(df = 61) = 395.7163, p-val < .0001          Model Results:          estimate       se     zval     pval    ci.lb    ci.ub                 0.8007   0.0853   9.3830   <.0001   0.6334   0.9679      ***       **_P.S. the fish-parasite research is a made-up story :)_**


title: Are there any *very* user friendly (preferably video) explanations of Baum-Welch?

text: I'd like to understand the Baum-Welch algorithm. I liked this video on the Forward Backward algorithm: http://www.youtube.com/watch?v=7zDARfKVm7s&feature=related I'm having trouble coming up with good resources for Baum-Welch. Any ideas? Thanks


title: Criteria for model selection for nested models fitted on nested data?

text: I consider a first model where the 6 observables (concentrations of metabolites) are fitted on the data set (the experimental measure of these 6 concentrations). I also have a second model, that is a subset of the first one and models only 4 observables. This second model is fitted on only a part (4 out of 6 concentrations) of the data set used for the first model. I cannot use the AIC to compare the two models because they are not fitted on exactly the same data. Is these an extension of the AIC or another information criterion that could help me in my case of "nested" models and "nested" data? (If you have a reference, would be great!)


title: Significance of time course difference over interval

text: I am dealing with some data which plotted looks like: ![enter image description here](http://i.stack.imgur.com/k118s.png) I would like to add a level of significance to my claim that the values of the two time courses are different over the latter half of the experiment (data is sampled at 60 Hz, standard errors of the mean are shaded). I have tried two approaches:   * Compare per-interval means over this interval - got $p = 2.8 \times 10^{-1}$    * Compare all values in the interval - got $p = 5.6 \times 10^{-53}$ Now, the first approach tells me that my claim is pretty much not significantly supported by my data. The second one tells me that it is - but the p value I get is ridiculously low, though the approach seems sensible and I would not call it inflated. Could you please tell me what you think is better or how you would advise me to go about this? I have tried fitting a mixed model to this, but the nonlinear time course makes this a bit difficult, and it will surely take me some time to test multiple models with multiple nonlinear components. I plan to do this eventually, but I am currently looking for a somewhat faster (even if less accurate) approach.


title: How to bootstrap confidence of a weighted measure?

text: I have a large data set of 20M observations with 20+ measures for each observation (people). Additionally I have a set of weights that reflect the racial population where each person lives. For example, from Census data this could be the percentage of population in the zipcode, or tract that is African American, White, Asian etc. These percentages will add to 100%. We lack direct measure of race on individuals. This information is not available. The profile data will then look like this:               ID,White,Asian,AfricanAmerican,...     1234    70% 20% 5%  5%  0%     2345    30% 50% 10% 5%  5%     3456    5%  90% 5%  0%  0%     4567    25% 25% 25% 15% 10%      The 20 measures are a combination of financial data and true/false flags. These finance data are heavily skewed, and also have a lot of zeroes. We will compute the weighted means of all the measures, e.g. average value of African Americans. My questions center around computing variance measures for these means. I think the most straightforward way to compute variance estimates is to bootstrap the IDs, preserving the profile of each ID and compute $N$ weighted means. Alternatively, we could randomly pick one race per person $N$ times according to the profile and compute $N$ weighted means. We could also do both, re- sample the IDs and pick 1 race in each iteration. I cant figure out what the "one race" simulation results actually measure. Because each person is assigned one race, it can't be used for a confidence interval of the original weighted means. _**Is there a name for this second method of simulating the random variable?_** Any references for using it?


title: Two questions on indexing and regression

text: 1. Is it OK to have causals indexed on different base years? Amongst many references, I think the answer is “no”. Is that true? _p.360 Wooldridge [Introductory Econometrics: A Modern Approach, 5th ed]._ says you should have a common base year. Just google… “common base year” regression and you will see loads of support.    2. Is it OK to have the Y variable as an index based on 100? For example 100 ,98, 103, 101, etc.


title: Standard error for the sum of weighted means

text: I am looking at a two-stage process of diagnosis. Stage 1 is a screen that identifies people with the diagnosis with some error, Stage 2 is a detailed examination (a gold standard) In Stage 1, N people were screened, and S were screened positively. When the S people screened positively were examined in detail it was discovered that G actually had the disease A 10% random sample of the those screened negative were also examined in detail, and it was found that J actually had the disease (the false negatives from the first stage) The prevalence of the disease is the sum of the true positives and the false negatives from the first stage, or (G/S)(S/N) + (J/(N-S))((N-S)/S)) = (G+J)/N My question is, what is the standard error of this prevalence estimate?


title: Using monte carlo integration with bisection method to find a true value in R

text: I try to (re)produce the results of Austin, PC. (2010).. The question is to obtain an approximated $\beta$, which leads to a specific risk difference(difference of proportions) $p_0 - p_1 = 0.02$, where $p_0$ & $p_1$ are proportions of subjects in the control and trt groups that have the response $Y = 1$, respectively. There are 4 covariates, 1 treatment indicator, and 1 outcome variable:$x_1, x_2, x_3,$ and $x_4$, which               set.seed(123)     N = 10000     x1=rnorm(N);x2=rnorm(N);     x3=rbinom(N,1,0.5);x4=rbinom(N,1,0.5)      where $N$ represents the number of subjects in each dataset and there are total $s = 1000$ simulated datasets. For each subject, we compute the $Pr(Y = 1)$ if a subject is in the trt $(t=1)$ or control $(t=0)$is $p_{i,1} = \frac{1}{1+e^{-(\alpha_0+\alpha_1x_{1i}+\alpha_2x_{2i}+\alpha_3x_{3i} + \alpha_4x_{4i} + \beta_{(k)})}}$, or $p_{i,0} = \frac{1}{1+e^{-(\alpha_0+\alpha_1x_{1i}+\alpha_2x_{2i}+\alpha_3x_{3i} + \alpha_4x_{4i})}}$ where the values of $\alpha_0=-1, \alpha_1=\alpha_2=\alpha_3=\alpha_4=1$ are fixed, but $\beta_{(k)}$ is approximated by monte carlo integrations with $k$th iterative process by using bisection method : $\bar{p}_1 = \quad \idotsint_{x_1} \frac{1}{1+e^{-(\alpha_0+\alpha_1x_{1i}+\alpha_2x_{2i}+\alpha_3x_{3i} + \alpha_4x_{4i} + \beta_{(k)})}} f(x_1)...f(x_4)dx_1...dx_4$ and $\bar{p}_0 = \quad \idotsint_{x_1} \frac{1}{1+e^{-(\alpha_0+\alpha_1x_{1i}+\alpha_2x_{2i}+\alpha_3x_{3i} + \alpha_4x_{4i})}} f(x_1)...f(x_4)dx_1...dx_4$. $\bar{p}_0^{mc(k)} = \frac{1}{10000}\sum_{i=1}^{10000}p_{i,0}$ and $\bar{p}_1^{mc(k)} = \frac{1}{10000}\sum_{i=1}^{10000}p_{i,1}$, where $mc$ represents monte carlo. The empirical marginal risk difference of $s^{th}$ dataset is $\gamma_{(s)}^{(k)} = \bar{p}_0^{mc(k)} - \bar{p}_1^{mc(k)}$. We repeated over 1000 datasets, the mean empirical marginal risk difference after $k^{th}$ iterative process (bisection method) is $\gamma^{(k)} = \frac{1}{1000}\sum_{s=1}^{1000}\gamma_{(n)}^{(k)} =0.02$ How can I write a R code to obtain approximated $\beta$ for $p_0 - p_1 = \gamma^{(k)}=0.02$? many thanks in advance.


title: Probing of 3-way-interaction with one continuous and two categorical variables

text: I have a model for multiple moderated regression with `payment` (continuous) as the DV, one continuous predictor (neuroticism), and two binary categorical between-subject-factors: stress vs. no stress and self-payment vs. other- payment. Essentially, I want to know if the increase in the difference between self-and other payment depending on neuroticism is larger under stress than without stress. So, I tested the 3-way interaction for significance and found a significant effect but as I understand 3-way interactions, a significance can either support my alternative hypothesis or the hypothesis that other two-way interactions vary. How do I test exactly the above mentioned hypothesis and how do I implement this in SPSS? I am also unsure as to which coding system I should use in this case and why in terms of interpretability.


title: Probability of an object of belonging to a group of objects

text: Say I have a data set composed of $N$ objects. Each object has a given number of measured values attached to it (three in this case): $x_1[a_1, b_1, c_1], x_2[a_2, b_2, c_2], ..., x_N[a_N, b_N, c_N]$ meaning I have _measured_ the properties $[a_i, b_i, c_i]$ for **each** $x$ object. The _measurement space_ is thus the space determined by the variables $[a, b, c]$, where each $x$ object is represented by a point in it. More graphically: I have $N$ objects scattered in a 3D space. What I need is a way to determine the probability (or likelihood?, is there a difference?) of a **new** object $y[a_y, b_y, c_y]$ of belonging to this cloud of objects This probability calculated for any $x$ object will of course be very close to 1. Is this feasible? * * * ## Add 1 To address AdamO's question: the object $y$ belongs to a set composed of $M$ _mixed_ objects (with $M$ > $N$). This means that some objects in this set will have a high probability of belonging to the first data set ($N$) and others will have a lower probability. I'm actually interested in these _low probability_ objects. I can also come up with up to 3 more data sets of $N1$, $N2$, and $N3$ objects, all of them having the same _global properties_ as those in the $N$ data set. Ie: an object in $M$ that has a low probability in $N$ will also have low probabilities when compared with $N1$, $N2$ and $N3$ (and viceversa: objects in $M$ of high probabilities of belonging to $N$ will also display high probabilities of belonging to $N1$, $N2$ and $N3$). * * * ## Add 2 According to the answer given in the question Interpretation/use of kernel density I can **not** derive the probability of a new object belonging to the set that generated the $kde/pdf$ (assuming I would even be able to resolve the equation for a non-unimodal $pdf$) because I have to make the a priori assumption that that new object was generated by the same process that generated the data set from which I obtained the $kde$. Could someone confirm this please?


title: Multi armed bandit for general reward distribution

text: I'm working on a multi-armed bandit problem where we do not have any information about the reward distribution. I have found many papers which guarantee regret bounds for a distribution with known bound, and for general distributions with support in [0,1]. I would like to find out if there is a way to perform well in an environment where the reward distribution does not have any guarantees about its support. I'm trying to compute a nonparametric tolerance limit and using that number to scale the reward distribution so I can use the algorithm 2 specified on this paper (http://jmlr.org/proceedings/papers/v23/agrawal12/agrawal12.pdf). Does anyone think this approach will work? If not, can anyone point me to the right spot? Thanks a bunch!


title: Mahalanobis distance on non-normal data

text: Mahalanobis distance, when used for classification purposes, typically assumes a multivariate normal distribution, and the distances from the centroid should then follow a $\chi^2$ distribution (with $d$ degrees of freedom equal to the number of dimensions/features). We can calculate the probability that a new data point belongs to the set using its Mahalanobis distance. I have data sets that do not follow a multivariate normal distribution ($d \approx 1000$). In theory, each feature should follow a Poisson distribution, and empirically this seems to be the case for many ($\approx 200$) features, and those that do not are in the noise and can be removed from the analysis. How can I classify new points on this data? I guess there are two components:   1. What is an appropriate "Mahalanobis distance" formula on this data (i.e. multivariate Poisson distribution)? Is there a generalization of the distance to other distributions?   2. Whether I use the normal Mahalanobis distance or another formulation, what should the distribution of _these_ distances be? Is there a different way to do the hypothesis test? **Alternatively...** The number of known data points $n$ in each class varies widely, from $n=1$ (too few; I'll determine a minimum empirically) to around $n=6000$. The Mahalanobis distance scales with $n$, so distances from one model/class to the next cannot be directly compared. When the data is distributed normally, the chi-squared test provides a way to compare distances from different models (in addition to providing critical values or probabilities). If there is another way to directly compare the "Mahalanobis-like" distances, even if it does _not_ provide probabilities, I could work with that.


title: Accuracy vs F1 Measure in Multilabel Classification

text: I'd like to evaluate a multilabel classification algorithms and I was thinking of using both Accuracy and F1-Measure with:               Accuracy = #Intersection(suggestions,correct_labels) / #Union(suggestions,correct_labels)     F1 Measure = 2 * (P * R) / (P + R)      The above is for a single example and would then be averaged according to the test size. Does it make sense to consider both of the measures above and why? Are there any circumstances under which one outperforms the other?


title: Why is boosting effective?

text: Boosting takes a bunch of weak learners and creates a strong learner. But why is it so difficult to create a strong learner right from the beginning without using boosting techniques? And therefore remove the need for boosting.


title: Why are there not obvious improvements over Tukey's method?

text: Is there a book that explains why there aren't better standard techniques than Tukey and ANOVA, for example? For comparison consider for example I read about the null hypothesis one- sample $t$ test and didn't even bother considering that there could be better tests. But that is probably just a biased belief that there isn't anything more mathematically sophisticated that can be done with the Gaussian to improve on the one-sample $t$ test. However, Tukey and ANOVA are mathematically more complicated and it didn't feel obvious to me why they should even be considered in the first place. For example in a previous question I asked about why Tukey's method is needed over all pairwise two-sample $t$-tests. The answer I got for that question is that all pairwise two-sample $t$-tests suffers from false positives. But its not intuitive to me how Tukey's method is the best way to evade that problem. How does one intuitively see that generically Tukey and ANOVA are reasonably very good techniques among all possibilities?


title: How to get confidence interval on population r-square change

text: For sake of a simple example assume that there are two linear regression models   * Model 1 has three predictors, `x1a`, `x2b`, and `x2c`   * Model 2 has three predictors from model 1 and two additional predictors `x2a` and `x2b` There is a population regression equation where the population variance explained is $\rho^2_{(1)}$ for Model 1 and $\rho^2_{(2)}$ for Model 2. The incremental variance explained by Model 2 in the population is $\Delta\rho^2 = \rho^2_{(2)} - \rho^2_{(1)}$ I am interested in obtaining standard errors and confidence intervals for an estimator of $\Delta\rho^2$. While the example involves 3 and 2 predictors respectively, my research interest concerns a wide range of different numbers of predictors (e.g., 5 and 30). My first thought was to use $\Delta r^2_{adj} = r^2_{adj(2)} - r^2_{adj(1)}$ as an estimator and bootstrap it, but I wasn't sure whether this would be appropriate. ### Questions   * **Is $\Delta r^2_{adj}$ a reasonable estimator of $\Delta \rho^2$?**   * **How can a confidence interval be obtained for the population r-square change (i.e., $\Delta\rho^2$)?**   * **Would bootstrapping $\Delta\rho^2$ be appropriate for confidence interval calculation?** Any references to simulations or the published literature would also be most welcome. ### Example code If it helps, I created a little simulation dataset in R which could be used to demonstrate an answer:               n <- 100     x <- data.frame(matrix(rnorm(n *5), ncol=5))     names(x) <- c('x1a', 'x1b', 'x1c', 'x2a', 'x2b')     beta <- c(1,2,3,1,2)     model2_rho_square <- .7     error_rho_square <- 1 - model2_rho_square     error_sd <- sqrt(error_rho_square / model2_rho_square* sum(beta^2))     model1_rho_square <- sum(beta[1:3]^2) / (sum(beta^2) + error_sd^2)     delta_rho_square <- model2_rho_square - model1_rho_square          x$y <- rnorm(n, beta[1] * x$x1a + beta[2] * x$x1b + beta[3] * x$x1c +                    beta[4] * x$x2a + beta[5] * x$x2b, error_sd)          c(delta_rho_square, model1_rho_square, model2_rho_square)     summary(lm(y~., data=x))$adj.r.square -              summary(lm(y~x1a + x1b + x1c, data=x))$adj.r.square      ### Reason for concern with bootstrap I ran a bootstrap on some data with around 300 cases, and 5 predictors in the simple model and 30 predictors in the full model. While the sample estimate using adjusted r-square difference was `0.116`, the boostrapped confidence interval were mostly larger CI95%(0.095 to 0.214) and the mean of the bootstraps was nowhere near the sample estimate. Rather the mean of the boostrapped samples appeared to be centred on the sample estimate of the difference between r-squares in the sample. This is despite the fact that I was using the sample adjusted r-squares to estimate the difference. Interestingly, I tried an alternative way of computing $\Delta\rho^2$ as   1. calculate sample r-square change   2. adjust the sample r-square change using the standard adjusted r-square formula When applied to the sample data this reduced the estimate of $\Delta \rho^2$ to `.082` but the confidence intervals seemed appropriate for the the method I mentioned first, CI95% (.062, .179) with mean of .118. Broadly, I'm concerned that bootstrapping assumes that the sample is the population, and therefore estimates that reduce for overfitting may not perform appropriately.


title: Frechet Mean of Multiple Curves

text: I'm working on determining the Frechet mean of several curves. I am using this equation: $$ \sum_{n=1} (P_a - P_b)^2 $$ but I'm not entirely sure how I would do this across multiple curves. Also, the answer should be, unless I am completely mistaken, should be a curve but I am getting a single number instead. Can someone help me out?


title: Wilcoxon signed-rank test. Is it necessary?

text: I have a sample with millions of points. Every point $x$ has two values associated, let's say $A_{x}$ and $B_{x}$. I calculated the mean values of A and B in the sample: $\overline{m}_{A} = 0.19 \quad \text{and} \quad \overline{m}_{B} = 0.21$ The distributions of the values are different and not normal. And the standard deviations of the means are negligible. Can I say that the means are statistically significantly different? Should I calculate the Wilcoxon signed-rank test to prove it?


title: How to choose an arrival rate distribution?

text: I have a system where arrival rate is assumed to be Poisson but there is a condition on accepting this arrival which is if $Z$ (uniform dist. random variable) $\le$ some value $X$ and $X$ can vary between the values $[0.5, 0.25, 0.125]$. Actually X takes value of 0.5 first and if Z>X for two times --> X became 0.25 and if Z>X for another two times --> X became 0.125 and continue with 0.125 till Z<=X occurred. How can I choose an arrival rate distribution under these conditions?


title: ANOVA-like test for bounded variables (percentages)

text: I'm interested in determining whether two or more groups of data share the same mean, and it seems like the ANOVA framework is a good way to approach this. However, ANOVA assumes residuals are normally distributed, while each of my data is a number between 0 and 100 (a percentage). Because normal distributions have support on all the real line, my data cannot adhere to the normality assumption. My question: **Is there an ANOVA-like test which is suited to predicting percentages** (i.e., a variables that is bounded)?


title: Obtain a Neyman-Pearson Test of H0: p1=p2=2/5 against the alternative p1=1/2 and p2=1/5 when the type 1 error is 0.05

text: I kinda have some trouble trying to start with this question. ![a busy cat](http://i.imgur.com/JKeSBhN.png) I know how to derive a Neyman-Pearson Test and its critical region. My main concern is deriving the Likelihood function of the sample given its null hypothesis and its alternative hypothesis. I'd like to solve this on my own but I'm not sure how to start it.


title: How compare AIC (akaike information criterior)?

text: I have two function. First function have a positive value for AIC and the second function have a negative value fir AIC. How to make comparisons? algebraic or absolute values?? Thanks :)


title: Should weights of neural network without hidden layer match logistic regression?

text: Should the weights of a neural network without hidden layer and a logistic activation function be the same as the parameters of a logistic regression? Mine are not the same?               nnallnohidden=nnet(         PartialPrepayzo~FIXPER+MEDSAL2+DREL+LEEFTIJD+HH2CRED+LTV_curr+         rate1Y+rate5Y+CIremFIRP+URB+WELSTAN2+OutNot+mover+SavRate+CRate,         data=test,         size=0,         skip=T)          log <- glm(PartialPrepayzo~FIXPER+MEDSAL2+DREL+     LEEFTIJD+HH2CRED+LTV_curr+rate1Y+rate5Y+CIremFIRP+URB+WELSTAN2+OutNot+mover+SavRate+CRate, data = test, family = "binomial")     summary(log)                           [,1]           [,2]      [1,] -1.029560622391 -1.26664018566      [2,] -0.078225500455 -0.06536644222      [3,]  0.410455341173  0.67036107254      [4,]  0.006961510972 -0.11463794856      [5,]  0.473629162069  0.70074482878      [6,]  0.614550199698  0.83536187570      [7,] -0.612837570442 -0.48086112696      [8,] -0.743739495966 -1.06994471577      [9,]  0.200419240204  0.83957097597     [10,] -0.166568966328 -0.50583277715     [11,]  0.017640270701  0.12678131085     [12,] -0.005947704128 -0.04248886193     [13,] -0.428175932694 -1.69521649738     [14,]  0.049657239050  0.26482261363     [15,]  1.602200661890  2.50479250068     [16,]  0.367771764513  0.96127873663


title: How to write an regression equation with mediation?

text: For my master's thesis I am examining a mediated relationship. I need to write out the equation model, but do not know exactly how to do this. In my case I have four independent variables, one mediated variable and one dependent variable. How do I add the mediated variable to the equation? Until now I got this: Exploitative innovation $= \alpha + \beta_1\ $Inspirational motivation $+\ \beta_2\ $ Individual consideration $ +\ \beta_3\ $ Intellectual stimulation $+\ \beta_4\ $ Idealized influence $+\ \beta_5$ Knowledge exchange and combination capabilities $+\ \epsilon, \epsilon \sim$ iin (0, $\sigma$). $\beta_5$ is my mediated variable. Can I just add this or do I have to multiply something?


title: ROC and multiROC analysis: how to calculate optimal cutpoint?

text: I'm trying to understand how to compute the optimal cut-point for a ROC curve (the value at which the sensitivity and specificity are maximized). I'm using the dataset `aSAH` from the package `pROC`. The `outcome` variable could be explained by two independent variables: `s100b` and `ndka`. Using the syntax of the `Epi` package, I've created two models:               library(pROC)     library(Epi)     ROC(form=outcome~s100b, data=aSAH)     ROC(form=outcome~ndka, data=aSAH)      The output is illustrated in the following two graphs: ![enter image description here](http://i.stack.imgur.com/kb5SH.jpg) ![enter image description here](http://i.stack.imgur.com/RA1gv.jpg) In the first graph (`s100b`), the function says that the optimal cut-point is localized at the value corresponding to `lr.eta=0.304`. In the second graph (`ndka`) the optimal cut-point is localized at the corresponding value to `lr.eta=0.335` (what is the meaning of `lr.eta`). My first question is:   * what is the corresponding `s100b` and `ndka` values for the `lr.eta` values indicated (what is the optimal cut-point in terms of `s100b` and `ndka`)? **SECOND QUESTION:** Now suppose I create a model taking into account both variable:               ROC(form=outcome~ndka+s100b, data=aSAH)      The graph obtained is: ![enter image description here](http://i.stack.imgur.com/LbN2p.jpg) I want to know what are the values of `ndka` AND `s100b` at which sensibility and specificity are maximized by the function. In other terms: what are the values of `ndka` and `s100b` at which we have Se=68.3% and Sp=76.4% (values derived from the graph)? I suppose this second question is related to multiROC analysis, but the documentation of the `Epi` package does not explain how to calculate the optimal cutpoint for _both_ variables used in the model. My question appears very similar to this question from reasearchGate, which says in short: > The determination of cut-off score that represents a better trade-off > between sensitivity and specificity of a measure is straightforward. > However, for multivariate ROC curve analysis, I have noted that most of the > researchers have focused on algorithms to determine the overall accuracy of > a linear combination of several indicators (variables) in terms of AUC. > [...] > > However, these methods do not mention how to decide a combination of cut-off > scores associated with the multiple indicators that gives the best > diagnostic accuracy. A possible solution is that proposed by Shultz in his paper, but from this article I'm not able to understand how to compute optimal cutpoint for a multivariate ROC curve. Maybe the solution from the `Epi` package is not ideal, so any other helpful links will be appreciated.


title: How do I handle measurement error in sparse data?

text: Under-educated non-statistician seeks short-term relationship for very one- sided benefit. System in question involves rocks and physical properties. Modelling bits of the earth typically means few real measurements compared to volume of model. Many estimates are required, and I have NO IDEA how to handle the statement of uncertainties. For illustration let's say the process flow of modelling is 'measure small sample of rocks: get property X (via mean of measurements)', then 'use X in simple model to determine Y'. E.g., Y = mX + b. If I have 20 rocks that I measure, and get a mean for property X, how do I show its uncertainty to start with, and then how do I propagate that through in my calculation of Y (assuming the uncertainty around m is insignificant in comparison)? I've looked at my data graphically, and they seem to have mostly Gaussian shape. E.g. this is a plot of the kernel densities of property X: ![image kernal](http://i.stack.imgur.com/7oM6l.png) Some skew and bumpiness evident, possibly affected by small sample size. OK so there's a fairly large spread of values for some of the rock types, but in essence we deal with large volumes, and have to idealise our models quite a bit, so tails like you see on the yellow curve far left, while real, are not going to be explored in our numerical simulations. I was once told by an older pro that if my model can explain about 85% of observations then I should buy champagne. This: Estimating error from repeated measurements seems to be a similar question, but I don't even really understand the accepted answer. Std dev/sqrt(n) is the 'standard error'?


title: Interpreting AIC and deviance values

text: I have fitted a glm Poisson to my frequency data and obtained the result:      Null deviance: 657.49 on 583 degrees of freedom Residual deviance: 575.00 on 571 degrees of freedom AIC: 1534.4 Is the high AIC value and relatively high residual deviance a serious matter to be taken into consideration? Thanks


title: interpreting scatter plots

text: I have plotted my predictors versus my response variable in order to be able to fit a good Poisson model to my data. However i am getting those scatter plots for nearly all predictors which make me rather worried since including these predictors is making my residual deviance a value of 575 on 571 df(high enough) Thanks ![scatterplot](http://i.stack.imgur.com/42Fbf.jpg)


title: Definition of likelihood ratio test

text: Do I understand this correctly: You fit a model with unknown parameters to a dataset. You choose the parameters so the likelihood of the dataset under the model is maximal. Let this be $L_{max, model1}$. You do the same with a more advanced model. Again you pick the coefficients of the model so that the likelihood of the dataset under the model is maximal. Let this be $L_{max, model2}$ To see which model is the best, you compare the two (maximum) likelihoods. Now what confuses me is the definition of the likelihood ratio test: > "The likelihood ratio is the ratio of the likelihoods of the two models... > ". Shouldn't that be "The likelihood ratio is the ratio of the maximum likelihoods of the two models..."?


title: How to handle missing data in all explanatory variables in linear regression

text: It's multivariate linear regression meaning multiple dependent variables (Y). Data A has both X (explanatory variables) and Y, but data B has only Y. I wonder if there is any way to incorporate data A and B into the regression model. If I remember correctly, Statistical Analysis with Missing Data by Little et al, suggested some kind of iterative approach like EM algorithm where you repeat following until convergence   1. Estimate missing X based on current regression model   2. Compute regression model with full X estimated in #1 But in this case, I have all columns of X matrix is missing. So my questions:   1. Is the approach above still effective? That is, would including data B into the model help in any way?   2. If so, how missing full columns of X can be estimated? Update: The reason I want to add data B is because data B is application- specific data while data A is more general data. I hope to have training data that has more "weights" for application-specific input space.


title: adk.test: sample size

text: I'm using the function `adk.test` (AD test, adk package) to estimate if a data sample follows a distribution family. My question is I have `data`with 100 values, should I simulate the distribution with 100 values or with as many as I want? `rnorm(100)` or `rnorm(1000000)` EDIT: using `fitdist` and `gofstat` functions of `fitdistrplus` package I get the following results: Ex: `Normal distribution`               ad = 1.4491191816923,      adtest = "rejected"      How can I check the validity of the function conclusion? I must compare `ad` with.... Also for the same `data` I get very different conclusions using `adk`:               1st quartile of 100 runs of `adk` with 1000 samples:     p-value = 0.20     min p-value = 0.05     conclusion: the null cannot be rejected


title: Clustering large and sparse datasets

text: I have a few datasets of "interactions" between pairs of elements like so:               element1 element2 1     element2 element3 1     element4 element5 1     ...     element505535 element4 2      where the value in the 3rd column is the "strength" of interaction. Almost all of these strengths are "1." A strength of 1 means that this interaction was observed one time. A strength of 2 means 2x, etc. I have actually gone one step further and normalized all of my datasets by the total # of interactions observed in the dataset so that datasets interaction values can be compared. There are 5-6 million interactions listed in each file and each dataset is obviously under-sampled since there are ~500k elements (making a square matrix of ~250 trillion positions). I would like to cluster these datasets so that I can make statements about which types of elements tend to cluster with which other types elements. Obviously, robusticity of clustering will be a factor—but this is partially ameliorated by the fact that I will make biological replicates of the data. I have tried a few different "naive" clustering approaches just to see what I could do easily with the data. I fully realize that these are problematic ways of clustering, either because they are not robust or because they rely on the data being very undersampled, but here is what I've done:   1. Clustering elements together as long as there is at least one interaction between each element in the cluster and at least one other element in the cluster. When I do this, all the elements end up in a single cluster. This was important to do because it tells me that there are no pairs of elements that are totally isolated from the rest of the group.   2. Finding "superclusters"—that is, clusters where every member of the cluster interacts with every other member of the cluster (e.g. a triangle for a cluster of 3 and a box with an X in the middle for a cluster of 4, etc). This yields almost exclusively clusters with 2 and 3 elements after about 10% of the data has been analyzed (this is still running). I would love to be able to do some sort of hierarchical clustering using my "interaction strength" values as the distance measure between each pair of elements (unobserved interactions have a strength of 0). Does anyone know of a way to do HC on this sort of large, sparse data—or know of a clustering method that might be more appropriate? I've used R up until now.


title: Convergence of a matrix product

text: Let $A=o_{a.s.}(1)$; $A:k\times k$ matrix and $Vu=O_p(1)$; $V:k\times k$; $u: k\times 1$.   Specifically, $Vu$ converges in distribution to $\mathcal N(0,I_k)$. Can we show that $VAu=o_p(1)$ or $\|VAu\|=o_p(1)$? * * * If $V=O(1)$ and $V^{-1}=O(1)$ then we can show that $VAV^{-1}Vu=o_p(1)$. However, $V^{-1}$ can be unbounded... * * * Related question: Can $Vu=O_p(1)$ and $V=O(1)$ imply $u=O_p(1)$?


title: How to calculate time series seasonality index in R?

text: In R, I use the decompose method on my time series object and it gives me seasonal + trend + random component. For seasonal component, it gives me absolute value which is good but I would also like to know the monthly seasonality index as well (like Jan .084, Feb 0.90, Mar 1.12, etc., for example). Is there a quick way to get this seasonality index in R? Also, how is it normally calculated?


title: Metropolis Random Walk Algorithm in Python

text: I am trying to implement random walk metropolis in Python for generate from a pdf $f$, which is defined in the domain: $$(0, \infty) \times [3, 4] \times (0, \infty)$$ I am using Metropolis Random Walk with a trivariated normal. Can someone just take a look and see if my algorithm is correct?               def MRW():     #Covariance Matrix     matrizCov = np.array([[5,0.0,0.0], [0.0,2,0.0], [0.0,0.0,0.6]])          #Array with will contain my chain     global chain          #Inicial state is the maximum likelihood point of my function     cadeia = np.array([[1.25,3.28,3.54]])          y=np.array([[1.0, 1.0, 1.0]])           #Generate 999 other chain states     for t in range(999):         y[0]=np.random.multivariate_normal(mean = chain[t], cov=matrizCov)         #Generates a new candidate for chain state while it is not on domain         if(y[0][0]< 0 or y[0][1] < 3 or y[0][1] > 4 or y[0][2]<=0 or f(y[0])<0 ):             a = 0         else:             a = min(1, f(y[0])/f(chain[t]) )              u = np.random.uniform() #Generates uniform number in (0,1)              #Verify if y[0] will be accepted         if(u <= a):#the case where y is accepted             chain = np.append(chain,y,axis=0)          else:#the case where y is not accepted             chain=np.append(chain,np.array([chain[t]]), axis=0)       I would be glad if someone could just take a look at my code Thanks!


title: Common weak learners for Adaboost

text: I'm looking for a **set of weak classifiers** that work with **Adaboost** to test on **popular datasets**. Most of the examples on the web use some kind of random weak learners which work on their own randomly generated dataset. Could you point me to any usable weak learners?


title: Sample size calculation for logistic regression using simulations

text: I want to find the sample size for logistic regression where I have a covariate with 15 levels and the covariates interacts with time, which means that the effect of the covariates is different for different periods of times. Can anyone help me find the sample size and effect size? I wish I could do this using simulations. Is there a code in SAS or R?


title: Testing overlapping binary and categorical data

text: I need to analyze some data obtained from "in vivo" software testing. Unfortunately I’m a SW engineer, not a statistician. The dataset contains observations for a number of messages sent between 2 nodes (actually, sent and returned). There are two independent possible outcomes: 1. “message was flagged” and 2. “message lost integrity”. (I’m also interested in analyzing “1 OR 2” and “1 AND 2.”) For each message, I have: a) Binary, overlapping data, describing system features activated or not (for instance, the system may have been setup with the “total ordering” feature only, or “total ordering and encryption X” features, or none, etc.). There’s a set of 10 such features, activated or not, for each message. These parameters are completely independent from each other. b) One parameter of categorical data: in which of 5 locations was the message analyzed. c) Finally, I have four parameters for message category (from automated semantics detection), which might be set to “yes”, “no” or “couldn’t tell”. For instance: (urgent, yes) (coordination, couldn’t tell) (critical, yes) (illicit, couldn’t tell). I’m trying to find out if there are combinations of these parameters that make it more likely for the possible outcomes to occur. For instance, activation of features 3 and 4 at location X for messages categorized as urgent and non- coordination predicts that messages loose integrity. I’ve looked at some techniques and already tried chi-squared for independence for each variable. This approach, however, can’t be right: I need combinations of parameters that explain the outcomes. How should I go about getting such thing?


title: Categorical variables clustering with R

text: I wonder whether it is possible to perform within R a clustering of mixed data variables. In other words I have a data set containing both numerical and categorical variables within and I'm finding the best way to cluster them. In SPSS I would use two - step cluster. I wonder whether in R can I find a similar techniques. Thanks in advance. I was told about poLCA package, but I'm not sure...


title: mathematics books for studying Econometrics PhD

text: I started my PhD study in Econometrics last month. Currently I am studying some compulsory course works (Probability theory and Statistical Inference), and also reviewing the literature in my field. I would also like to study some more mathematics, which I believe will make me more prepared for my formal research. Can anyone suggest some math (text)books ? Thanks a lot.


title: Resampling with SAS/IML

text: I am defining the basic matrix:               1         2         3         4     5         6         7         8      As follows:               PROC IML;         RESET NOPRINT;         matrix = {1 2 3 4, 5 6 7 8};     EXIT:      I'm looking to resample each row (for using the bootstrap) with replacement, so we generate matrices like:               1         3         4         2     6         6         5         5      and               1         1         3         2     5         6         6         8      How can this be done?


title: Transformation of random variables

text: What is the probability distribution function of the variable $y$ given by $$y=\frac{x_1}{x_1-x_2},\quad \: x_i\ge 0,$$ given that $x_1$ and $x_2$ are independent and identically distributed and $$x_i = c+z_i,\quad i=1,2 $$ where $c$ is a real nonnegative constant and $z_i\sim\chi^2_\nu$.


title: Test for one category

text: I am doing a research in linguistics and I am not to sure what test to choose in the following problem. There is a group of 26 students. They are asked a question and two possible variants are suggested. 12 of them choose the first answer, and 14 of them choose the second variant. What test should I use to statistically determine whether students favor one answer over another?


title: How to fit a standardized generalized hyperbolic distribution to my data using optim?

text: I want to fit a standardized generalized hyperbolid distribution to my data using the optim command. Regarding my post here, I think I have to try it manually, by using the optim command. My problem is, that the optimization command leads to parameter values, which are not defined, so I get the error message:               Fehler in if (delta <= 0) stop("delta must be greater than zero")       See the following code (my data can be found here):               library(fBasics)     startvalue<-c(1,0,1)     loglikstandghyp <-function(zetapar,rhopar,lambdapar){     if(par>0) return(-sum(log(dsgh(mydata,zeta=zetapar,rho=rhopar,lambda=lambdapar))))     else return(Inf)     }          #shape parameter zeta is positive, skewness parameter rho is in the range (-1, 1).      optim(startvalue, fn=loglikstandghyp, method="BFGS")     param = optim(startvalue,loglikstandghyp, method="BFGS")$par      This leads to the result, that optim sets the parameters to a value, which are not allowed. So my question is, how can I to this correctly? R says: shape parameter zeta is positive, skewness parameter rho is in the range (-1, 1). I cannot use Brent-method to give the range allowed for the parameters, because Brent is only for one-dimensional problems. So the problem is to use the optim command, but I have to give some restrictions, how can I do this? Afterwards, I would use the .paramGH(zeta,rho,lambda) command to change the parametrization into the desired alpha, beta, delta, mu, lambda parametrization. Do you think what I am doing is correct? EDIT: Even if I use other starting values and the L-BFGS-B method like in the following, I still get the error message:               startvalue<-c(20,0,1)     optim(startvalue, fn=loglikstandghyp, method="L-     BFGS-B",lower=c(0.00000001,-1,-Inf),upper=c(Inf,1,Inf))      EDIT: Ok, I now tried the following:               optim(startvalue, fn=loglikstandghyp, method="L-     BFGS-B",lower=c(0.00000001,-0.999999999,-Inf),upper=c(Inf,0.99999999,Inf))      This seems to be better, but my problem now is, that the solution is highly dependent on the startingvalues and I don't know, what starting values are appropriate? E.g. If I set the starting values to `startvalue<-c(100,-0.5,90)` I get:               > optim(startvalue, fn=loglikstandghyp, method="L-BFGS-B",lower=c(0.00000001,-0.999999999,-Inf),upper=c(Inf,0.99999999,Inf))     $par     [1] 100.0189990   0.3125539  90.0142759      And If I set them to `startvalue<-c(50,0.5,1)`:               > optim(startvalue, fn=loglikstandghyp, method="L-BFGS-B",lower=c(0.00000001,-0.999999999,-Inf),upper=c(Inf,0.99999999,Inf))     $par     [1] 0.00000001 0.09105067 8.27463528      Totally different, so what should I do?


title: Locating a set of 5 geometrically constrained points within a cloud of points

text: What is the most efficient way to find a set of 5 rigid, coplanar, non- collinear points within a 3D point cloud of, for example, 100 points? (1) The configuration is specified by 2D coordinates (or in 3D with z=0). The points in the configuration and the ones in the point cloud we are intending to search are specified in metric coordinates. (2) The search has to be repeated for multiple point clouds for the same configuration. (3) The goal is to find the quickest method. (I tried this same question in stackoverflow, mathematics, and theoretical computer science forums, so I thought I'd try here in case this was a problem related to data visualization.)


title: Is ARIMA model appropriate for this dental research?

text: Please, I have a doubt in my study for doctoral thesis… Title: “Development of a thermal cycling protocol for dental materials”. Objective: to create a protocol for thermal aging, from measurements that were found in people mouth. Doubt: What analyses could I use for this variation temperature data? Experimental unit: six volunteers Variable: temperature A student of PhD in dentistry did a work to describe an average pattern of temperature (°C) variation in the oral cavity of six volunteers. She collected oral temperature data (°C) for 3 days, in 6 patients at each 4 of 4 seconds. The aim is to establish an average pattern of temperature variation, according to the time, temperature variation and the number of changes per day. The patients have different routines of each others and per day, so, is there any statistical analyses that I can establish a standard protocol with her data? Basically, she would like to find an answer like as: A person, on average, per day, has three temperature increases of 40°C , which was maintained for 40 seconds, and 3 temperature decreases from 16°C maintained for 40 seconds. Is ARIMA test appropriate for this study? If not, is there another solution it?


title: Test if differences between frequencies is significant

text: I have the following frequency table:               35  0   4   3   7   6   5   4     39  1   9   6   7   7   6   8     36  0   7   10  11  11  10  16     41  0   9   8   8   7   6   7     41  0   8   9   10  9   12  11     55  2   12  9   11  12  11  13     55  1   10  10  11  10  12  11     47  1   14  8   12  15  12  12     45  1   10  11  10  10  9   18     56  0   13  16  12  12  12  11      The Kruskal-Wallis ANOVA test returns:               Source      SS      df     MS     Chi-sq   Prob>Chi-sq       Columns  25306.8    7   3615.26   47.16    5.18783e-008       Error     17083.2   72   237.27     Total     42390     79      According to a multiple comparison of mean ranks:   * Six groups of mean significantly different from group 1 (column 1)   * Six groups of mean significantly different from group 2 (column 2) * * * Now the Kruskal-Wallis and multiple comparison tests make sense, however the Chi Square Test returns a chi square value of 31.377 and a p-value of 0.9997, which leads us to accept the null hypothesis that the frequencies are independent. I understand that an assumption of ANOVA is independence, but... I want to see test if the frequencies are statistically independent, was the Kruskal-Wallis and multiple comparison tests the correct methodology? Note: I am not trying to be subjective, but for a given set of frequencies, how do you test that the differences between groups are significant?


title: How to do rank prediction using R?

text: I'd like to do rank prediction (something very similar with regression) based on historical data, is there any package I can use in R ? Here's my problems: I have a historical data of sports games, and all the rank of each team and some statistic of these teams, I'd like to use these data to predict these teams' next year rank. Thanks


title: Statistical faults of the article relating GDP growth with certain non-sensical variable

text: Colleague draw my attention to the following discussion paper. In my opinion this discussion paper is clearly intended as a critique of the routinely published articles in the field. The discussion paper relates the growth of gross domestic product with the length of certain male body part. The relationship is clearly non-sensical, but the author lays out quite convincing argument backed by data and reasonable (at first glance) statistical analysis. My question is are there any statistical faults in author's arguments? This paper was mentioned in Freakonomics blog.


title: With regard to ARMA time series, what exactly is eacf (extended auto-correlation function)?

text: I'm going through my copy of Analysis of Financial Time Series, 2nd Edition, and I'm at the ARMA portion. One of the techniques for model selection is computation of the extended auto-correlation function, which indicates that, in a table of EACF's with MA coefficients listed across the top and AR coefficients listed down the side, the top-leftmost EACF that is less than the absolute value of 2 times the standard error of the EACF is at the position that indicates a good choice for the orders of the model. Can someone please explain to me exactly what the EACF is?


title: How to compare features and classifiers which achieve perfect accuracy?

text: So I'm looking to compare different combinations of features and classifiers. But I'm getting a lot of combinations that achieve 100% cross validation accuracy. I'm trying to figure out how I would compare the usefulness of each combination. For example I can both train an SVM using Features 1, 10, 15 to get 100% accuracy. But at the same time I can train a logistic regression classifier only using Feature 7 to get 100% accuracy. Also this is a binary classification problem.


title: What are good data visualization techniques to compare distributions?

text: I am writing my PhD thesis and I've realized that I rely excessively in box plots in order to compare distributions. Which other alternatives do you like for achieving this task? I'd also like to ask if you know any other resource as the R gallery in which I can inspire myself with different ideas on data visualization.


title: Whether to use r-square or adjusted r-square with a small sample size that may represent the entire population?

text: I read online that it is only necessary to use adjusted-$R^2$ when you are working with a sample rather than the entire population. The data I'm working with is information on a series of live educational seminars. Each datapoint represents a single seminar that was held in the past, and contains various information on that program's characteristics. In trying to decide whether to use $R^2$ or adjusted-$R^2$, I can see two different sides to the coin.   1. Since my dataset contains every seminar we've held to date, I'm working with the entire population, so I should go with regular old $R^2$.   2. The population of interest is really all _possible_ seminars, including those that haven't happened yet, especially since my goal in this model is to better understand the relationship of factors going forward. Therefore I am looking at a sample, and I should use adjusted-$R^2$. Which logic is correct, and which measure of correlation should I use?


title: Investigation to assess interaction of two predictors using R software

text: In preparation for my end of semester R exam, I've been working on the following CSV file which was modelled on R, which I'm currently having problems figuring out. ![enter image description here](http://i.stack.imgur.com/cdRNC.png) ![enter image description here](http://i.stack.imgur.com/B6D61.png) To investigate this issue I decided to introduce a bilinear variable (SxPR) which is the multiplication of the variable S (power plant’s capacity) and the categorical variable PR (whether LWR existed on the same site previously) and this can be seen in the excel data. My current model which I've been asked to 'add suitable terms' to in order to investigate whether there is good reason to believe that the impact of a power plant's capacity on construction costs might not be the same for plants that did have a LWR on the same site previously and those that did not is: ![enter image description here](http://i.stack.imgur.com/hZn14.png) This gives the equation $C = -9752+ 140.3*D + 4.868*T2+ 0.4180*S- 86.49*PR+ 153.2*NE- 8.38*N$ So basically, I need to confirm that what I've currently done to solve this problem is correct or if not, what else can I do because I've had a lot of problems with this question. ![enter image description here](http://i.stack.imgur.com/DYxo7.png) I concluded that the results show the bilinear term (SxPR) has a p-value>0.05 which suggests there is no interaction between S and PR (since SxPR is not a significant predictor to the model). This causes me to believe the impact of power plant capacity on construction cost is not the same for the plants that did have a LWR existing on the same site previously and those that did not). To support this argument, the Scatter Plot below shows the relationship between the Construction Cost and Power Plant Capacity depending on whether or not a LWR existed on the same site previously. The crosses are the plots of the Power Plant Capacity where a LWR did not exist on the same site previously and the diamonds where a LWR did exist on the same site previously. ![enter image description here](http://i.stack.imgur.com/I338E.png) ![enter image description here](http://i.stack.imgur.com/M2o13.png) Would anyone agree I'm on the right lines. Apologies for the lengthy question. Hope it makes sense!


title: Data Conversion to Standard data format in hierarchical Dirichlet process

text: I'm trying to test the performance of posterior inference on a set of documents with hierarchical Dirichlet process for topic modeling. How can i convert my data (document) to standard data format ? such as this command :               [M] [term_1]:[count] [term_2]:[count] ...  [term_N]:[count]      where               [M] is the number of unique terms in the document     [count] associated with each term is how many times that term appeared      in the document.


title: Is there a way to explain a prediction from a random forest model?

text: Say I've got a predictive classification model based on a random forest (using the randomForest package in R). I'd like to set it up so that end-users can specify an item to generate a prediction for, and it'll output a classification likelihood. So far, no problem. But it would be useful/cool to be able to output something like a variable importance graph, but for the specific item being predicted, not for the training set. Something like: Item X is predicted to be a Dog (73% likely)   Because:   Legs=4   Breath=bad   Fur=short   Food=nasty You get the point. Is there a standard, or at least justifiable, way of extracting this information from a trained random forest? If so, does anyone have code that will do this for the randomForest package?


title: The weight updating in adaboost

text: 1.AdaBoost updates the weight of the sample By the current weak classifier in training each stage. Why doesn't it use the all of the previous weak classifiers to update the weight. **(I had tested it that it converged slowly if I used the previous weak classifiers to update the weight )** 2.It need to normalize the weight to 1 after updating(just need to multiply factor). I think this step can be omitted in implementing. Right?


title: Interpretation of TBATS Components

text: I've searched everywhere for this answer, and come up empty handed. I'm building a forecasts model for interface traffic reported in 5 minute intervals over a 1 year period. To account for multiple seasonality, I came to the conclusion that the TBATS function in R would be most useful (with seasons of 60/5*24=288 and 7 for 7 days/wk). Ploting the TBATS object produced the following output: ![Tbats test decomposition](http://i.stack.imgur.com/eGxNT.png) And while I think I have a good understanding of what the seasonality means, I'm not sure how to interpret the slope and level. I think that the slope is the first derivative of the trend (so since its consistently positive, the trend is increasing at various rates), but I'm completely confused as to what level represents. Is there anything comparable to the remainder graph in STL analysis Also I'm not entire sure what the time values along the bottom represent. Clarification on that would be greatly appreciated. ![enter image description here](http://i.stack.imgur.com/VQM6d.png)


title: Cross Validation for Ridge Regression

text: I'm using ridge regression for calculating optimal weights of a set of scores. These scores are correlated so the usage of ridge regression is used for penalizing large values of weights. So the purpose of ridge regression is to find beta that minimize the following: $$ \sum_i{(y_i - x^T_i\beta_i)^2} + \lambda \sum_j{\beta^2_j} $$ **My question is** : How do I choose an optimal value for lambda, _in the sense of cross validation_? I'm having trouble grasping this conceptually. In classification cross validation is straightforward- split the data into k folds, train on k-1 folds, predict on the last fold and average the prediction error over all folds. How does this work for regression? I can measure the sum of squared distances over each fold, but this is prone to noisy outliers. The reason for using ridge regression instead of standard regression in the first place was _not_ to minimize this. I looked into the following article but I still don't understand the general approach of using cross validation for choosing an optimal ridge regression model.


title: Which is the best accuracy measuring criteria among rmse, mae & mape?

text: I have created training set and test set from my data. Then I performed auto.arima() and ets() in R on the training set to predict one-step ahead forecasts. These were then compared with the test set values to measure error, namely RMSE, MAPE & MAE. This is the output of both ets and auto.arima                RMSE.ets      [1] 3767.561      RMSE.ar      [1] 3776.308      MAE.ets      [1] 2885.112      MAE.ar      [1] 2624.482      MAPE.ets      [1] 0.04232065      MAPE.ar      [1] 0.03857747      Which criteria should be ideally used to select one of the two models (ets or auto.arima) for future predictions. Or is there any other criteria that I am missing out on. Kindly help.


title: Re-parameterization of an asymmetric s-shaped function

text: Is anyone aware of a re-parameterization of any asymmetric s-shaped function (like, but not necessarily the 5 parameter logistic curve), where one of the parameters is the first inflection point of the first derivative (i.e. the maximum of the second derivative). I mean **point 1** in the upper figure: ![enter image description here](http://i.stack.imgur.com/29SIF.jpg) (The picture shows the point mentioned above for the of case of a symmetric logistic function.) So far, I had a look at various references (among others Ratkowsky, D. A. (1990). Handbook of Nonlinear Regression Models. New York, Dekker). However, I have only found parameterizations, where one of the parameters is the inflection point of the function (point 2 in the upper figure). Unfortunately, calculating this point after the estimation is not a possible solution in my case as this parameter should be integrated in another equation that is estimated simultaneously.


title: R regression output - Factors vs numeric variables

text: Let's say I have the following logistic regression models:                df=data.frame(income=c(5,5,3,3,6,5),                       won=c(0,0,1,1,1,0),                       age=c(18,18,23,50,19,39),                       home=c(0,0,1,0,0,1))          > md1 = glm(factor(won) ~ income + age + home,      +           data=df, family=binomial(link="logit"))     > md2 = glm(factor(won) ~ factor(income) + factor(age) + factor(home),      +           data=df, family=binomial(link="logit"))     > summary(md1)          Call:     glm(formula = factor(won) ~ income + age + home, family = binomial(link = "logit"),          data = df)          Deviance Residuals:            1        2        3        4        5        6       -1.0845  -1.0845   0.8017   0.4901   1.7298  -0.8017            Coefficients:                  Estimate Std. Error z value Pr(>|z|)     (Intercept)  4.784832   6.326264   0.756    0.449     income      -1.027049   1.056031  -0.973    0.331     age          0.007102   0.097759   0.073    0.942     home        -0.896802   2.252894  -0.398    0.691          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 8.3178  on 5  degrees of freedom     Residual deviance: 6.8700  on 2  degrees of freedom     AIC: 14.87          Number of Fisher Scoring iterations: 4          > summary(md2)          Call:     glm(formula = factor(won) ~ factor(income) + factor(age) + factor(home),          family = binomial(link = "logit"), data = df)          Deviance Residuals:               1           2           3           4           5           6       -6.547e-06  -6.547e-06   6.547e-06   6.547e-06   6.547e-06  -6.547e-06            Coefficients: (3 not defined because of singularities)                       Estimate Std. Error z value Pr(>|z|)     (Intercept)      2.457e+01  1.310e+05       0        1     factor(income)5 -4.913e+01  1.605e+05       0        1     factor(income)6 -2.573e-30  1.853e+05       0        1     factor(age)19           NA         NA      NA       NA     factor(age)23   -1.383e-30  1.853e+05       0        1     factor(age)39   -3.479e-14  1.605e+05       0        1     factor(age)50           NA         NA      NA       NA     factor(home)1           NA         NA      NA       NA          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 8.3178e+00  on 5  degrees of freedom     Residual deviance: 2.5720e-10  on 1  degrees of freedom     AIC: 10      So depending on the mode of the predictors, R produced different outputs. For factors, R splits out the coefficients into separate categories for the levels, but not for the model with numeric predictors. I'm wondering about a couple things.   1. Is it ever useful to have the response categories expressed as individual rows?   2. To express the general regression equation, how does one go from a model with the categories expressed in an individual equation to an equation with a single B_i. So, for example, if gender has two coefficients, 3.5 for Male and 2.3 for Female, how does one use that in an equation such that (besides converting them into numeric values): Y = B0 + B1 (Gender)


title: Is the selection of ARMA(0,0) wrong?

text: From my model selection based on information criteria such as SBIC returns me a ARMA(0,0) model. Should I choose it or reject it? and Why? Need some guidance on this.


title: Distribution of an estimated parameter

text: I'm studying by myself, and the text says « If $\epsilon_\alpha$ is the upper $\alpha$-quantile of the distribution of $\frac{\hat{\theta}-\theta}{\hat{\sigma}}$, then $P(\hat{\theta}-\epsilon_\beta\hat{\sigma}\leq \theta \leq \hat{\theta}-\epsilon_{1-\alpha}\hat{\sigma} )\geq 1-\beta-\alpha$ .» However, from wiki, I think we have this $P(\hat{\theta}-\epsilon_\beta\hat{\sigma}\leq \theta)=P(\frac{\hat{\theta}-\theta}{\hat{\sigma}}\leq \epsilon_\beta)\geq \beta$, together with $P(\hat{\theta}-\epsilon_{1-\alpha}\hat{\sigma}\geq \theta)=P(\frac{\hat{\theta}-\theta}{\hat{\sigma}}\geq \epsilon_{1-\alpha})=1-P(\frac{\hat{\theta}-\theta}{\hat{\sigma}}< \epsilon_{1-\alpha})\geq 1-(1-\alpha)=\alpha$, imply that $P(\hat{\theta}-\epsilon_\beta\hat{\sigma}\leq \theta \leq \hat{\theta}-\epsilon_{1-\alpha}\hat{\sigma} )\geq \beta-\alpha$. So where did I go wrong? I'm reading Asymptotics Statistics, by Van der Vaart, on page 326


title: How can I display a survival curve with right censoring marks?

text: I have the following survival data and have constructed a survival plot however cannot mark the right censored points (which makes me think my survival graph is also incorrect) - how can I manage this please? Now answered! - ![survival plot](http://i.stack.imgur.com/V9099.png) ![enter image description here](http://i.stack.imgur.com/EETFq.png)


title: Uncertainty modeling for samples

text: I have instances of data that are composed of 6 attributes. Each attribute takes values between 1-10. Any instance has a utility calculated by some utility function: Example: {1,1,2,5,5,2} , utility = 0.89 In some readings for these attribute values, undefined values may appear in the instance, I mark them as question mark '?' Example: {1,1,?,5,?,2} These question marks represents uncertainty in the instance. How can I model this uncertainty in a statistical or mathematical way ? Thanks


title: Draw balls from a set of color balls: the probability of drawing a color seen before?

text: Recently I am studying a probability problem related to the multivariate hypergeometric distribution. The problem is stated as: > Given well-mixed $n$ balls of $m$ colors, and assume that $n_i$ is the > number of balls with color $i$, where $i \in {1, ..., m}$. So we have > $\sum_{i=1}^{m}n_i = n$ > > Do $k$ random draws without replacement. What's the probability that the > $k$-th draw is a color which has been drawn from the previous $(k-1)$ draws? My naive solution is like this: the total possible ball sequences for $i$ draws would be: $P(n, k)$ where $P(n, k)$ is the permutation of $k$ items from $n$ items. And among these sequences, only the ones with at least one occurrence of the last ball's color satisfy the requirement. So the possible sequences would be the all sequences except for the cases where none of the $(k-1)$ draws contain the last ball's color: $P(n, k) - \sum_{i=1}^{m}n_i*P(n-n_i, k-1)$ So the probability can be calculated as: $1 - {{\sum_{i=1}^{m}n_i*P(n-n_i, k-1)}\over{P(n, K)}}$ Now my questions are:   1. Is my naive solution correct, or there are cases missing or overlapping?   2. I tried to write a computer program to calculate the probability given $n, m$ and $k$, however it turns out to be very slow for large inputs, thanks to the permutations. Is there any simpler (so fast on computation) formula for this probability, like a closed form formula? Any hints would be appreciated, and let me know if the question is too hard to understand~


title: Problem with classifying new observations with discriminant analysis

text: I have a data set of 40,000 individuals which I clustered using k-means. I used 30 variables, each ordinal from `1=minimum` to `5=maximum`. I reduced these 30 variables to 10 factors and ran K-means on these new variables. I kept a clustering of 12 clusters. I now have 7000 new observations, I want to classify these using a discriminant function. I built one in SPSS using the 30 previous variables. I get a good clasification on the original dataset (40,000), but when I run it on the new observations, one segment doesn't show at all. 0 cases are of this type. This makes no sense to me. The original clusters seem healthy and apart, I did them in R with multiple different initial center methods and ran them many many times, and stuck with the best one. The discriminant function has no problem finding this "disappearing" cluster in the original dataset. What could be causing this?


title: randomForest and variable importance bug?

text: I dont get the difference between the `rfobject$importance` and `importance(rfobject)` in the MeanDecreaseAccuracy column. Example:               > data("iris")     > fit <- randomForest(Species~., data=iris, importance=TRUE)     > fit$importance                       setosa  versicolor   virginica MeanDecreaseAccuracy MeanDecreaseGini     Sepal.Length 0.027078501 0.019418330 0.040497602           0.02898837         9.173648     Sepal.Width  0.008553449 0.001962036 0.006951771           0.00575489         2.472105     Petal.Length 0.313303381 0.291818815 0.280981959           0.29216790        41.284869     Petal.Width  0.349686983 0.318527008 0.270975757           0.31054451        46.323415     > importance(fit)                    setosa versicolor virginica MeanDecreaseAccuracy MeanDecreaseGini     Sepal.Length 1.277324   1.632586  1.758101            1.2233029         9.173648     Sepal.Width  1.007943   0.252736  1.014141            0.6293145         2.472105     Petal.Length 3.685513   4.434083  4.133621            2.5139980        41.284869     Petal.Width  3.896375   4.421567  4.385642            2.5371353        46.323415     >       I get different MeanDecreaseAccuracy values but have the same order for the importance Variables (for `fit$importance` as well as for `importance(fit)`):   1. Petal.Width   2. Petal.Length   3. Sepal.Length   4. Sepal.Width But in other datasets I sometimes get different orders. Can someone explain what is happening here? Is this possibly a bug? Thank you in Advance. * * * **Edit** (in response to Martin O'Leary) Okay Thanks! I noticed something else. Taking a Look at the `rfcv()` function i noticed the line:               impvar <- (1:p)[order(all.rf$importance[, 1], decreasing = TRUE)]      with this line we choose the first column of `all.rf$importance` which gives us the order of the class-specific (for the **first factor** ) measures computed as mean descrease in accuracy only. This has not always the same order as the mean descrease in accuracy over all classes (`MeanDecreaseAccuracy`). Wouldn't it be better choosing either the `MeanDecreaseAccuracy` or `MeanDecreaseGini` column, or better using the `importance()`-function for the scaled values? So we would have a sequentially reduced number of predictors ranked by variable importance (over all classes) and not only ranked by variable importance for the first class.


title: Which distribution is likely to be the source of the largest value?

text: Suppose I sample many times independently from each of two normal distributions. Given the means, standard deviations, and number of samples from each distribution, how can I calculate the probability that the highest sampled value came from the first distribution? I can solve this question by simulating it numerically, but haven't figured out an analytic expression.


title: What is the distribution of an average of Poisson random variables?

text: If I have random variables $X_1,X_2,\ldots,X_n$ that are Poisson distributed with parameters $\lambda_1, \lambda_2,\ldots, \lambda_n$, what is the distribution of $Y=\left\lfloor\frac{\sum_{i=1}^n X_i}{n}\right\rfloor$ (i.e. the integer floor of the average)? A sum of Poissons is also Poisson, but I am not confident enough in statistics to determine if it is the same for the case above.


title: How to reconstruct database from cross-tabulated data?

text: I am given a set of almost 100 crosstabulations (most 2D, but some three dimensional) based on a single, yet unavailable dataset. My task is to perform basic statistical tests ($\chi^2$, Mann-Whitney, Kruskal-Wallis, etc.). I could calculate them by hand, by it would be over one order of magnitude faster to use my "streamlined" procedures which automatically make beautiful Excel figures, tables and stuff. Unfortunately my tools use database as input and I cannot feed them with summaries. The crosstabs do overlap, i.e. if one crosstab is between var1 and var2, the could well be between $var1 \times var3$. Is there any existing software (e.g. in R) that I can use to create any database (dataset) which agrees with all given crosstabulations? (I know, that reconstructing database from cross tabulations is not a unique problem. I know that everyone will be better off having a proper database to start with.)                                                        O  O  O      I suspect, that there is no ready solution, so I set writing one: Every crosstab can be written as a set of linear equations (one equation for every known frequency, coefficients equal to 1 or 0), so from mathematical point of view I see the problem equivalent to finding integer solution to set of linear equations. It is a problem from integer programming domain, just like mr. whuber suggested. Unfortunately there is a problem in prohibitive size of the equations; the only way I can think up is by partition the dataset into $\text{no}(var1) \times \text{no}(var2) \times \text{no}(var3)$ groups, where $\text{no}(var)$ denotes number of distinct levels (groups) of variable var. It might be fine for this simple example, but in my case, when I multiply all the numbers of levels of each variable I come up with astronomical number in order of $10^{25}$. And I suspect that sparse matrices optimization wouldn't help, so all I wrote might just be the dead end. Does anyone have any idea about how to solve this problem?


title: PCA and Linear Regression

text: I want to run PCA on a set of variables and then regress the scores on my dependent variable. I have the following questions:   1. Should I scale and center my variables?   2. If yes, should I also standardize my dependent variable before running linear regression?   3. What if I don't standardize my dependent variable?


title: Directional Derivative of a function containing an Indicator function - Optimiality condition for quantile regressors

text: I'm trying to understand a passage in Koenker's Quantile regression book (p.33). It says: (note that y,x, are vectors and w is the direction vector) ![enter image description here](http://i.stack.imgur.com/KZVlI.png) With the first part of the outcome no problem: I apply the product rule for derivatives and the derivative of the indicator function is 0 since the discontinuity of the indicator function doesnt happen in y-x'b different from zero. But what happens in y-x'b=0? My logical mind tells me that there is no derivative for the indicator function if not some complex approximation, so how came Koenker up with that result? Thanks in advance!


title: Uncertainties on parameters from maximum likelihood fitting

text: I am trying to fit a model to binomially distributed data, which I have done via the maximum likelihood method. Normally I'm working with normally distributed data (or data which I can convince myself and occasionally others is normally distributed) and I'm a little confused about how to find asymmetric errors on my parameters. Usually I find uncertainties by scanning the space of the fitted parameters and finding the surface where the X^2 is one greater than the best fit value. I have found some unproven statements suggesting that this method is applicable regardless of whether one is using Gaussian maximum likelihood or not. Is this the case?


title: Simulation Analysis of a Cox Survival Model with Change Point.

text: I wish to simulate survival data for the following probabilistic model which will be analyzed using a Cox model: An exposure $X$ is modeled as binary having $X \sim_{iid} \mbox{Bernoulli}(p)$ where $p = \frac{1}{3}$. The hazard is a function of $X$ with: \begin{eqnarray} \log \left( \lambda(t|X=0) \right) &=& \beta_0\\\ \log \left( \lambda(t|X=1) \right) &=& \beta_0 + \Delta\cdot\mathcal{I}(t \geq \tau) \end{eqnarray} I also wish to generate a censoring variable having an independent exponential distribution.