title: Where to find graphical explanations for various classification models?

text: There are several good books about machine learning and there are many mathematical explanations for several classification models (in most books personally read commonly explain Naive Bayes, Logistic Regression and SVM). Besides the Naive Bayes, Logistic Regression, SVM and several decision trees for which I understand both mathematically and graphically, by observing data instances as scatterplot, there are other existing models for which (so far) haven't found graphical explanation. When dealing with classification a good starting point is to draw scatterplot for data instances and then decide which algorithm to use by observing how the points are distributed in space and which fuction would split observing data into classes the best, therefore I would be glad to find any good resource with the collection of gpraphically explained classification algorithms, perhaps on some real examples. In case anyone of you knows for some useful resource of this kind I would be very thankful for sharing the address/book or article title.


title: Testing median, and distribution of Binomial

text: Given the median is $m>0$ We have $n$ random variables: $Y_1,Y_2,..,Y_n$ If we observe a value above $m$ then $Y_i=1$ if we observe a value below $m$ then $Y_i=0$ The probability of getting an observation above $m$ is "$p$" * * * $(i)$ Then what distribution is $Y_i$ describing? I think $Y_i$ is describing Bernoulli trials. $(ii)$Does that mean that if we want a model of 'number of observations above median' this would be $\sum^{n}_{i=1}Y_i$ which is $Binomial(n,p)$? $(iii)$Is the assumption p=1/2, valid? Since half the observations are above $m$ and half below. * * * $(iv)$I now wish to test if $m$ is the true median, and require a 95% CI. If I know '$n$' then I should calculate exact binomial probabilities until I have cumulative probability in the lower tail below 2.5% and cumulative probability in the uppertail of below 2.5%. Alternatively I may use the normal approximation. * * * Essentially, I am looking to understand steps from $(i)$ to $(iv)$: 1) the validity/ accuracy of my method. 2) assumptions I may have made. 3) any other more accurate/ efficient tests out there? * * * Any input: thoughts/ comments/ hints would be greatly appreciated. Thank you


title: How to deal with death in disease-free survival analysis?

text: If I have disease free survival data (defined as whether or not a particular disease has been diagnosed or not along with the time to that event or loss to follow up) and also overall survival data, how do I deal with deaths that occur without the disease event? Are these censored or should I exclude such patients from the disease-free survival (dfs) analysis? I plan to run dfs analyses for several particular types of disease separately.


title: What is the number of subjects needed to be representative for a small, finite population?

text: I have been studying an organisation with 58 people with a number of Likert style questionnaires. The data are not random samples but rather incomplete census's. For the 6 questionnaires I managed to get between 49 & 57 respondents. However when I come to combine questionnaire data, the number can fall to say 43 subjects. The responses to the questionnaires have been uniformly non-normal and generally leptokurtic. NP bootstrap regression on this number of subjects gives significant results, however I am unsure whether this enables me to say anything about the organisation more generally. If you have any guidance on the limits of what I can say with this sort of data it would be helpful. Many thanks for any advice.


title: Self-study question on Margin of Error

text: I was attempting a self-study question on Margin of Error and managed to get the answer right. The answer key states 2.77, but I managed to get at 2.76. Close enough I guess. But I thought that it seems pretty weird that I did not have to use many of the provided variables such as mean or even the confidence interval. Although I managed to get the right answer, I am not sure if I did really answer the question correctly. Appreciate some guidance please. **Question** > A portable DVD manufacturer based in Shanghai provides battery perfor- mance > information in order to inform customers when making a purchasing decision. > One of the important measures of battery performance is the length time > average batteries last compared with a standard battery with a base of 100. > This mean that a set of batteries with a grade of 200 should last consumers > twice as long, on average, as a set of batteries graded with a base of 100. > A Chinese manufacturer’s organization wants to estimate the actual battery > life of brand-name batteries that claim to be graded 200. A random sample of > n = 400 indicates a sample mean battery life of 195.3 and a standard > deviation of 21.4. A 99% confidence interval for the population mean battery > life for batteries produced by this manufacturer under this brand name is > constructed. Find the margin of error in this estimation. **My Attempt** > ![enter image description here](http://i.stack.imgur.com/Ajv7a.jpg)


title: Bootstrapped regression with total data or bootstrap with matched data?

text: I'm investigating the effect of a continuous variable A on a measurement variable M stratified by another factor variable C in an observational dataset. Due to heteroscedasticity I decided to use a bootstrapped regression analysis. However looking at the data, the background set of variables are not evenly distributed if I dichotomise A (present or not). I've just finished running another analysis where I do the same analysis after having matched the dataset for confounders (using CEM in R). Now the problem is which analysis to trust: the bootstrapped regression approach on the entire dataset or the bootstrapped version of the matched data? Under one of the factors in C the results diverge. Any ideas how this can be analyzed?


title: Reporting results of linear mixed-effects model

text: Linear mixed-effects models aren't commonly used in my corner of biology, and I need to report the statistical test I used in a paper I'm trying to write. I know that awareness of multilevel modeling is starting to appear in some areas of the biosciences (A solution to dependency: using multilevel analysis to accommodate nested data), but I'm still trying to learn how to report my results! My experimental design, in brief:   *Subjects were assigned to one of four treatment groups   *Measurements of the dependent variable were taken on various days after the start of treatment   *The design is unbalanced (unequal numbers of subjects in the treatment groups, and missing measurements for some subjects on some days)   *Treatment A is the reference category   *I centered the data at the final day of treatment  I want to know whether Treatment A (the reference category) yields significantly better results than the other treatments (at the end of treatment). I did my analysis in R, using nlme:               mymodel <- lme(dv ~ Treatment*Day, random = ~1|Subject, data = mydf, na.action = na.omit,      + correlation = corAR1(form = ~1 |Subject), method = "REML")      And the output (in part; truncated for brevity) is:               >anova(mymodel)                   numDF denDF  F-value p-value     (Intercept)      1   222 36173.09  <.0001     Treat            3    35    16.61  <.0001     Day              7   222     3.43  0.0016     Treat:Day       21   222     3.62  <.0001          >summary(mymodel)     Fixed effects: dv ~ Treatment * Day                             Value Std.Error  DF  t-value p-value     (Intercept)         7.038028 0.1245901 222 56.48945  0.0000     TreatmentB          0.440560 0.1608452  35  2.73903  0.0096     TreatmentC          0.510214 0.1761970  35  2.89570  0.0065     TreatmentD          0.106202 0.1637436  35  0.64859  0.5208      So, I know that the effect of Day differs by Treatment, and that, on the final day of treatment (where the data are centered), dv is significantly different in Treatment A than Treatments B or C. What I want to say is: "As predicted, we found that Dependent Variable was significantly lower in subjects receiving Treatment A (mean +/- SE) than in subjects receiving Treatment B (mean +/- SE, p=0.0096) or Treatment C (mean +/- SE, p=0.0065), as measured on the final day of treatment." But, I have to indicate what statistical test was done. Would this be an acceptable way to describe the analysis? "[Measurement Method] was performed on days indicated and Dependent Variable (units) was determined; we analyzed the log-transformed data using a linear mixed-effects model centered at [final day of treatment]. Symbols represent mean dv; error bars are standard error. On the final day of treatment, dv was significantly lower in Treatment A (mean +/- SE) than Treatment B (mean +/- SE, p=0.0096)..." Specifically,   *Does that say enough about the statistical test used? (Readers are used to seeing something more like "mean +/- SE, p=0.0096, Student's t-test," but it seems weird to write "p=0.0096, coefficient for Treatment B vs. Treatment A from linear mixed-effects model at [final day of treatment].")   *Is there a better way to put this? (The methods section will include more information about the stats: "[Measurement Method] data were analyzed using R and the R packages... We analyzed the log-transformed Dependent Variable data by using linear mixed effects models using Subjects as random effects and an autocorrelation structure of order 1 (AR1). As fixed effects, we included Treatment and Day, and the interaction of Treatment and Day. We checked for normality and homogeneity by visual inspections of plots of residuals against fitted values. To assess the validity of the mixed effects analyses, we performed likelihood ratio tests comparing the models with fixed effects to the null models with only the random effects.") Any advice on how to report results of a linear mixed effects model for an often stats-averse audience (and written by a relative stats novice) would be greatly appreciated!


title: Non normal distribution for 4-way mixed ANOVA

text: I am analysing a study using a 4-way mixed ANOVA. It has one between-subjects factor and 3 within-subjects factors (with 2, 2, and 3 levels). My sample size is 63 - 29 in one group and 34 in the other. I have 9 dependent variables, all which asked participants to choose their answer on a scale of 0 to 10. I ran normality tests on my data, and since it is a 2 x 2 x 2 x 3 ANOVA, I had to run 24 tests per dependent variable. The histograms from most of the data don't look very normal, but it is a particular problem from 2 groups of data, since participants were more likely to enter '0'. This means that transforming the data has minimal effect, since the value remains the same at '0'. Since there are no non-parametric tests for mixed ANOVAs, what can I do? Thanks for any help you can give me, I really appreciate it. I am really stuck and have spent so much time looking at different questions but have not been able to find my answer yet.


title: How to deal with a non-gaussian distribution & heteroscedasticity

text: I am working on my thesis project and have come across a problem with the statistics which I am looking for a bit of guidance. I am running ANOVA tests to determine significance between groups but i have come across a problem with a couple of the results which do not have equal variance (p<0.05 & p<0.001). This alone wasn't particularly a problem (Welch correction and unequal variance post-tests) But one of my assays shows non-normal distribution & unequal variance:   * I ran D’Agostino-Pearson omnibus test with P<0.05 in one of the groups (however my sample sizes are small 3 datasets with n=7,15,&16).    * In the literature it appears there is some skew in this distribution from the 'healthy' adult population so it is reasonable to assume it is not normal in my groups.   * Bartletts test for homogeneity of variance has p>0.01. It was recommended I run Kruskal-Wallis ANOVA, but: \- This test assumes equal variance, which is not the case. \- Also my sample size is a bit small in the first patient group (n=7), which to my understanding is that Kruskal-Wallis has little power in this case as well? I essentially do not know what to do with this data set! Any help would me most appreciated!


title: Training Restricted Boltzmann Machines with Bias Units

text: I have written some code that basically does the following:- 1: creates a matrix of binary features and a matrix of scaled features ( in the range 0 to 1 ) 2: RBM trains separately on each of the above features matrices to get weight matrices 3: horizontally stacks a random row vector for bias unit weights and the two weight matrices from step 2 to create a single, initial input layer to hidden layer weight matrix 4: matrix multiplies the input features by the matrix from step 3 and feeds forward through the logistic function hidden layer to create a hidden layer output 5: RBM trains on the hidden layer output from step 4 to get a weight matrix 6: horizontally stacks a random row vector for hidden layer bias unit weights and the weight matrix from step 5 to create a single, initial hidden layer to Softmax weight matrix 7: uses the initial weight matrices from steps 3 and 6 instead of random initialisation for backpropagation training of a Feedforward neural network However, I am now wondering whether it's wise/correct to exclude the weights of the bias units from the RBM training?


title: How to model a conditional probability without estimating joint PDF?

text: I've around 3e3 two-dimensional data points,                  x,   d     1  1, 0.1     2  3, 0.1     3  2, 0.2     4  1, 0.5          range(x) = [-600, 600], range(d) = [0, 1]      We are trying to model `d` with a probability model. I believe the distribution of `d` should be dependent on `x`. So I'm up to model $Pr(d|x)$ rather than $Pr(d)$. To do this, I am now using kernel density distribution to estimate the joint distribution of $Pr(d, x)$ and divide it by $Pr(x)$. However, the number of data is enough to model $Pr(d)$ with kernel density estimation while insufficient to model $Pr(d,x)$. So I'm wondering about this: **Is there any model you suggest that can model`d` in a conditional way without modeling the joint distrubution?** I have the following thoughts: I can use SVM to deal with the short of data. However is there a way to train a "conditional SVM"?


title: Is there a way to measure the viability of an experiment?

text: Is there a way to calculate a number that evaluates the viability of an experiment? By that I mean the probability that the observed probability is just chance. I read about the $p$ value, but it doesn't seem to be that, as it doesn't take at all into account the number of experiments. Let's say I have an experience that I know should produce precisely a certain result, and that my experience proves it wrong. How many "tries" should I have done to consider my experience valid, and is there a measure for that? Edit for a bit of clarification: My experiment consists of taking a (computer generated) random integer number between 1 and 4 - assuming it distributing the numbers evenly. Then, people are placed in front of four cards, representing one of these numbers, and have to choose one of them. Once chosed, the computer generates a number and if it as the same it is a success, otherwise a fail. So the expected result is a 25% success rate. What I am wondering is if there exist a standard measurement for the fiability of an experiment, depending on the number of tests executed, so that if for example I see a 30% success rate I can assume that people are actually right 5% more than pure random luck


title: "Zero-inflated" predictors in regression?

text: I know that zero-inflated models (e.g. zero-inflated Poisson or negative binomial models) can be used for dependent variables. I also know that in general there are no assumptions for the independent variables (i.e. predictors) in regression analyses. However, I have a quantitative (continuous or count) predictor which has many (say, 40-60%) zeros. When I used it as a quantitative predictor in regression (linear or logistic) I got a small P value (i.e. P<0.01), but when I used it as a binary predictor (zero or not) I got a P value>0.05. Why did it happen? How do I interpret this result?


title: Relationship between Bernoulli and Normal CDF

text: Is there any relationship between the draw from a Bernoulli with parameter $p$ and the Normal CDF. Specifically is the condition $p>\Phi(x)$, where $x$ is drawn from the standard normal distribution, equivalent to a draw from Bernoulli$(p)$ ?


title: Understanding SVM regression: objective function and "flatness"

text: SVMs for classification make intuitive sense to me: I understand how minimizing $||\theta||^2$ yields the maximum margin. However, I don't understand that objective in the context of regression. Various texts (here and here) describe this as maximizing "flatness." Why would we want to do that? What in regression is the equivalent to the concept of "margin"? Here are a few attempted answers, but none that really helped my understanding.


title: Categorize continuous data effectively (taking into account a response variable)

text: I wonder what are the better approaches to **categorize continuous data** (e.g. _age_ ) than dividing them with the use of quantiles and `cut` function (in `R`). I have heard about using trees to divide data in the way which takes into consideration how a division would differentiate a response variable, but I cannot find any quick reasonable explanation for that. I want to categorize my data with the aim of using them in **multinomial logit model**. Is there any other approach to do it? (Little off-topic: I use `R` so I would be grateful for some package references or something like this.)


title: What to do with p-values when standard errors are obviously biased

text: What am I supposed to do when I want to interpret significances, although I know that standard errors are biased because of wrong error term assumptions? I know that there is the possibility to use White estimators, weighted OLS. But my prof told me not to do so. Maybe some extra information: 1.) I am analyzing an OLS with a whole bunch of dummy variables. 3.) the assumption of normal distributed error terms is wrong (they are t-distributed) and heteroskedasticity occurs. 4.) I am doing cross sectional analysis 5.) I have a lot of Interaction in the model. And most of the variables are non significant (p-values around 0.8, so that the coefficients are close to zero). My prof doesn't want me to get rid of these variables, although they are not significant (He said that this is not a good way, because stepwise elimination is trouble because of deleting wrong variables and choosing the right criteria). On the one hand I understand why there is no way to interpret the significances. But on the other hand it makes interpreting not easier. Sure I can change the model, but I have to do OLS first, before I am allowed to switch!


title: Analysing grade distributions: difference in proportions across years

text: I have data on grades obtained in a course over three consecutive years. The grade distribution includes the number of students who received a fail, pass, credit, distinction and high distinction. I need to work out if the proportion of students who obtain a pass (or fail, credit, etc.) is the same, or different across the three cohorts. How can I do this? What type of test should I use?


title: A model that can fit binary data in R

text: I have data that indicates if the computers in my network are up or down. Data is collected every minute and has a seasonality of one week, looks like this for 3 hosts. Servers are supposed to go down less often than workstations.               ID TIME STATE TYPE       1  60   0     workstation       2  60   1     server       3  60   0     server       1  120  1     workstation       2  120  1     server       3  120  0     server        I am using R and I'm looking for a model that can give me the probability of a host being down or up in a future time interval. Which model do you think would fit the best these requirements? I have thought of using duration models with the survival package, considering the host going up and the host going down as events. What do you think of this ? Are there better candidates?


title: Performing PCA for normal score transformed data

text: In an attempt to improve the results of Bayesian NNT, I transformed the 7 variables that I have into normal scores (subtracted the mean and divided by SD).   Then I used a PCA on the transformed variables to generate new 7 PCs. I used these PCs to run the Bayesian NNT and the classification results improved a little bit. My question is is this statistically valid and applicable or does it violate some fundamental statistical rules?


title: How do I run a two-part model of health care expenditures in Stata?

text: Does anyone have Stata code they could share with me so that I can run a two- part model to look at health care expenditures? The first part of the model will determine whether or not an individual had any visit and the second part will determine how much individuals spent (estimated on the subset of those who had a visit). I know that I need to take the product of these to get the expected cost of health care for individual i and also need to do a transformation for the first part but I am not sure how to do this in Stata. Thank you!


title: Creating an index based on a set of measurements without a target for purpose of rank ordering

text: The problem I'm trying to solve here is very simple but the available data is very limited. That makes it a hard problem to solve. The available data are as follows:   1. I have 100 patients and I need to rank order them in terms of how healthy they are.   2. I only have 5 measurements for each patient. Each of the five readings is coded as a numeric value, and the rule is that the bigger the reading the healthier is the patient. Should I have some sort of doctor's "expert judgement based ranking" I could use that as the target variable and fit some sort of an ordinal logistic regression model trying to predict doctor's assessment. However, I don't have that. The only thing I have is (1) and (2). How would you come up with a simple "scoring" algorithm which would combine those five measurements into a single score which would be good enough (not perfect) in rank ordering patients?


title: Calculating standard deviation from log-normal distribution confidence intervals

text: I have the results of a meta-analysis of 10 studies that reports a combined random effects odds ratio (computed using Woolf's method) and 95% confidence interval of an event happening in one group relative to another: $OR = 7.1\ (95\%\ CI\ 4.4-11.7)$ I'm now building a model that needs to sample around this odds ratio (for the purposes of a probabilistic sensitivity analysis). Given that it's an odds ratio, I'm assuming that it's log-normally distributed and that 7.1 is the mean, but what's the best way to convert the confidence interval to a standard deviation so I can sample the distribution using Excel's `LOGNORMDIST` function? (I've found similar questions for the normal and gamma distributions (From confidence interval to standard deviation - what am I missing? and How to calculate mean and standard deviation in R given confidence interval and a normal or gamma distribution?) and also questions calculating the confidence interval for a log-normal distribution (How do I calculate a confidence interval for the mean of a log-normal data set?), but I can't seem to find how to go the other way round.)


title: Multiplying multivariate Gaussian with a univariate Gaussian

text: Not sure if this question makes sense but I have been spending a lot of time looking at information about expectation propagation and one of the key operations in it as explained by Tom Minka is the projection of product of distributions to the best Gaussian approximation. So, I have a setup where my joint distribution takes the following form: $$ P(w)\times \prod_{i}P(y_i|w) $$ Now the likelihood can be written as: $$ P(y_i|w) = \frac{\phi}{2\pi} \exp^{0.5e_i\phi e_i} $$ where let us assume for the sake of simplicity that $\phi$ is given and $e$ can be calculated from the observation $y_i$. $P(w)$, however, is a multivariate Gaussian with a covariance structure. So, now when I approximate each of these factors in turn, I will have a product over such univariate distribution and then a product over my current approximation to $p(w)$, which is a multivariate Gaussian. I am utterly confused as to how this can be done. I know simple expressions for mean and variance exist when one is multiplying univariate Gaussians but not even sure what to do in this case. Do I need to look at the ith component of this multivariate Gaussian as well? I would appreciate any help/suggestion with this.


title: Different powers of coefficient - solvable within GLM?

text: I have model where coefficient is of different powers: $$\mbox{log} ( \mu_{i} ) = \alpha + \beta x_1 + \beta^2x_2 + \beta^3x_3 + ... + \beta^nx_n \\\ \\\ N_{i} \sim \mbox{Poiss} ( \mu_{i} ) $$ $N_i$ and $x_i$ is the data, $\alpha$ and $\beta$ are model coefficients.   1. **Can this problem be somehow transformed so that it can be solved using GLM model?** I.e. so that I can use the `glm()` function in R.    2. **If not, how can I solve this problem using some frequentist package in R?** I can write model in WinBUGS but prefer to have "frequentist" solution in R, because it is much faster and inference is easier (one has p-values, t-tests, F-tests...). But which function or package use to fit it? * * * If you are interested how this problem was created: actually taking model from this question and adding coefficient $\beta$ for $\mbox{log} ( \mu_{i,j} )$ on the predictor side of the equation. While the transformation of the original model for use in GLM was elegant and without problem, adding $\beta$ coefficient complicated things and leads to the model presented here.


title: Predicting binary dependent gives non-binary predictions

text: I am trying to predict the result of an experiment (binary dependent variable) based on a number of continuous independent variables. When I do this using a largish model (9 main effects + 2 factor interaction) it seems to "work" meaning I get a plot where the predicted values are 0/1. However, if I reduce the model to say just two main effects + interaction, I get predicted values ranging from 0 - 1 (and several points in between). I have tried using probit and logit with `glm()` in R. With the smaller model I get no warnings / errors. With the larger one I get:               Warning messages:     1: glm.fit: algorithm did not converge      2: glm.fit: fitted probabilities numerically 0 or 1 occurred       What causes this behaviour? How can I determine which of the main effects and/or interaction terms are the most important (they all have `***` in the summary)?


title: How to find optimal neural network structure and parameters values?

text: I tried to create a neural network for predicting temporal data. Having selected a different nuber of inputs (number of previous known values), number of hidden units and regularization parameter value, I calculated and visualized the training and test values. As I haven't managed to tune network structure and regularization parameters such that the predicted values would be close enough to actual values, I would be very thankful if anyone came across similar problem and found a way to solve it. Below are visualized values ('x' -predicted, 'o' - actual) for different number of units in hidden layer and regularization parameter lambda: ![enter image description here](http://i.stack.imgur.com/gugW7.png) lambda = 0.02, number of hidden units = 30 (if it is not visible, the predicted are below actual) ![enter image description here](http://i.stack.imgur.com/uXVWO.png) lambda = 0.02, number of hidden units = 35 (predicted are below actual) ![enter image description here](http://i.stack.imgur.com/e6pN0.png) lambda = 0.08, number of hidden units = 35 I wrote the cost function in Matlab The error cost (J) is defined with the following logic:               cost = (yActual - yPredicted).^2;     J = (2*m)^(-1)*sum(sum(cost));      where               yPredicted = [a2]*[Theta2];      [a2] are the output values of hidden layer units and [Theta2] is the weight matrix The output values of hidden layer are calculated with sigmoid function:               a2 = sigmoid([input values]*[Theta1]);      where sigmoid is thefined as follows:               function g = sigmoid(z)     g = 1.0 ./ (1.0 + exp(-z));      I would be very thankful if anyone could suggest an approach or any modifications in order to improve the prediction results. Thank you!


title: with SVM, are incorrectly classified training points support vectors?

text: I always thought $\alpha_i \neq 0$ implied that $i$ was a support vector. However, now that I think about it, I am confused. If a training point is incorrectly classified, it will have a non-null $\alpha$, does it make sense to use it for classifying new points? Do we use it to compute the optimal separating hyperplane? Is the classifier different from one trained without that point?


title: optimization with two factors

text: I am trying to figure out how can I create an optimization function that will select from a data.frame or matrix, the row with the most drawn value (the numbers are not integer, so I need to maybe look for the number with the highest probability of occurrence) and at the same time minimize the second column (that represents the error columns). This is how my data looks like:               dataset=data.frame(data=abs(rnorm(200,0,2)),error=runif(400,.3,.5))      How cand I do this?


title: Logistic regression using 2 time points

text: I am performing logistic regression and would like some help with interpretation. My data was collected at 2 time points; baseline and follow- up. My outcome variable is a psychosocial functioning variable (as an example working or not working) at follow-up. In step 1 I entered the functioning variable at baseline. In step 2 and 3 I entered a number of covariates e.g sex, age, SES, psychopathology variables (dichotomous). In step 4 I entered a continuous variable which is the BPD difference score and it is dimensional. So my problem is with interpretation. If I find that step 4 was significant and BPD difference score significantly predicted "working at time 2". If this was a linear regression I could say that there was an improvement in working as working at baseline was controlled for. But I am unsure of how to interpret it when I control for the dichotomous variable in step 1 (working at baseline) in a logistic regression. Does that make sense? I have done the analysis in SPSS. Thank you for any responses and I apologise for my inexperience and possible simple question, I am trying my hardest to understand it but I cannot manage to figure it out without help.


title: Can I use Empirical Cumulative Distribution Function to derive standardized values of various data sets?

text: I have an idea to use ECDF conversion of data to their uniform distribution of equal sample sizes. So, for say 1000 pieces of data, each value should more or less correlate with a .1% representative value, if a value is repeated, that repeated value is translated into another correlated 1/1000%. So if the number 2 appears 2 times out of 1000 entries. The overall % of that value is .2% However, whatever value that was before it, say 1 , that appeared once would be .1%. so 1 - Maps to .1% and 2 - Maps to .3% (i.e. 2 counts of 2 out of 1000 = .2%, add prior values for cumulative distribution function. aka .1% and .2% = .3%) I also have a formula that flags for skewed distributions by testing the ECDF converted mean of a distribution for a max error of .275 from .5 mean. If it is, we do a frequency check on the dataset for the median value, and if the frequency is >50% of the values, then we treat the distribution differently. We don't want the 0% to be valued at a high value. Instead we remove all 0's from the list, convert to ECDF rating. Then re add back 0's. That way we normalize the tail end of the skew. So would this work to say compare varying sets of data to each other on an additive scale by normalizing them in this method?


title: Selecting number of clustering classes automatically

text: I am working in text clustering. I would like to find a way to identify the number of classes for the clustering process automatically rather than proving the number of class manually. Is their any proof or model for this?


title: How do you build dynamic linear models based on very frequent time series data?

text: I like to be able to use dynamic linear modelling (DLM) on a time series data to perform forecast. I am currently using lm function but I like explore dlm to see the difference. My data looks like this:                   structure(list(DATE = structure(c(1389287553, 1389287583, 1389287613,      1389287643, 1389287673, 1389287703, 1389287733, 1389287763, 1389287793,      1389287823, 1389287853, 1389287883, 1389287913, 1389287943, 1389287973,      1389288003, 1389288033, 1389288063, 1389288093, 1389288123, 1389288153,      1389288183, 1389288213, 1389288243, 1389288273, 1389288303, 1389288333,      1389288363, 1389288393, 1389288423, 1389288453, 1389288483, 1389288513,      1389288543, 1389288573, 1389288603, 1389288633, 1389288663, 1389288693,      1389288723, 1389288753, 1389288783, 1389288813, 1389288843, 1389288873,      1389288903, 1389288933, 1389288963, 1389288993, 1389289023, 1389289053,      1389289083, 1389289113, 1389289143, 1389289173, 1389289203, 1389289233,      1389289264, 1389289294, 1389289324, 1389289354, 1389289384, 1389289414,      1389289444, 1389289474, 1389289504, 1389289534, 1389289564, 1389289594,      1389289624, 1389289654, 1389289684, 1389289714, 1389289744, 1389289774,      1389289804, 1389289834, 1389289864, 1389289894, 1389289924, 1389289954,      1389289984, 1389290014, 1389290044, 1389290074, 1389290104, 1389290134,      1389290165, 1389290195, 1389290225, 1389290255, 1389290285, 1389290315,      1389290345, 1389290375, 1389290405, 1389290435, 1389290465, 1389290495,      1389290525, 1389290555, 1389290586, 1389290616, 1389290646, 1389290676,      1389290706, 1389290736, 1389290766, 1389290796, 1389290826, 1389290856,      1389290886, 1389290916, 1389290946, 1389290976, 1389291006, 1389291036,      1389291067, 1389291097, 1389291127, 1389291157, 1389291187, 1389291217,      1389291247, 1389291277, 1389291307, 1389291337, 1389291367, 1389291397,      1389291427, 1389291457, 1389291487, 1389291517, 1389291547, 1389291577,      1389291607, 1389291637, 1389291667, 1389291697, 1389291727, 1389291757,      1389291787, 1389291817, 1389291847, 1389291877, 1389291907, 1389291937,      1389291967, 1389291997, 1389292027, 1389292057, 1389292087, 1389292117,      1389292147, 1389292177, 1389292207, 1389292237, 1389292267, 1389292297,      1389292327, 1389292357, 1389292387, 1389292417, 1389292447, 1389292478,      1389292508, 1389292538, 1389292568, 1389292598, 1389292628, 1389292658,      1389292688, 1389292718, 1389292748, 1389292778, 1389292808, 1389292838,      1389292868, 1389292898, 1389292928, 1389292958, 1389292988, 1389293018,      1389293048, 1389293078, 1389293108, 1389293138, 1389293168, 1389293198,      1389293228, 1389293258, 1389293288, 1389293318, 1389293348, 1389293378,      1389293408, 1389293438, 1389293468, 1389293498, 1389293528), class = c("POSIXct",      "POSIXt"), tzone = ""), pc = c(7.29, 7.4, 7.15, 7.63, 7.46, 7.35,      8.1, 8.05, 7.9, 7.84, 7.86, 7.76, 7.63, 7.52, 7.21, 7.58, 7.98,      7.99, 7.16, 7.28, 7.07, 7.37, 6.93, 7.16, 7.14, 7.18, 6.97, 7.71,      7.19, 7.16, 7.29, 7.11, 7.08, 7.24, 7.31, 7.2, 8.23, 9.18, 8.13,      7.6, 7.53, 7.49, 7.18, 7.76, 7.03, 7.23, 7.23, 7.5, 7.18, 6.96,      7.08, 7.03, 7.2, 7.11, 7, 7.07, 7.92, 8.56, 8.05, 7.29, 7.21,      7.06, 6.94, 7.19, 7.02, 7.24, 7.76, 7.2, 7.32, 7.58, 7.68, 7.45,      7.52, 7.21, 7.22, 7.16, 7.57, 7.25, 7.32, 7.67, 7.23, 7.09, 7.12,      7.01, 7.08, 7.37, 7.25, 7.07, 7.26, 7.34, 7.01, 7.02, 7.15, 6.83,      7.1, 7.26, 7.58, 7.69, 7.75, 8, 7.8, 7.94, 7.64, 7.07, 7.13,      7.19, 7.39, 7.15, 7.41, 7.1, 7.11, 7.36, 7.04, 7.14, 7.04, 7.33,      8.37, 9.05, 7.85, 7.33, 7.62, 7.88, 7.29, 7.29, 7.64, 7.65, 8.18,      8.34, 7.67, 7.74, 7.49, 8.11, 7.55, 7.77, 7.57, 7.26, 7.31, 7.77,      6.92, 6.98, 7.09, 7.56, 7.16, 6.99, 6.92, 7.33, 7.03, 7.34, 7.27,      7.16, 7, 7.06, 6.92, 7.1, 6.83, 7.03, 7.64, 8.59, 7.27, 7.39,      7.21, 7.62, 7.34, 7.96, 7.03, 7.13, 7.06, 7.11, 7.05, 7.05, 7.11,      6.81, 6.78, 7.02, 6.82, 7.03, 7.79, 8.26, 7.25, 7.06, 7.1, 6.93,      6.96, 7.01, 7.24, 7.03, 7.42, 7.09, 7.31, 7.14, 7.15, 6.93, 7.52,      7.51, 7.16, 7.09, 7.45, 7.12, 7.08, 7.06)), .Names = c("DATE",      "pc"), row.names = c(1L, 3L, 5L, 7L, 9L, 11L, 13L, 15L, 17L,      19L, 21L, 23L, 25L, 27L, 29L, 31L, 33L, 35L, 37L, 39L, 41L, 43L,      45L, 47L, 49L, 51L, 53L, 55L, 57L, 59L, 61L, 63L, 65L, 67L, 69L,      71L, 73L, 75L, 77L, 79L, 81L, 83L, 85L, 87L, 89L, 91L, 93L, 95L,      97L, 99L, 101L, 103L, 105L, 107L, 109L, 111L, 113L, 115L, 117L,      119L, 121L, 123L, 125L, 127L, 129L, 131L, 133L, 135L, 137L, 139L,      141L, 143L, 145L, 147L, 149L, 151L, 153L, 155L, 157L, 159L, 161L,      163L, 165L, 167L, 169L, 171L, 173L, 175L, 177L, 179L, 181L, 183L,      185L, 187L, 189L, 191L, 193L, 195L, 197L, 199L, 201L, 203L, 205L,      207L, 209L, 211L, 213L, 215L, 217L, 219L, 221L, 223L, 225L, 227L,      229L, 231L, 233L, 235L, 237L, 239L, 241L, 243L, 245L, 247L, 249L,      251L, 253L, 255L, 257L, 259L, 261L, 263L, 265L, 267L, 269L, 271L,      273L, 275L, 277L, 279L, 281L, 283L, 285L, 287L, 289L, 291L, 293L,      295L, 297L, 299L, 301L, 303L, 305L, 307L, 309L, 311L, 313L, 315L,      317L, 319L, 321L, 323L, 325L, 327L, 329L, 331L, 333L, 335L, 337L,      339L, 341L, 343L, 345L, 347L, 349L, 351L, 353L, 355L, 357L, 359L,      361L, 363L, 365L, 367L, 369L, 371L, 373L, 375L, 377L, 379L, 381L,      383L, 385L, 387L, 389L, 391L, 393L, 395L, 397L, 399L), class = "data.frame")      Can somebody start me out about how would you go about building a DLM based on this type of time series data? Any insight/help/guide is greatly appreciated.


title: How to Obtain “Right” Parameters of Multinomial Logit Model (or Other Conditional Models) in R?

text: I started to use the function `multinom` of `R` package `nnet` in order to fit several conditional probability distributions with the multinomial logistic model. I need the parameters of the fittings in order to pass them to a Java program, which will compute the probabilities and use them. My problem is that the probabilities computed with the parameters returned by `multinom`, following the usual definition of multinomial logistic model, are not the same as those directly computed in `R`, which are the correct ones. On Stack Overflow I have already asked a question about this issue, but I do not still know how the `R` function `multinom` computes these probabilities; my guess is that it relies on neural networks, since this function belongs to `R` package `nnet`, but I do not have any idea about the details, and an inspection of the code led to nowhere. Do you know an `R` package which fits conditional probabilities and returns the corresponding parameters of the model, so that we may easily compute the probabilities in another program? E.g., using the MARS model (`R` package `earth`) or Projection Pursuit Regression (`R` package `ppr`) is not feasible, since computing the probabilities from the parameters of these models would be a mess. Besides, the function `mlogit` from the `R` package with the same name is not applicable as well, since the dataset should be in a certain format (we would also need the predictors corresponding to the alternative, "non-chosen" response variable).


title: How to choose the proper supervised machine learning algorithm

text: I have following data: ultra-sound audio recording, and I have the raw data of those samples, and of course the Fourier Transform data for those samples. Each dataset is assigned a value (for example 7.2 units), and I need some sort of algorithm which chooses which frequencies to use to get to that value. I have about 100 data samples, and I am looking at frequencies from 1 MHz to 80 MHz (in 100 kHz steps, so in total 790 features I guess). Why would I choose what algorithm?


title: How to specify the validation/training and test set in nntool in matlab

text: Is there any way to specify manually, the validation/test and training set in nntool in matlab. Actually, I wanted to test the performance of nn with increasing number of hidden neurons. So what I do is call this function               net = newfit(train_data', targets', n);      everytime with increasing number of hidden neurons. As a result, I guess each time it has different training/validation and test set. What I expected was the training error to decrease with increasing number of hidden layer neurons. But what I observed what initially the training error started increasing but after like 5 neurons it started increasing which is strange. Any suggestions?


title: Power calculation for case-control study

text: My study is assessing whether there is a correlation between two medical conditions (endometriosis) and bleeding disorders. The case participants have endometriosis and the controls do not. In both case and control we are measuring a number of different haemostatic variables (clotting factor levels, platelet number, platelet function) which are continuous numerical values usually from between 0-150. We are measuring roughly 8 different variables. I want to be able to demonstrate a signifant difference in the mean and range between cases and controls. Do you know, roughly what the sample size estimation would be? Many thanks, Joanna


title: Detecting and extracting end points of signal in larger signal

text: I have the following time series current measurement recordings at 5,000 samples per seconds for 10 minutes: https://docs.google.com/file/d/0ByIzsQrgb6ILakdaVVJkcEFUUGc/edit?usp=sharing I am trying to find a way to extract the recordings of the periodic intervals that are higher in current consumption as shown in the picture below: https://docs.google.com/file/d/0ByIzsQrgb6ILNXRJemtSTkhKRlk/edit?usp=sharing These periodic portions come every 120ish seconds. However, their interval duration is not constant, i.e. you can notice how the width changes. So I want to get the starting and end points to be able to extract that portion and do some analysis on it. Here are three intervals extracted manually and plotted: docs.google.com/file/d/0ByIzsQrgb6ILd241d2poMjc2SFU/edit?usp=sharing I have been trying various things attempting to extract this portion but wasn't successful. Maybe I am following the wrong methods overall. I am doing this in MATLAB and I am not a DSP expert. The latest thing I have tried to do is a simple moving average with a lag of 5,000, i.e. the sampling rate, using the MATLAB function `tsmovave(data, 's', 5000)` and I got the following: docs.google.com/file/d/0ByIzsQrgb6ILT0dpbWpIWDE0MHc/edit?usp=sharing However, the signal doesn't look clean enough to apply simple thresholds on it and pull out what I need. I also tried to increase the lag time to 120 seconds (i.e. the period interval) and then attempted to extract minimums and maximums using the MATLAB function `findpeaks()`, but that didn't work because there are a lot of local minimums and maximums in the data. So I tried to resample at a lower rate but that still didn't work. Someone suggested that I should use matched filter and provide a manually extracted interval to the filter. Another method I came across in my research is implementing some sort of step-finding algorithm. From your experience, what's the best way to go about doing this DSP. It seems like a simple problem but I spent a lot of time trying to solve it and couldn't. I am not an expert in DSP. I appreciate any feedback and direction.


title: Method to remove bad values in time series (bad values known to take on a particular value)

text: This sounds easy, but I don't know of a good statistical method for it. I have a time series that has (good) data points that range from ~3.5 to 30. The data are collected by an automated sensor. However, there are flawed measurements in the time series-- the sensor will sometimes read values that are typically at exactly 4.7. Example (in R):               Data <- c(10.7,  4.7, 10.7,  4.7, 11.5,  4.7,  8.4,  4.7, 11.0,  4.7, 10.8,  4.7, 12.0,  10.3,  4.7, 10.4, 10.9,  4.7, 10.8,  4.7, 11.1,  4.7,  8.8, 4.7,  6.4,  4.7,  9.1,  4.7,  9.1,  4.7,  9.7,  4.7,  7.8,  4.7,  7.2,  4.7,  6.2,  4.7,  8.7,  4.7,  9.0,  4.7,  9.6,  4.7,  7.9,  4.7,  9.2,4.7,  8.1,  4.7,  7.4,  4.7,  7.6,  9.3,  4.7,  9.2,  4.7,  8.8,  4.7,  9.3,  4.7,  9.2,  4.7,  7.0,  4.7,  9.2,  4.7,  8.5,  4.7,  6.2,4.7,  7.1,  4.7,   7.4,  4.7,  8.0,  4.7,  7.3,  4.7,  6.6,  4.7,  6.9,  4.7,  7.2,  7.9,  4.7,  9.0,  4.7,  8.7,  4.7,  8.2,  4.7,  5.1,4.7,  5.6,  4.7,  7.0,  4.7,  9.4,  4.7,  7.6,  4.7,  8.6,  4.7,  9.3,  4.7,  9.7,  4.7, 10.4, 4.7,  10.6,  4.7, 10.9,  4.7, 10.2,  4.7, 10.0,4.7,  8.3,  10.0, 8.7,  4.7, 10.2,  9.7,  4.7, 10.2,  4.7, 10.6,  4.7, 10.5,  4.7,  9.7,  4.7, 10.6,  4.7, 10.6,  4.7, 11.6,4.7, 11.4,  4.7, 10.2,  4.7, 10.3,  4.7, 10.0,  4.7,  9.3,  4.7,  9.1,  4.7, 10.2,  4.7,  8.5,  4.7,  7.2,  4.7, 10.4,  4.7, 10.4, 10.0,4.7,  9.9,  4.7, 11.2, 9.5,  4.7, 11.5,  4.7, 11.5,  4.7, 11.1,  4.7, 11.4,  4.7, 12.0,  4.7, 11.4,  4.7, 11.2,  4.7,  9.9,  4.7, 11.6, 4.7, 14.4,  4.7, 11.5,  4.7, 11.1,  4.7, 11.3,  4.7, 10.9,  4.7, 11.1,  4.7, 11.1,  4.7, 10.9, 11.4,  4.7, 12.6,  4.7, 11.0,  11.3)          plot(1:length(test), test, type="o")      As you can see, almost every other point is 4.7. I don't want to simply use the rule of "throw out every value that is equal to 4.7" for two reasons: 1) sometimes there are real values that approach 4.7, which can become apparent as you see the slope of the time series change and enter the 4.7 region; 2) the "bad" values aren't always _exactly_ 4.7 (they almost always are, but there is no guarantee). This seems like a fairly simple problem, and I'm hoping that someone has heard of a statistical approach that would be suitiable for this problem. The desired outcome would be a time series where the "bad" values were statistically flagged (e.g., with a probability of being an error). Suggestions? Thanks! EDIT: I'm still mulling over the first two responses. My original inclination was to define two states of the system: an erroneous state and a measured state. The measured state would reflect a modeled process (e.g., AR). The erroneous state would reflect my knowledge that bad values tend to be near 4.7. Both states could ~N(), with appropriate means and sd's. At each time step I could calculate the likelihood of either state, and define a threshold ratio for declaring an obs erroneous. This is similar to the 2 answers. For this approach, I'm stuck at finding a way to let the model for the measured state make estimates w/o being influenced by bad values. However, as stated above, I'm still considering the other two approaches.


title: OLS vs. logistic regression for exploratory analysis with a binary outcome

text: In the idealized logistic model, we obtain an S-shaped curve linking each continuous IV to the DV. But in practice this infrequently happens, making the logistic approach seem a little less superior for such types of data. Of course predicted probabilities that each observation will be "1" on the DV are usable in logistic and not in OLS regression, since in the latter these probabilities can exceed the bounds of [0,1]. But, for exploratory purposes, how sound is it to use OLS to see which IV have strong vs. moderate vs. weak relationships with the DV? Wouldn't this amount to a sort of multivariate version of point-biserial correlation? (Standardized regression coefficients, not to mention collinearity statistics and partial plots, are all I think more easily obtained in OLS than in logistic.)


title: How to compare whether the coefficients of two independent variables statistically different from each other?

text: If I have two independent variables and they are dummy variable along with other independent variables and I run a linear probability model, I want to compare whether the coefficients of two dummy variables are statistically different from each other. I do not know how to compare that. Could someone help me on that?


title: Do Beta weights from regression have error terms?

text: I am looking at standardized regression weights (i.e., Beta weights). I was thinking of reporting the errors next to the weights in a figure, but upon some thought I was debating whether such errors even exist. I know un- standardized regression weights have standard error parameters, but do standardized regression weights have error terms?


title: How to extract uniform noise from the mixture distribution

text: Let $X$ be a sample from the distribution which is the mixture of the useful signal with the distribution $\xi(\theta)$ and a uniform noise $U[a, b]$. The probability of observing $U[a, b]$, the parmeters $a$ and $b$ are not known. I would like to estimate $\theta$. How do I do it? By the mixture I mean weighted average of probability density functions


title: A failure of 'regressing out'?

text: I sample two timecourses:   * $X=shardNoise+noise_X$   * $Y=signal+\beta*sharedNoise+noise_Y$ ($signal$, $sharedNoise$, $noise_X$ and $noise_Y$ are independent from each other) My goal is to 'regress-out' $sharedNoise$ from $y$, so my estimation of $signal$ won't be affected by it. A straightforward approach would be to use a linear regression: $\hat{y}=\hat{\beta}x$ and then 'regress-out' by calculating $y-\hat{y}$. Unfortunately, the outcome of this approach is poor - there are clear 'leftovers' of $sharedNoise$ in $y-\hat{y}$. I think that this underfitting is caused by the presence of $noise_X$. For example, consider increasing the amplitude of $noise_X$ to infinity - that would drive $\hat{\beta}$ to zero.   * Is there a better way to solve this problem? I guess that many will think of ICA (Independent Component Analysis) but I prefer a supervised approach in which I won't need to guess which component stands for $sharedNoise$.   * Can such supervised approach extended to a case in which there are multiple Xs and Ys? (all consisting of different combinations of $sharedNoise$ plus some local signal and/or noise).


title: What is parameter fine tuning means in SVM?

text: I got this sentence in one of paper, but I dont understand what does it mean?? "Training a learningbased classifier such as an SVM on an imbalanced dataset often requires parameter fine-tuning, making the system less generic "


title: Obtain the covariance matrix with the correlations and the variances

text: i have the following problem that makes me crazy. COnsider a trivariate normal variable, let $\sigma_1^2$, $\sigma_2^2$ and $\sigma_3^2$ be the variance of the three components and $\rho_{1,2}$, $\rho_{1,3}$ and $\rho_{2,3}$ be the correlation coefficients between the components. I want to obtain the covariance matrix of the trivarate normal variable. For me it is a really simple task but there is something wrong.... The covariance matrix is $$ \left( \begin{array}{ccc} \sigma_1^2 & \sigma_1 \sigma_2 \rho_{1,2} & \sigma_1 \sigma_3 \rho_{1,3} \\\ \sigma_1 \sigma_2 \rho_{1,2} & \sigma_2^2 & \sigma_2 \sigma_3 \rho_{2,3} \\\ \sigma_1 \sigma_3 \rho_{1,3} & \sigma_2 \sigma_3 \rho_{2,3} & \sigma_3^2 \end{array} \right) $$ but if i use the following value $\sigma_1^2=1$, $\sigma_2^2=1$, $\sigma_3^2=3$, $\rho_{1,2}=0.9$, $\rho_{1,3}=-0.9$ and $\rho_{2,3}=0.1$, the covariance matrix $$ \left( \begin{array}{ccc} 1 & 0.9 & -1.5588457 \\\ 0.9 & 1 & 0.1732051\\\ -1.558846 &0.1732051 & 3 \end{array} \right) $$ is not positive definite. It seem strange to me...


title: Can the curse of dimensionality be solved by changing the metric tensor?

text: My (not too deep) understanding of the curse of dimensionality that affects a classification algorithm, such as k-Nearest Neighbor, is that at higher dimensions the 'sparsity' of euclidean space kicks in (this can be seen, for example by comparing the volume/content/measure of the unit ball with respect to the unit box) I wanted to know if researchers have considered working on riemannian manifolds other than euclidean space (or spaces with other $L_p$ norms ($p \neq 2$)$^1$), say with metric tensor $\sum g_{ij} dx^i \otimes dx^j$ where the coefficients $g_{ij}$ are non-constant or maybe even dependent on the data? 1.Aggarwal C. C., Hinneburg A., Keim, D. A. (2001), ”On the Surprising Behavior of Distance Metrics in High Dimensional Space”


title: Distance Transform In Image?

text: I understand how to compute Distance Transform for BW image In BW image we can calculate the distance of the pixel with non zero pixel value using Eucledian Distance. But I dont undertand how to calculate it for GrayScale and RGB (Colour) Image?? And what is purpose by calculating the DT?? What we can use with information of DT??


title: Tsdiag results analysis help

text: What do you think about result from my arima model? Should I need to change the model. I do not want to change it because it is the best I have but p values do not seem fine ![tsdiag](http://i.stack.imgur.com/VqrhR.png)


title: Interraction in Poisson regression

text: I'm trying to understand what is happening in my Poisson regression. I have the following variables: Dependent variable: transactions(number of transactions(continuous)). Independent variables: education(years of education (continuous)),sex(binary),health(health status), income(continuous) other than salary, loan(binary(yes no)). i create the following model(model 1 see below). What i know: there is over dispersion and 2 variables are not significant at a 95% confidence interval. i know that there is no confounding! and very low correlation between terms i create a second model with an interaction term. And suddenly everything becomes very significant, except sex. The estimates change. Is this because the interaction term causes confounding? or what would be the explanation? how could i check which model best fits the data? \---------MODEL 1----------------------------------------------------------------------               glm(formula = transactions ~ salary + education + factor(sex) +      factor(health) + income + factor(loan), family = poisson)     Deviance Residuals:          Min       1Q   Median       3Q      Max       -3.9632  -1.9344  -0.6717   0.5504  12.6591       Coefficients:                      Estimate Std. Error z value Pr(>|z|)         (Intercept)     -0.413667   0.269088  -1.537 0.124221         salary           0.004367   0.001303   3.352 0.000801 ***     education       -0.010717   0.006005  -1.785 0.074302 .       factor(sex)1     0.041353   0.027840   1.485 0.137450         factor(health)1  1.133173   0.030301  37.397  < 2e-16 ***     income           0.149229   0.036047   4.140 3.48e-05 ***     factor(loan)1   -0.140154   0.026550  -5.279 1.30e-07 ***     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1     (Dispersion parameter for poisson family taken to be 1)         Null deviance: 8848.8  on 2226  degrees of freedom     Residual deviance: 7419.9  on 2220  degrees of freedom     AIC: 11899     Number of Fisher Scoring iterations: 5      \-----------------------MODEL 2----------------------------               `Call:     glm(formula = transactions ~ salary + education + factor(sex) +          education:income + factor(health) + income + factor(loan),          family = poisson)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -3.9828  -1.9270  -0.6606   0.5475  12.6271            Coefficients:                       Estimate Std. Error z value Pr(>|z|)         (Intercept)      -6.663294   1.370830  -4.861 1.17e-06 ***     salary            0.004822   0.001305   3.694  0.00022 ***     education         0.548942   0.120117   4.570 4.88e-06 ***     factor(sex)1      0.043472   0.027839   1.562  0.11840         factor(health)1   1.130069   0.030289  37.310  < 2e-16 ***     income            0.951831   0.176518   5.392 6.96e-08 ***     factor(loan)1    -0.137008   0.026557  -5.159 2.48e-07 ***     education:income -0.071936   0.015443  -4.658 3.19e-06 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for poisson family taken to be 1)              Null deviance: 8848.8  on 2226  degrees of freedom     Residual deviance: 7397.5  on 2219  degrees of freedom     AIC: 11879          Number of Fisher Scoring iterations: 5}`


title: How to calculate probability distribution for the duration of recurring tournament?

text: I've hit an interesting general problem and cannot work it out. I would be grateful for any pointers: The problem: 5 players are randomly assigned a number x=1-5. A player compares his number with another random player. The player with the lower score leaves and the player with the higher score goes into the next round and repeats the tournament. What is the probability that player i goes out in round 1, 2, 3 or 4 Current strategy: Calculate the probability of winning with number x, multiply by the probability of getting that number and sum to find the probability of winning that round. For round 1, P(winning with x=1)=0, P(winning with x=2)=1/4, P(winning with x=3)=2/4 etc. The probability of getting any number x is 1/5, so the probability of succeeding in round 1 is (1/5 x 0)+(1/5 x 1/4)+(1/5 x 2/4)+(1/5 x 3/4)+(1/5 x 4/4) This sums to 0.5 so the probability of going out in round 1 is (1-0.5) = 0.5. For the next round, I have to generate the probabilities of having numbers 2,3,4,5 and this is where I am getting stuck. It will no longer be uniform as there will e.g be a higher probability that x=4 reached round 2 than x=2 but I cannot work this out.


title: how to identify salient words after LSA?

text: I used Latent Semantic Analysis (LSA) to extract latent topics (i.e., the polynomials coefficient1*word1 + coefficient2*word2 + ...) from a certain corpus. I know that the larger the (absolute value of the) coefficient, the more important the word for that topic. But is there a better metric? Python's gensim package uses this metric, but I couldn't find any textbook or paper that discusses the issue. Any pointers will be greatly appreciated!


title: Autocorrelation problem in multiple regression on SPSS

text: What is autocorrelation and how it can effect the results of multiple regression in _**SPSS**_. Durbin-Watson value of **1.04** is acceptable or not? If this value effects the results, then how to **correct** it in SPSS.


title: two questionnaires, two variables, two populations but one correlation

text: I have two variable in my research. the first is managers commitment & the second is personnel satisfaction. so we have two population & two questionnaire, one for managers & the other with "different" questions for personnel. moreover the number of managers is different from the number of personnel. how can I calculate correlation & regression between these two population as one population with two variable? thanks


title: Why do we stabilize variance?

text: I came across variance stabilizing transformation while reading Kaggle Essay Eval method. They use a variance stabilization transformation to transform kappa values before taking their mean and then transform them back. Even after reading the wiki on variance stabilizing transforms I can't understand, why do we actually stabilize variances? What benefit do we gain by this?


title: Understanding factors returned by RPart classification

text: This is a very basic question on using R for classification. I'm trying to use rpart for classification task and would like to have a class label as a result, i.e. I use type="class" in predict method               predict(tree, data, type = "class")      The class labels are "0" and "1". When I try to classify a single observation I get something like this:               1      0      Levels: 0 1      As I understand this is a factor with 2 levels: 0 and 1. However, why there are 2 numbers reported and what is exactly the class assigned to this observation?


title: How do I interpret logistic regression output for categorical variables when two categories are missing?

text: I am using binary logistic regression; the dependent variable is 1 or 0; the independent variables are two groups: the first group includes continuous variables (`LNTA: logarthim of total assets`, `ROA: return on assets`, and `Leverage`; the second group includes categorical variables (`Type of auditor: 1 or 0`, `Industry sector: 1,2,3, and 4`, `country: 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10`), and finally `region: 1 or 2`. The problem is that when I put all these independent variables together, I get results only for 8 countries not for 9 countries. I know that I have 10 countries and in the results table only 9 countries will appear as the country number 10 is reference category; but this case is different because both countries 9 and 10 are not included in the table.


title: What's the meaning of the class indicator matrix when transforming the class label matrix into it in canonical correlation analysis?

text: When using canonical correlation analysis (CCA), we can integrate the dataset and label information via transforming the class label matrix _Y_ into the class indicator matrix _T_. Such as: $T = (YY^T)^½Y$ in this article on LS CCA. While for the numeric dataset and the binary dataset, the binary dataset can be viewed as the class label matrix and transformed into the class indicator matrix. After that, should I use the CCA? Though without this transformation, the CCA can be used, I'd like to know how to explian this kind of transformation once using it on the binary dataset (not label matrix). Thank you.


title: Autocorrelation of discrete time series

text: I am currently planning on calculating the autocorrelation for various lags given a time series. However, my elements of the time series are "discrete" and abstract classes; i.e., no integers. For example, my series could look like: class A, class B, class A, class A ... Of course I could now say that class A = 1 and class B = 2 and then use the common Pearson correlation way for determining the autocorrelation. However, this seems to be problematic to me, because the value of the number does not express anything. Is there some way around this?


title: What is the preferred way to give asymmetric uncertainties?

text: I have seen something like $4.1^{+2.1}_{-1.5}$ or $4.1 [+2.1, -1.5]$ but I would like to know, what is the preferred way or something like a standard in scientific work (especially in physics/engineering)?


title: Obtain p-values from t-values

text: I used Tobit regression and wanted to test if the coefficients are significant. The problem is that when I run the R code I get t-values instead of z-values. I found code on how to calculate p-values for z-values and I'm wondering if you can use the same code for the t-values that you use for the z-values to calculate the p-values. **This is the code for when you have z-values:**               ctable <- coef(summary(m))     pvals <- 2 * pt(abs(ctable[, "z value"]), df.residual(m), lower.tail = FALSE)     cbind(ctable, pvals)      **This is the code I used for the t-values:**               ctable <- coef(summary(m))     pvals <- 2 * pt(abs(ctable[, "t value"]), df.residual(m), lower.tail = FALSE)     cbind(ctable, pvals)                              Value  Std. Error   t value         pvals     (Intercept):1   209.559678 32.50165463  6.447662  3.327237e-10     (Intercept):2    4.184758  0.05227107 80.058783 7.628269e-246         read             2.697963  0.62024382  4.349843  1.737355e-05        math             5.914589  0.70502334  8.389210  8.842612e-16     proggeneral     -12.714341 12.36481457 -1.028268  3.044547e-01     progvocational  -46.143265 13.65615960 -3.378934  8.002457e-04


title: Correlation coefficient of noise

text: A correlation coefficient is computed as $$r=Cov(X,Y)/(\sqrt{Var( X)}*\sqrt{Var (Y)})$$ What is the correlation coefficient of a number with itself then. I have a noise function,assuming with a known value and variance; and I perform a correlation of noise with itself what should I get as the correlation coefficient. The noise follows a $\chi^2$ distribution.


title: Fast missing data imputation in R for big data that is more sophisticated than simply imputing the means?

text: I need a package for missing data imputation in R. But since I am dealing with big data, the number of missing data entries can also be high. The packages which impute using mean or median are of course working fast, but more complicated packages which impute using regression or PCA take too long for a high number of missing values. I tried `missMDA` and `missForest`, but as I said, they look like taking forever. There is a package named `FastImputation`, but I could not figure out how to use it when I have no patterns from some training data. Any suggestions of packages which would impute fast?


title: Confidence interval for the "shift parameter" of a non-central exponential distribution

text: Suppose that $X_1, X_2, \ldots, X_n$ (with $n > 0$) is a random sample from a non-central exponential distribution with probability density function: $$f(x | \lambda, k) = \lambda * e^{-\lambda*(x-k)}$$ Both the scale parameter ($\lambda$) and the **_shift parameter_** ($k$) are unknown, with: $k < x < \infty$ and $\lambda > 0$ . For the previously mentioned sample ($X_1, X_2, ... , X_n$) and $k$, the maximum likelihood estimate (MLE) is: $$K = min_i X_i$$ **Question** : how would a $95\%$ confidence interval (for the shift parameter, $k$) look like ? I'd imagine it must be something like $[K-\epsilon, K]$, but how can we figure out the value for "$\epsilon$". Thank you !


title: Is it possible for a regression to have the right functional form but still suffer from omitted variable bias?

text: Suppose that $$ Y = b + aX + e$$ where you know that $E[Y|X] = b + aX$. Is it true that the model cannot suffer from omitted variable bias? If this is true, then it follows that omitted variable bias can always be considered as a funcional form bias (namely, changing the functional form $f(X)$ will remove the bias.) Here is a proof of why I think this is true. Since $b + aX = E[Y|X]$, we have that $$ E[e|X] = E[Y-b-aX|X]= E[Y - E[Y|X] |X] = E[Y|X] - E[E[Y|X] |X] = E[Y|X] - E[Y|X] = 0.$$ Therefore, $e$ and $X$ are not correlated and there is not omitted variable bias.


title: How to estimate GARCH in R? (Exogenous variables in mean equation)

text: What I'm trying to do is estimate the following GARCH(1,1) model in `R` with the garchFit function from the `fGarch` package: * * * Mean equation: $Y_t = a + bX_t + e_t$ $e_t = z_th_t^{0.5}$ $e_t \sim N(0,h_t)$, $z_t \sim N(0,1)$ Variance equation: $h_t = \omega + k_0e^2_{t-1} + k_1h_{t-1}$ * * * In the above, $Y_t$ is the response of my mean equation and $X_t$ is the predictor. I will estimate this equation within the GARCH framework because of heteroscedasticity of residuals. I want estimates of both the mean equation and the variance equation (similar to what `EViews` would give). However, the garchFit function has two inputs that I'm concerned about:   * _formula - ARMA(m,n) + GARCH/APARCH(p,q) mean and variance specification_ As you can see in my equations, my mean equation is _not_ of the ARMA(m,n) form.   * data - any univariate time series which can be converted into a timeSeries using the generic function as.timeSeries I'm modeling both $Y_t$ and $X_t$ together; thus it shouldn't be univariate. If I put this model through `EViews`, It'll give me parameter estimates of $a$ and $b$. This is what I want, but in `R`. In `R`, I can do `garchFit(data=Y)` (where `Y` is the column of $Y_t$), and it'll do a GARCH(1,1) over $Y_t = \mu(Y) + e_t$; but it doesn't give me $a$ or $b$ in the mean equation that I want!


title: Rules of algebra with Normal distribution

text: As part of my self study for Financial Economics, I came across simplifying the following expression, $$\sigma Z_{0.25} + (2\sigma)(Z_{0.5} - Z_{0.25})$$ where $Z$ are independent Brownian motion ($Z_0 = 0$ and $Z_t-Z_s \mathtt{\sim} N(0, t-s)$. It is simplified as, $$\sigma N(0, .25) + 2\sigma N(0, .25) = N(0, .25\sigma^2) + N(0, .25 \times 4\sigma^2)$$ $$= N(0, 1.25\sigma^2) = \sqrt{2.5}\sigma Z_{0.5}$$ However, if I _simplify_ the first expression to $\sigma (2Z_{0.5} - Z_{0.25})$, I get a different Normal distribution, $N(0, 1.75\sigma^2)$ which is $\sqrt{3.5}\sigma Z_{0.5}$. **Edit** : Here's how I did the above simplification, $$\sigma Z_{0.25} + (2\sigma)(Z_{0.5} - Z_{0.25})$$ $$= \sigma Z_{0.25} + 2\sigma Z_{0.5} - 2\sigma Z_{0.25}$$ $$= 2\sigma Z_{0.5} - \sigma Z_{0.25}$$ $$= \sigma (2 Z_{0.5} - Z_{0.25})$$ using basic algebra. So, obviously, the Brownian motion variables and Normal distributions don't support all of algebra. So, what subset of the algebraic rules are allowed when dealing with these? Is there somewhere these are set out that I can refer to? Thank you.


title: How to find the correlation coefficient between two technologies when those are sub correlated?

text: Suppose I have two power generating technologies, Coal and Oil, of which the generating cost components(total generating cost = Capital cost+Fuel cost+Variable O&M cost+Fixed O&M cost) are correlated in the following manner:                                                        Oil                           Fuel           Variable O&M             Fixed O&M                Fuel       0.48              0                       0     Coal    Var O&M       0                 0.7                     0.1           Fixed O&M       0                 0.1                     0.7      Suppose Technology A is COAL and Technology B is OIL. How do I find the correlation between the total costs of the two technologies, that is COAL and OIL (one correlation coefficient) when three categories of costs are already correlated as given above? I have the weights of cost components of each total technology cost and the standard deviation of each cost component of the two technologies as well.


title: Most natural generalization of covariance/correlation to model dependence of extreme events

text: One of the most serious shortcomings of covariance/correlation are the assumptions of linearity and normality. What is the most natural generalization of these measures of dependence when you want to model the dependence structure of extreme events using heavy- tailed distributions, e.g. the Generalized extreme value distribution? With "most natural generalization" I mean that the classical covariance/correlation is included as a special case when the usual assumptions hold. (Disclosure: After having received no answers for nearly two weeks I posted this question also at Quantitative Finance)


title: How to get pooled p-values on tests done in multiple imputed datasets?

text: Using Amelia in R, I obtained multiple imputed datasets. After that, I performed a repeated measures test in SPSS. Now, I want to pool test results. I know that I can use Rubin's rules (implemented through any multiple imputation package in R) to pool means and standard errors, but how do I pool p-values? Is it possible? Is there a function in R to do so? Thanks in advance.


title: comparing proportions for non-mutually exclusive data

text: I have an experiment where I am interested in comparing proportions. Basically, participants write 4 different essays to 4 different open ended questions. Each answer is rated on a categorical variable of 1, 2, and 3. However, it is possible to write an essay containing more than 1 category. I don't need to predict, I want to look at proportions. All I want to be able to do is to test if there is a greater proportion of response 3's than 1 or 2. If they were not mutually exclusive I think I could just do a couple z-tests for proportions in dependent groups. But because the categories are not mutually exclusive I don't think I can. Very confused at this point. Suggestions on type of analysis?


title: The null and alternate hypotheses

text: The null and alternate hypotheses are: $H0 : μ_d ≤ 0$   $H1 : μ_d > 0$ The following sample information shows the number of defective units produced on the day shift and the afternoon shift for a sample of four days last month.                     Day        1       2       3       4     Day shift       10       12      15      18        Afternoon shift  8       9       13      16         At the .005 significance level, can we conclude there are more defects produced on the afternoon shift? Hint: For the calculations, assume the Day shift as the first sample. (a) State the decision rule. (Round your answer to 2 decimal places.) Reject H0 if t > (b) Compute the value of the test statistic. (Round your answer to 3 decimal places.) Value of the test statistic


title: Forecasting nonstationary time series

text: I would like to forecast the non-stationary time series, involving several crucial a-priori assumptions following from studying of instances of such series.   1. I've constructed time-averaged one-point probability distribution function approximated by normal distribution. $$\hat p(x) = \frac{1}{\sqrt{2\pi \sigma^2_{\infty}}} \exp\left(-\frac{x^2}{2\sigma^2_{\infty}}\right)$$ From this point of view, I want the forecast $z_t(l)$ not to exceed this when $l \to \infty$. To put it in other words, variance of $z_t(l)$ must be bounded.   2. ~~The average two-point probability distribution function $\hat p(x_i,i;x_j,j)$ also has been constructed, which led to identification of autocorrelation function. $\rho(j) \approx A j^{-\alpha} $ provided $0 <\alpha<0.5$.~~ At first, Box-Jenkins identification process led me to $ARIMA(0,1,3)$ model, however   1. I can't have bounded variance until $d \ne 0$ (which follows from equations for BJ weights $\psi_j$). At the same time, I can't use $d=0$ since initial autocorrelation decreases slowly (which is probably evidence of non-stationarity according BJ). This is the main obstacle to me.   2. Visually, simulation of $ARIMA(0,1,3)$ does not coincide with behaviour of my samples. And correlations of first difference of the series are in the bad agreement with correlations following from the model.   3. The analysis of residuals shows significant correlations starting lag 3. This is why my initial statement about $ARIMA(0,1,3)$ is incorrect. Trying to fit different $ARIMA(p,0,0)$ models, I see that there is significant residual correlations close to the lag $p$ for every $p$. It may assume that I need $ARIMA(\infty,0,q)$ model (as limiting choice), for instance fractional ARIMA. From [1] I've learned about Fractional $ARIMA(p,d,q)$ models which are $ARIMA(\infty,0,q)$ in effect.   1. I've not found any GNU R packages with support of missing values for this. Missing values seems to be a kind of challenge.   2. The publications on fractional ARIMA are quite rare. Are such fractional models really used? Maybe there is a good replacement of ARIMA models for my needs? The forecasting is not my major, I have only pragmatic interest.   3. From different literature (for instance [2]), I learned that it is practically impossible to decide between fractional ARIMA and models with "level shift". However, I have not found the package for GNU R to fit 'level shift' models. [1]: Granger, Joyeux.: J. of time series anal. vol. 1 no. 1 1980, p.15 [2]: Grassi, de Magistris.: "When long memory meets the Kalman filter: A comparative study", Computational Statistics and Data Analysis, 2012, in press. **Update:** to render my own progress and to answer @IrishStat My statement about two-point probability distribution is incorrect in general. Constructed in this way function will depend on full series length. So, there is a little to extract from this. At least, parameter named $\alpha$ will depend on full series length. Lists 2 and 3 also have been updated. My data is available as dat file here. At the current moment, I doubt between FARIMA and level shifts, and I still can't find appropriate software to check this options. This is also my first experience with model identification, so any help will be appreciated.


title: Replicating simulation results from a paper

text: I’m reading R. Kreps' paper _Parameter uncertainty in (log)normal distributions_ and trying to figure out how the simulations were done. In order to generate Figure 1, Eqn (2.41) was used. So this is what I did in R: This is the equation they used: $Z_{eff} = v + z\sqrt{\frac{n(1+v^2)}{w}}$ where Note: This code is for the case when $n=3$               v <- rt(10000000, 1)/sqrt(3-2)     w <- rchisq(10000000,2)     z <- rnorm(n=10000000, m=0, sd=1)      z_eff_3 <- v + z * sqrt((3*(1+v*v))/w)     plot(density(z_eff_3),ylim=c(0,1))      However, I got a completely different graph. I would truly appreciate it if someone can take a look and explain to me what I did wrong.


title: EM algorithm R code on Cox PH model with frailty

text: Let say I have a 'kidney catheter' data set. Data are about the recurrence times to infection, at the point of insertion of the catheter, for kidney patients using portable dialysis equipment. Catheters may be removed for reasons other than infection, in which case the observation is censored. Each patient has exactly 2 observations. If I want to fit a Cox PH model with random effect (here Gamma frailty) using EM algorithm. By using built in survival package R code `coxph()`, I can easily do it like this > >     library(survival) >     data(kidney) >     fit<-coxph(Surv(time, status) ~ age + sex + >                frailty(id, dist='gamma', method='em'), kidney) >   But if I want to write a step by step function for EM algorithm, how can I proceed?


title: Relative variable importance with AIC

text: I am confused and just need some confirmation about calculating the relative variable importance value for the co-variates I used in AIC model selection procedures. I know that there is this one discussion but it doesn't quite confirm explicitly enough what I should do. Burnham and Anderson (2002) describe a simple way to quantify variable importance. > Page 168: Estimates of the relative importance of predictor variables xj can > best be made by summing the AIC weights across all the models in the set > where variable j occurs. However, to use this method, one must have an equal number of models for each variable; otherwise, some variables will be over represented or under represented resulting in biased relative importance values. > Page 169: When assessing the relative importance of variables using sums of > the AIC weights, it is important to achieve a balance in the number of > models that contain each variable j. Does this mean if I have a set of models with their model weights from AIC procedure (these are not ranked by weight, just the order I created them):               1   INTERCEPT                 2   REPRO   TIME              3   REPRO   TIME    R*T       4   REPRO   TIME    WR       5   REPRO   TIME    WR  WR*R      6   REPRO   TIME    WR  WR*T      7   WR      To calculate the relative variable weight I would sum up the weight for each incident that TIME was in the models and I would do so for each of the other variables. However, this is not completely correct right? Because there is not a balance in the number of models that contain each variable right? So, to correct for this I would then divide the sum of these weights by the number of models that had that variable. (Kittle et al 2008 "the scale-dependent impact of wolf predation risk...." does this). So, for instance if the sum of the weights for time was .75 I would divide it by 5 because it was in 5 of the models, likewise, WR would be divided by 4. It seems like a silly question but it really changes the results and interpretation from my analysis. Because for instance, WR*T is only in 1 model and it comes out to be in one of the top models so it has a high model weight, but Time and Repro are also in this top model but also in 4 other candidate models. So dividing the weight of T & R by 5 reduces the importance of T or R from (0.999), giving them a RVI of 0.2 and the RVI to WR*T value of 0.7. Is that right? In addition to this, my next question would be - do you do this over JUST the "BEST" (within 2AIC or what ever criteria) models or over all 7 regardless of what surfaced to the top? I used MuMIn package and use the importance command, but then when you use the get best models, it asks if you want to recalculate the importance which then it recalculates for just the top models. Which is more appropriate to use? This doesn't make sense when only 1 model is the best. I would then assume it should be calculated over all models.


title: probability distribution of output value with regression tree methods

text: If I have a regression problem where I try to estimate the value of $y$ as function of $x_1 \dots x_d$: $$ y = f(x_1,\dots,x_d) $$ using a Boosted Regression Tree or a Random Forest Regression, is it possible to estimate the probability distribution of the output value (instead of only an estimate of $y$) : $$ p(y\ |\ x_1,\dots,x_d) $$ from the distribution of the outputs of the individual nodes ? (i.e. for fixed $x_1,\dots,x_d$ I'm interested in having a probability distribution for $y$ from which one then could obtain the mode, mean and also the width which is an indication of 'how well' this value is predicted).


title: Estimating error from repeated measurements

text: Say we use a ruler to make a measurement of the width of a block of wood. We get some value like `4.2 +/- 0.1 cm`, where the error is our estimated error of our ruler's precision. Now we have four new individuals measure the same block of wood, and two of them obtain `4.5 +/- 0.1cm` and two obtain `4.2 +/- 0.1 cm`. If we take the mean of these five measurement we obtain `4.32 +/- 0.16 cm`, where I'm ignoring significant figures in this example and the uncertainty is the standard deviation of the samples. Now, in stating the total uncertainty for the width of the block of wood, should I combine the 0.16 in quadrature with the 0.1, or should I propagate the error of each measurement through the standard deviation formula and obtain an uncertainty on my uncertainty, like `0.16 +/- alpha`, where alpha is the propagated uncertainty value? Or perhaps I should ignore the standard deviation of the measurements and use the classic 0.1/sqrt(5)---but this is silly given the spread of our measurements. This seems like a simple question but I can't seem to find two sources that agree on the matter (har har). Thanks for any help.


title: Is it possible to calculate variable confidence intervals, conditional on $\hat{Y}$ to address heteroscedasticity?

text: Estimating confidence intervals for non-normally distributed residuals can be accomplished using bootstrapping procedures, sandwich estimators or quantile regression. But is there a way to calculate $\beta$ confidence intervals given the estimated value? That is, for each $\hat{Y}$ we know the residuals variance for this predicted value _region_ , so rather than applying the same (robust) intervals to all predicted values, isn't it possible to adjust the intervals conditional to $\hat{Y}$? In other words, the residuals could be split in different bins with their own distribution as the value of $\hat{Y}$ increases, then specific confidence intervals could be applied to each bin. **UPDATE** I'm going to elaborate a bit more because I'm surprised there isn't a simple answer to this question (or an obvious flaw). Let's take the residuals analysis from this thread: ![Residuals vs. Predicted](http://i.stack.imgur.com/s2mAj.png) Here we have a loss of predictive power at the sides of the chart, that is for $\hat{Y} ~ (0.694, 2.23], \hat{Y} ~ (3.51, 9.53]$. Clearly the confidence intervals are not uniform across all values of $\hat{Y}$. Wouldn't make sense to adjust the $\beta$ confidence intervals given $\hat{Y}$?


title: How to test which of five proportions based on Likert scales is highest?

text: I have used a Likert scale in my survey, where respondents answer statements in terms of 1=strongly agree to 5=strongly disagree. There are 5 items measuring a single type of 'need gratification' and there are 5 different needs. After summing the five items for each need gratification, a score is obtained between 5 to 25. These scores have been categorised into four "quartiles": 5-10 is Q1 (first quartile) , 10-15 is Q2 , 15-20 is Q3 and 20-25 is Q4 . **How can I test which of the five needs scales has the largest proportion of responses in the 20 to 25 range (Q4)?**


title: Estimated normal density is fat-tailed?

text: I wanted to check what is the relative error of an estimated N(0, 1) density function. By relative error I mean: $$\alpha(x) = E \bigg ( \frac{f(x|\hat{\mu}, \hat{\sigma})}{f(x|0,1)} \bigg ) - 1$$ where $f(\cdot, \mu, \sigma)$ is the normal pdf. To get an estimate of $\alpha(x)$ I simulated $10^3$ data points from a N(0, 1), estimated sample mean and variance and evaluate the relative error at 100 points $x \in (-6, 6)$. I repeated these steps $reps$ times to calculate the average relative error. This is the R code:               nX = 100     reps <- 1000     X <- seq(-6, 6, length.out = nX)     myM <- matrix(NA, reps, nX)          for(ii in 1:reps)     {      data <- rnorm(1000)      myM[ii, ] <- dnorm(X, mean = mean(data), sd = sd(data))/dnorm(X) - 1     }          plot(X, colMeans(myM), type = 'l', ylab = expression(hat(alpha)),          xlab = "x")      I expected the relative error to be close to 0 for every $x$, but to my surprise it looks like the estimated normal has fatter tails that the real thing. So my questions are: 1) Is this tail behaviour happening because of Jensen inequality (since $f()$ is non-linear)? 2) Suppose I want to use another density estimator, and I want to check how well it does in terms of relative error. I thought that for a normal I could just calculate mean and variance and use the resulting density as a "gold standard" against which I could compare other density estimators. Is there a sensible way to compare how well different density estimators are doing for different values of $x$? Thanks!


title: IRT in R: Does anyone know of an IRT item calibration function that can cope with NA's?

text: I have a dataset that looks like this.               test.takers item1   item2 item3 item4 item5 item6 item7 item8 item9 item10 item11 item12 item13 item14 item16 item17 item18 total_score     tt1             1       1     0     1     0     1     1     1     0     NA      1      1     NA      1      1      1      1          12     tt2             0       1     0     0     0     0     0    NA     1      1     NA      0      1      0     NA      1      0           6     tt3             1       1     1     1     0     0    NA     1    NA     NA      1      1      1     NA      0     NA      0           8     tt4             1       1     1     0     1    NA     1     1     0      0     NA      1     NA      1      0      0     NA           8     tt5             0       1     1     0     1     1    NA    NA     0      1      1      1     NA      0      0      1      1           9     tt6             0       0     0     1     1     1     1    NA     1      1      1     NA      1      1     NA      0      0           9     tt7             1       0     0     1     1     1     1     1    NA      0      1      1      1      0      1      1     NA          11      The dataset consists of 3000 test takers with their responses on an ability test. Not all test takers could respond on all 18 items. So some test takers saw only 10 items, others only 12, etc. The result is a data frame with a lot of NA's. Now I want to calibrate the item parameters using a 2pl irt model. And after that I want to calibrate test taker abilities on the same dataset. The problem I face is that I can't find an R package that can handle calibration on a data frame with missing data. Does anyone know of an R package or R function that could do the job for me or knows how to deal with this problem in some other genious way?


title: Error trying to reduce my data dimentions

text: I'm trying to produce a linear regression model, but I only have 25 observations and 34 predictors. I'm trying feature selection,               library(MASS)     full.m <- lm(fmla, data=mydata)     fsel.m <- step(full.m, direction = "both")      but I get this error,               Error in step(full.m, direction = "both") :        AIC is -infinity for this model, so 'step' cannot proceed      Well, then I tried PCA and factor analysis but I get this weird image from PCA.               r <- prcomp(formula=pca.fmla, mydata, scale=FALSE)      ![enter image description here](http://i.stack.imgur.com/gg45b.png) In factor analysis,               f <- factanal(x=pca.fmla, mydata, factors=2)      I get,               Error in solve.default(cv) :        system is computationally singular: reciprocal condition number = 1.09137e-18      Any help please? EDIT: I can answer this.


title: How to rank the elements for use in Spearman rank correlation?

text: If I have 50 items and I want to compare the top 30. I have several cases to do the ranking , so which one is correct? For simplicity I will replace in the example below 50 items with 5 ( A,B,C,D,E) , and top 30 = 3. So the cases: 1- Rank the 5 items using 2 different functions independently and then compare the top 3 items and remove the items that are not within the 3. Example               Item - Function 1 Rank - Function 2 Rank      A - 1 - 2     B - 2 - 5     C - 3 - 1     D - 4 - 4     E - 5- 3      Now to compare the top 3 we have                  A - 1 - 2         B - 2 - 5         C - 3 - 1      And since 5 is not within 3 range then just remove it and compare the 2 remaining items               A - 1 - 2     C - 3 - 1      2- Rank the 5 items using 2 different functions and the compare the top 3 items from first function with the corresponding items in the 2nd function , even if they are not in the top 3 initially.               Items - Function 1 Rank - Function 2 Rank      A - 1 - 2     B - 2 - 5     C - 3 - 3     D - 4 - 1     E - 5 - 4       Now we select the 3 top items done by function 1. And do the ranking again using function B. So A will change from 2 to 1 since it is the highest, B will change from 5 to 3 since it is the lowest, C will change from 3 to 2               A - 1 - 2     B - 2 - 5     C - 3 - 3      The new ranking will be               A - 1 - 1     B - 2 - 3     C - 3 - 2      Many Thanks !


title: L1 regression versus L2 regression

text: I am trying to develop an intuition for why L1 regression is more expensive than L2 regression. Can somebody point me to some material that explains why this is the case


title: Help with complex model formula in lmer (lme4) for R

text: Most examples about `lmer` formula description in R target rather simple study designs. However, sometimes one is confronted with more complex designs and there is no documentation about a general rule to formula-construction. For instance, how would I formulate the following design:   * 15 subjects (`ID`) treated as between random effect in 3 groups (`Group`) treated as between fixed effect (but would be nice to treat as another between random effect)   * 2 within effects (`Position`, `Size`)   * 1 between random factor primary eye (`Eye`) To elaborate about the design as requested: 15 subjects participated in an HCI user study and were exposed to targets at different Positions of different Sizes. The Time to recognition was measured as the dependent variable. Between Groups we also randomized the visual appearance of targets and delays for display, which was necessary to randomize but only between groups to limit the number of overall factors. However, Group may thus have introduced a random effect. I tried:               lmer(Time ~ Group*Position*Size + (1|ID) + (1|Eye) + (1|Position:ID) + (1|Size:ID))      However, with this formula I'm not sure:   * whether `Group` should also appear in the definition of `1|Position` and `1|Size` as it is actually nested not only within `id` but `id/Group`   * similar question as above for the random factor `Eye`   * how I would go about defining `Group` as random between factor instead of fixed between factor as the group sampling happened randomly


title: Computing non-central moments and normalizer of a quartic exponential distribution

text: Consider a random variable $X$ which has quartic exponential distribution: $$X \sim P(x)=\frac{1}{Z}e^{ax + bx^2 + cx^3 + dx^4}$$ How can one compute $Z$ or non-central moments $E X^k$ given that they exist? As far as I understand, there are no closed-form formulas for these quantities, but is there a good numerical procedure for estimating them? Since I'm very far from being any kind of expert in numerical integration, I'm open to any suggestions that will get the job done with reasonable precision.


title: How to forecast multivariate time-series 'accurately' with a large number of unknown factors using R?

text: I am relatively new to statistics and not formally trained but have been given a complex problem to solve and need some guidance. I realise that I am out of my depth a bit here but would appreciate whatever help I can get bearing in mind that there is no budget for this and as a result it is not possible to purchase software or hire consultants. **The Problem** The business I work for has a large number of mobile representatives that can be dispatched to a variety of different jobs. There are ~100 different job types and each job can be broken up into 4 different final outcomes. Each of these 400 outcomes requires an allocation of man hours to complete. I have a count of how many times each one of these outcomes occurred in each hourband for the past 5 years. I have been asked to forecast how many of each outcome will occur in each hourband for the 28days from the present. The resulting forecast will be used to anticipate staffing requirements on an hour-by-hour basis. As a result the forecasts for each hourband need to been fairly accurate. **Factors** In my data there are clearly some yearly, weekly, and daily seasonal effects. In general each outcome is more likely to occur at certain times of the day on certain days of the week and with some yearly trends. Each different outcome is likely to be related to the frequency of a number of different outcomes. i.e. if _x_ happens then _y_ and/or _z_ are likely but _a_ and/or _b_ are not. There are a large number of environmental factors that contribute to the frequency of each outcome. These can include, but are not limited to weather, sociopolitical, financial trends, one off events. **What I have tried** So far I have tried using simple auto.arima, holtwinters and ets forecasts. holtwinters ended up producing a flat line (i.e. 5 and hour for the next 672 hours). ets doesnt work because the seasons are longer than 24 intervals. auto.arima produced the best results but they were still a long way off being accurate. It was then suggested that I try tbats() and provide it with multiple seasonal lengths. I achieved best results by giving it seasonal lengths of 8760 (1yr) and 168 (1wk). Frustratingly, these results are within 1% when viewed as a sum of all hourbands in a 1 month block but are anything up to 300% (avg 20%) off when considering each individual hourband. Both of these approaches have been applied over an individual outcome rather than considering all possible outcomes (and their correlation to each other). **My thoughts so far** At this stage I feel like my two options are to either to find a way to use something similar to tbats() that will look at the relationships between the multiple different outcomes as well as the seasonality and forecast based on that information. or Abandon that approach for a Neural Network model. My understanding (limited) is that using the Neural Network approach I may be able to 'factor' for the multitude of unknown environmental factors without having to actually identify them. I know this is lazy but my feeling for the data is that there are going to be a fair few unknown factors to identify and forecasting them may end up being a job in itself (i.e. weather conditions) **The Question** (finally) What I am looking for is some guidance. Considering the information above and the fact that I am pretty much limited to R, what is the best approach?? and What are the basic steps I need to follow? While I cant post my data online (due to my employers restrictions) I can send it an individual or two if someone was interested in giving us a hand to find a solution.


title: Neural Network, dependence among outputs?

text: Is there a way to train a neural network in the following manner: You have $n$ observations in the training set. The neural net will start with random weights, and produce $n$ outputs. I want to then apply a function, $f$, on these $n$ outputs. $f$ will output a single real valued number $\alpha$. The goal of the neural net is to then maximize $\alpha$. That is, I don't care what the neural net outputs for any single observation. I want it to maximize some property of the outputs as a whole (e.g. I want to maximize the variance of the outputs). Are there neural networks like this, or other techniques for what I am looking for? Any insight is appreciated!


title: Finding estimator for a one-parameter Weibull distribution

text: I'm doing some practice problems on methods of moments from a textbook. I am stuck on the following question: The pdf of a one-parameter Weibull distribution is given by: $f(x) = \begin{cases} 2\alpha xe^{-\alpha x^2}, &\mbox{if }\: x > 0 \\\ 0 & \mbox{otherwise}. \end{cases} $ Using a random sample of size $n$, obtain a moment estimator for $\alpha$. My Question: I understand that we need to match population moments with sample moments to get the estimator for $\alpha$. In trying to find the first population moment E[X], I set up the following integral: $E[X]=\int^{\infty}_0x2\alpha xe^{-\alpha x^2}$ $E[X]=2 \alpha\int^{\infty}_0x^2e^{-\alpha x^2}$. My question is: how would I best evaluate this integral? Do I need to do integration by parts from here? Or is there an easier "trick" I can use to evaluate this? Thanks in advance.