title: Histogram with uniform vs non-uniform Bins

text: This question describes the basic difference between a uniform and a nonuniform histogram. And this question discusses the rule of thumb for picking the number of bins of a uniform histogram that optimizes (in some sense) the degree to which the histogram represents the distribution from which the data samples were drawn. I can't seem to find the same kind of "optimality" discussion about uniform vs non-uniform histograms. I have a clustered nonparametric distribution with far away outliers so a non-uniform histogram intuitively makes more sense. But I would love to see a more precise analysis of the following two questions:   1. When is a uniform-bin histogram better than a non-uniform bin one?   2. What is a good number of bins for a non-uniform histogram? For a non-uniform histogram, I am considered the simplest case where we take $n$ samples from an unknown distribution, order the resulting $n$ values, and separate them into $k$ bins such that each bin has $\frac{k}{n}$ of these samples (assuming that $n \equiv c k$ for some large integer $c$). The ranges are formed by taking the midpoint between the $\max$ of the values in bin $i$ and the $\min$ of the values in bin $i+1$. Here and here are links that describe these type of non-uniform histograms.


title: Large non-significant effect in a generalized linear mixed model

text: In my experiment I measure accuracy (0-1) in two groups, on two different types of stimuli. I use the generalized linear mixed model because I think this is the only model that can deal with binary DVs and with nested data (trials within the same participant). SPSS syntax:               GENLINMIXED       /DATA_STRUCTURE SUBJECTS=id REPEATED_MEASURES=trial COVARIANCE_TYPE=AR1       /FIELDS TARGET=accuracy TRIALS=NONE OFFSET=NONE       /TARGET_OPTIONS DISTRIBUTION=BINOMIAL LINK=LOGIT       /FIXED  EFFECTS=group stimulusfactorA stimulusfactorB group*stimulusfactorA group*stimulusfactorBstimulusfactorA*stimulusfactorBgroup*stimulusfactorA*stimulusfactorB USE_INTERCEPT=TRUE       /RANDOM USE_INTERCEPT=TRUE COVARIANCE_TYPE=VARIANCE_COMPONENTS       /BUILD_OPTIONS TARGET_CATEGORY_ORDER=ASCENDING INPUTS_CATEGORY_ORDER=ASCENDING MAX_ITERATIONS=100 CONFIDENCE_LEVEL=95 DF_METHOD=RESIDUAL COVB=MODEL       /EMMEANS_OPTIONS SCALE=ORIGINAL PADJUST=LSD       /SAVE PREDICTED_VALUES(PredictedValuere)  PREDICTED_PROBABILITY(PredictedProbabilityre) MAX_CATEGORIES(25) CONFIDENCE(Confidence1re) CONFIDENCE_METHOD(PREDICTEDPROBABILITY) CONFIDENCE_INTERVALS(CIre) PEARSON_RESIDUALS(PearsonResidual1re).      In the fixed effect table with the F-tests, I see a very large effect of Stimulusfactor B (which has 4 categories). However, when I look at the simple contrasts, I see no significant effects at all. The graphs don't show the P-value, but just a red line (instead of a gold one), indicating that it is not significant. I feel something is wrong with this model but I am not sure what it is. Can someone help out? PS When I analyze the averages per condition (percentages correct) in a simple repeated measures ANOVA, the pairwise comparisons do show significant differences.


title: Equivalent to K-S test on discrete data with uneven quantization

text: I have some data resulting from a simulation that consists of several groups, each containing a single real datapoint and a variable number of matched controls. I take the rank of each real value within its distribution of controls, and normalize the ranks to between 0 and 1. Now I would like to test whether, across all groups, the ranks of my real datapoints are different to what I would expect by chance. However, I'm unsure about how to do this, since the distribution of possible normalized ranks across groups is discrete with uneven quantization (due to the variable number of control datapoints in each group). I expect the CDF for the null distribution to look a bit like this: ![Simulated CDF for normalized ranks](http://i.stack.imgur.com/wQj7s.png) I have considered the possibility of bootstrapping a null distribution by drawing randomly from each set of possible normalized ranks, then doing a 2-sample K-S test on the real and null distributions. Is this the way to go, or would there be a more appropriate test? **Edit:** I've also posted a more complete description of the problem I'm trying to solve here.


title: How important is using the exact method for ties in a Cox model, and how long should this take?

text: I'm analyzing the results of some simulation work using a Cox proportional hazard model, and I have what I perceive are a great many ties in the data, representing when a particular individual in the simulation was infected. For example, one of the (many) runs of this model:               Time, Freq     1, 8     2, 9     3, 5     4, 5     5, 6     6, 4     7, 1     8, 6     9, 4     10, 1     11, 3     12, 2      Followed be sporadic ties of two or three individuals all the way up to Time = 55. Worried about ties, I've analyzed the data using `coxph()` using both the Breslow and Efron methods. The Breslow method for this particular data has a log(HR) of 1.95 with a standard error of 4.73. The Efron method produces a log(HR) of 1.75 with the same standard error. Ignoring for a moment the _quality_ of those results in a general sense, they are slightly different, so I'd like to be able to check by using the Exact method. However, this seems to be _extremely_ computationally intensive, and a small data set of ~200 individuals has taken over 24 hours to run. How important is it to use the exact method for handling ties in this analysis? Is there another method someone would suggest? Do parametric models suffer from the same problems with tied data?


title: Representation of ARMA processes

text: **Question** Consider the following process: $$2y_t-3y_{t-1}+y_{t-2}=\epsilon_t-\theta\epsilon_{t-1}$$ What is the model for the process $w_t=\Delta y_t = y_t-y_{t-1}$? **Attempt** I have solved the question in two ways, which seem to be giving conflicting answers.   1. Substitute the form of $y_t$ into the expression of $w_t$: $$w_t=\frac{3}{2}y_{t-1}-\frac{1}{2}y_{t-2}+\frac{1}{2}\epsilon_t-\frac{\theta}{2}\epsilon_{t-1}-y_{t-1} \\ \implies w_t=\frac{1}{2}y_{t-1}-\frac{1}{2}y_{t-2}+\frac{1}{2}\epsilon_t-\frac{\theta}{2}\epsilon_{t-1} \\ \implies w_t = \frac{1}{2}w_{t-1} +\frac{1}{2}\epsilon_t-\frac{\theta}{2}\epsilon_{t-1} $$ This is an ARMA(1,1) model.   2. Use the expanded forms of $y_t$ and $y_{t-1}$, and then substitute into the expression for $w_t$: $$y_t=\frac{3}{2}y_{t-1}-\frac{1}{2}y_{t-2}+\frac{1}{2}\epsilon_t-\frac{\theta}{2}\epsilon_{t-1} \\ y_{t-1} =\frac{3}{2}y_{t-2}-\frac{1}{2}y_{t-3}+\frac{1}{2}\epsilon_{t-1}-\frac{\theta}{2}\epsilon_{t-2} \\ \implies w_t = \frac{3}{2}w_{t-1} - \frac{1}{2}w_{t-2}+\frac{1}{2}\epsilon_t-\frac{\theta+1}{2}\epsilon_{t-1}+\frac{\theta}{2}\epsilon_{t-2}$$ where the facts that $w_{t-1}=\Delta y_{t-1}= y_{t-1} - y_{t-2}$ and $w_{t-2}=\Delta y_{t-2}= y_{t-2} - y_{t-3}$ were used. This is an ARMA(2,2) model. As can be seen, both these approaches yield different models. The first is stationary, while the second is not (by checking the roots of the AR characteristic polynomials). Is this possible, or is there a mistake in one approach?


title: Iterative solving of ML estimators

text: I have derived this likelihood function \begin{equation} \begin{split} &-\frac{1}{N}\log L(\eta,\beta,\mathit{\Omega})\\ &=\frac{1}{2}\log (\sigma^2_{\varepsilon}+T\sigma^2_{c})+\frac{T-1}{2}\log\sigma^2_{\varepsilon}+\frac{1}{2}\frac{1}{\sigma^2_{\varepsilon}+T\sigma^2_{c}}\frac{T}{N}\sum_{i=1}^{N}\left(\bar{Y_i}-\bar{\eta}-\bar{X}\beta\right)^2\\ &+\frac{1}{2}\frac{1}{\sigma^2_{\varepsilon}}\frac{1}{N} \left(Y_i-\boldsymbol{1}\bar{Y}_i-\left(\eta-\boldsymbol{1}\bar{\eta}\right)-\left(X_i-\boldsymbol{1}\bar{X}_i\right)\beta\right)'\left(Y_i-\boldsymbol{1}\bar{Y}_i-\left(\eta-\boldsymbol{1}\bar{\eta}\right)-\left(X_i-\boldsymbol{1}\bar{X}_i\right)\beta\right). \end{split} \end{equation} where $Y_i=(Y_{i1},...,Y_{iT})'$, the time effect $\eta=(\eta_1,...,\eta_T)'$, and define the $T\times k$ matrix $X_{i}=(X_{i1},...,X_{iT})'$, and the $T-$vector with ones $\boldsymbol{1}=(1,...,1)'$. Time averages are here defined as $\bar{Y}_i=T^{-1}\boldsymbol{1'}Y_i$ and $\bar{X}_i=T^{-1}\boldsymbol{1'}X_i$, so that $\boldsymbol{1}\bar{Y}_{i}=\boldsymbol{1}(\boldsymbol{1}'\boldsymbol{1})^{-1}\boldsymbol{1}'Y_i$ and $\boldsymbol{1}\bar{X}_{i}=\boldsymbol{1}(\boldsymbol{1}'\boldsymbol{1})^{-1}\boldsymbol{1}'X_i$. It is now convenient to define $\psi=\sigma^2_{\varepsilon}+T\sigma^2_{c}$. The partial derivatives are \begin{align} \frac{\partial L(\eta,\beta,\mathit{\Omega})}{\partial \beta}=& \sum_{i=1}^{N}(Y_i-\eta- X_i\beta)'\mathit{\Omega}^{-1}X_i\text{ set}=0\\ \frac{\partial L(\eta,\beta,\mathit{\Omega})}{\partial \psi}=& \frac{1}{2\psi}-\frac{T}{2N\psi^2}\sum_{i=1}^{N}(\bar{Y}_i-\bar{X}_i\beta)^2\text{ set}=0\\ \frac{\partial L(\eta,\beta,\mathit{\Omega})}{\partial\sigma^2_{\varepsilon}}=& \frac{T-1}{2\sigma^2_{\varepsilon}}-\frac{1}{2N\sigma^4_{\varepsilon}}\sum_{i=1}^{N}\Big[\left(Y_i-\boldsymbol{1}\bar{Y}_i-(\eta-\boldsymbol{1}\bar{\eta})-(X_i-\boldsymbol{1}\bar{X}_i)\beta\right)' \nonumber\\ & \left(Y_i-\boldsymbol{1}\bar{Y}_i-(\eta-\boldsymbol{1}\bar{\eta})-(X_i-\boldsymbol{1}\bar{X}_i)\beta\right)\Big]\text{ set}=0\\ \frac{\partial L(\eta,\beta,\mathit{\Omega})}{\partial\eta}=& \sum_{i=1}^{N}(Y_i-\eta-X_i\beta)'\mathit{\Omega}^{-1}\text{ set}=0 \end{align} Solving for the parameters yields the following equations \begin{align} \hat{\beta}(\mathit{\Omega})&=\left(\sum_{i=1}^{N}Y_i\mathit{\Omega}^{-1}X_i-\eta'\mathit{\Omega}^{-1}X_i \right)\left(\sum_{i=1}^{N}X_i'\mathit{\Omega}^{-1}X_i \right)^{-1} \\ \hat{\psi}(\beta)&=\frac{T}{N}\sum_{i=1}^{N}(\bar{Y}_i-\bar{\eta}-\bar{X}_i\beta)^2\\ \hat{\sigma}^2_{\varepsilon}(\beta)&=\frac{1}{N(T-1)}\sum_{i=1}^{N}\Big[\left(Y_i-\boldsymbol{1}\bar{Y}_i-(\eta-\boldsymbol{1}\bar{\eta})-(X_i-\boldsymbol{1}\bar{X}_i)\beta\right)' \nonumber\\ &\left(Y_i-\boldsymbol{1}\bar{Y}_i-(\eta-\boldsymbol{1}\bar{\eta})-(X_i-\boldsymbol{1}\bar{X}_i)\beta\right)\Big]\\ \hat{\eta}(\beta)&=\sum_{i=1}^{N}Y_i-X_i\beta \end{align} I can see that the equations have to be solved iteratively, but I'm not sure on the order. Any ideas on how to proceed from here?


title: Curve fitting multivariate data for maximal correlation with univariate data?

text: I have multivariate time series data of the EURUSD financial vehicle. In this data each variable represents a different metric. There are ~200,000 rows and ~20 variables. There are no NULL values for any variable at any row. All data is numerical. Alongside this data, at each time point I have the univariate data "Profit." I want to curve-fit a function to transforms my multivariate data set into a new univariate data set which having the MAXIMAL correlation to my "Profit" variable. In other words, I want to iterate through different mathematical transformations of my multivariate data set until I find the one that is optimally correlated with my "Profit" data. What is the best way to do this? From what I understand, a genetic algorithm should work well.


title: Neural network weights

text: This is going to be a long question : I have written a code in MATLAB for updating the weights of MLP with one hidden layer . Here is the code : weights_1 : weight matrix for input to hidden layer weights_2 : weight matrix for hidden to output layer               function [ weights_1,weights_2 ] = changeWeights( X,y,weights_1,weights_2,alpha )     %CHANGEWEIGHTS updates the weight matrices      %   This function changes the weight of the weight matrix     %   for a given value of alpha using the back propogation algortihm          m = size(X,1) ;     % number of samples in the training set          for i = 1:m         % Performing the feed-forward step          X_i  = [1 X(i,1:end)]    ;            z2_i = X_i*weights_1'    ;              a2_i = sigmoid(z2_i)     ;              a2_i = [1 a2_i]          ;              z3_i = a2_i*weights_2'   ;              h_i  = sigmoid(z3_i)     ;               % Calculating the delta_output_layer          delta_output_layer = ( y(i)' - h_i' )...             .*sigmoidGradient(z3_i')  ; % 3-by-1 matrix          % Calculating the delta_hidden_layer          delta_hidden_layer =  (weights_2'*delta_output_layer)...             .*sigmoidGradient([1;z2_i']) ; % 5-by-1 matrix          delta_hidden_layer = delta_hidden_layer(2:end) ;           % Updating the weight matrices         weights_2 = weights_2 + alpha*delta_output_layer*a2_i ;          weights_1 = weights_1 + alpha*delta_hidden_layer*X_i  ;     end          end      Now I wanted to test it on the `fisheriris` dataset given in MATLAB which can be accesed by `load fisheriris` command . I renamed `meas` to `X` and changed `species` to a `150-by-3` matrix where each row depicts the name of species (as for example first row is `[1 0 0]`) I compute error of the output layer using the following function :               function [ g ] = costFunction( X,y,weights_1,weights_2 )       %COST calculates the error        %   This function calculates the error in the        %   output of the neural network             % Performing the feed-forward propogation       m = size(X,1) ;        X_temp  = [ones([m 1]) X]   ;  % 150-by-5 matrix        z2 = X_temp*weights_1'       ; % 150-by-5-by-5-by-4        a2 = sigmoid(z2)       ;        a2 = [ones([m 1]) a2]            ; % 150-by-5       z3 = a2*weights_2'     ; % 150-by-3       h  = sigmoid(z3)       ; % 150-by-3            g = 0.5*sum(sum((y-h).^2)) ;        g = g/m ;      end      Now in the course the prof gave an example of toy network with 3 iterations , I tested this on that network and it gives the right values but when I test it on the fisheriris data the cost keeps on increasing . And I am not able to understand where it is going wrong . Here is the toy network for which it runs fine : ![enter image description here](http://i.stack.imgur.com/5V1ZQ.png) there is only training example for this set . PS : Ignore the comments ( they are size of matrix used for checking the validity of matrix multiplication for a sample case ) Finally here is the test-bench execute.m , sigmoid.m and sigmoidGradient.m which I have shared just in case to run the functions and test them


title: Low pass filter to maintain edge information

text: I am looking for a kernel as low pass filter that satisfy as:I must find a kernel that statisfies as follows: ![enter image description here](http://i.stack.imgur.com/gVBYW.png) In the my reference paper, the author suggest gaussian kernel that is ![enter image description here](http://i.stack.imgur.com/uWltv.png) The gaussian kernel is very good to approximate that properties. But the kernel reduce edege information when kernel size is big. Could you suggest to me any kernel that can statisfy three above conditions and more robust than gaussian about noise reducing (maintain edge information). I found a kernel as modify kernel. But it very difficult to implement. Thank you so much.


title: Could it be shown statistically that cars are used as murder weapons?

text: I heard a story recently in which someone said if they wanted to kill somebody (and get away with it) they would do it with their car. They cited various statistics about the number of auto-related deaths (including car-on- pedestrians) coupled with additional stats about the number of drivers actually sentenced to any sort of crime...blah, blah, blah. My question is this: Is it statistically feasible to demonstrate that cars ARE (statistically speaking) actually used as weapons to commit murder? In other words, I realize it may not be possible to demonstrate that any single car-on-pedestrian 'accident' was actually an attempted/committed murder. Rather, I'm wondering if a method might be imagined in which it could be demonstrated that some percentage of those 'accidents' are actually, in all likelihood, not accidents at all...


title: General Gage R&R Question

text: I have been doing some reading about gage R&R and I have a simple application for it. I will be using ANOVA for performing my gage R&R. From my standpoint, most people seem to be using gage R&R to determine the relative contribution that various sources make to the total variation. However, when I read about ANOVA I usually think in terms of hypothesis testing and answering a question. In gage R&R, what is the question my ANOVA will be answering?


title: To use or not to use a variable that contains information w.r.t. uncertainty

text: The data I'm looking at is concerned with percentages of people who would recommend a hospital.               E.g.                     NOT         PROBABLY    DEFINITELY          FREQ     hospital_1     10.0        80.0        10.0                100     hospital_2     20.0        50.0        30.0                200      We'd like to do the analysis on the hospital level, so basically I converted the percentages in some kind of weighted score (for example the arbitrarily chosen NOT * 0 + PROBABLY * 1 + DEFINITELY * 2, does not matter), which will yield some kind of single variable called "average_recommendation_score". I'm expecting that this score is unbiased regardless of the number of people who filled in the survey (FREQ), well... in the limit. However, in doing this kind of analysis, where this might actually be a dependent variable, should we also use the FREQ variable for the estimation of parameters (it feels like this given uncertainty about the average_recommendation_score contains extra information; you can imagine how a maximum score is probably not the "real" score if only based on 1 participant's vote)? _Is it OK to disregard the FREQ variable or should it / could it be incorporated somehow yielding a single variable? If it should be, how?_ Note: The goal is inference using some variables that are assumed to be predictors.


title: Which models require stationarity?

text: Regression models (at least up to GLMs) do not traditionally require stationarity (although the requirement for the residuals is even stronger than stationarity). ARMA-style time series models seem to always require stationarity. **What other models always require stationarity?** Vague area: Using a linear regression `y ~ x` when both `x` and `y` are time series. In special cases, the residuals could be i.i.d. etc even though `x` and `y` are highly nonstationary. Of course, when both variables have a clear time trend, time becomes a confounding variable, a problem that matters only if you are interested in interpreting coefficients instead of simply making predictions. Motivating example: Distributed Lag Nonlinear Models (see the `dlnm` package in `R`) use GLMs to estimate the effect of one time series on another time series . I'm having a hard time deciding whether it ever makes sense to use this tool on nonstationary series.


My Dissertation analysis method involves a Multivariate multiple moderated regression. Where I have 4 Dependent Variables (DVs) and 4 Independent Variables (IVs). My model 1 is testing main effects and model 2 is testing interaction effects. Here is what my code looked like:               mod1<-lm(cbind(PerceivedGuilt_Mean, FalseAccusation_Mean, Severity_Mean,      JobOutcome_Mean) ~ Z.FWM_MEAN+ Z.MEAN_RWA+ Gender+Evidence, data=      Unwanted_Sexual_Attention.SAV)          summary(mod1)          mod2<-lm(cbind(PerceivedGuilt_Mean, FalseAccusation_Mean, Severity_Mean,      JobOutcome_Mean) ~ RWA_Evidence_Interaction+ RWA_Gender_Interaction+      FWM_Evidence_Interaction+ FWM_Gender_Interaction,      data=Unwanted_Sexual_Attention.SAV)      My question is for each DV R is giving an omnibus test result. For one of them one interaction term is showing significant. However, when I am running the simple slope for the same analysis, the interaction term is no longer significant. Am I reading the output incorrectly? What is going on? Any help is appreciated.


title: GMM adaptation to new data

text: I have been using the GMM cluster package by Bouman, for which I did not find any adaptation module online. Before I start off reading up on the GMM adaptation theory and implementing it, I did like to know if there are any other opensource GMM projects online which does all of training, testing and adaptation to new data.?


title: Monte Carlo Simulation...?

text: I just have 2 questions: (1) If we can obtain samples from the posterior distribution, is there any need to try to compute posterior expectations and intervals analytically...? (2) Also, I know that Monte Carlo methods will, in general, produce different results every time they are run, but how might I convince someone that my results obtained using MC are reliable? Any help would be appreciated.


title: Bootstrapping in Bayesian Statistics

text: I am reading the book Data Analysis: A Bayesian Tutorial. The book mentions it would be wrong to bootstrap a single data element and use the posterior pdf as the prior. I understand it would be wrong to do this but can you show me either mathematically or intuitively why this bootstrap idea is folly.


title: Cramer's theorem for a precise normal asymptotic distribution

text: I am working on a homework problem for my probability class: (Cramer Application) A. Let $X_1, X_2, ... X_n$ be a sample from a distribution with pdf $f(x;p) = q^xp$. Determine the MLE of $p$ Solution: I was able to do this straight from the definition of of the distribution and the log-likelihood. I get that $\hat{p} = (\bar{x} + 1)^{-1}$, no problem. B. Use Cramer's theorem to obtain the precise normal asymptotic distribution of $\sqrt{n}(\hat{p} - p)$ This is where I am struggling. Is this the delta method? If so, is my function $(\bar{x} + 1)^{-1}$? $\bar{x}^{-1}$? Or just classic, taking $x^2$? And if it isn't the delta method what is the theorem being referenced? Has anyone ever been exposed to a problem like this? EDIT: So it seems to me that this must be the Cramer-Rao variance identity, where $Var(p) = \frac{1}{I(p)}$


title: Why doesn't correlation of residuals matter when testing for normality?

text: When $Y = AX + \varepsilon$ (i.e., $Y$ comes from linear regression model), $$\varepsilon \sim \mathcal{N}(0, \sigma^2 I) \hspace{1em} \Rightarrow \hspace{1em} \hat{e} = (I - H) Y \sim \mathcal{N}(0, (I - H) \sigma^2_{})$$ and in that case residuals $\hat{e}_1, \ldots, \hat{e}_n$ are correlated and not independent. But when we do regression diagnostics and want to test the assumption $\varepsilon \sim \mathcal{N}(0, \sigma^2 I)$, every textbook suggests to use Q–Q plots and statistical tests on residuals $\hat{e}$ that were designed to test whether $\hat{e} \sim \mathcal{N}(0, \sigma^2 I)$ for some $\sigma^2 \in \mathbb{R}$. How come it doesn't matter for these tests that residuals are correlated, and not independent? It is often suggested to use standardised residuals: $$\hat{e}_i' = \frac{\hat{e}_i}{\sqrt{1 - h_{ii}}},$$ but that only makes them homoscedastic, not independent. **To rephrase the question:** Residuals from OLS regression are correlated. I understand that in practice, these correlations are so small (most of the time? always?), they can be ignored when testing whether residuals came from normal distribution. My question is, why?


title: How to have a captured response table for a 2-year prediction base Cox model in R?

text: My study is about "Churn Risk Evaluation of College Students". I decided to use R software, already had a result for the survival probability of an old data (all the students already left or graduated), and also have identified the significant variable using Cox model analysis.   1. I want to have an output of a table called "Captured Response"   2. I want to have a 2-year-ahead prediction using the new data (presently enrolled students who haven't churned yet    3. I want to identify when and what specific students will churn. Is it possible in R, and what will be the codes that I will use? The available covariates I used are `age`, `gender` (female =`0`, male =`1`), `permanent residence` (also numerically coded), `high school` (private =`0`, public =`1`), etc.). Thank you.


title: Does the GB2 distribution have a conjugate prior?

text: I did a bunch of searches on this but did not find anything conclusive. Does the generalized beta distribution of the second kind, also known as the generalized beta prime distribution, have a closed-form conjugate prior? If so, could someone tell me what it is, with a literature citation if you happen to have one handy? If not, I’d like to know that too Any assistance would be greatly appreciated.


title: Is it correct to say the gamma distribution is completely determined by its first two moments?

text: I have heard this said by one of my professors, but he is not available for me to ask anymore. But if this is the case then why are there others defined on the wikipedia page, and why do authors refer to higher order moments in papers?? Many thanks.


title: How to apply Maximum Entropy principle to classify a continuous data set?

text: Maximum entropy classifier is mostly used for Natural Language Processing, in which the datas are discrete. I learned the formalism from the paper _A maximum entropy approach to natural language processing_ and understood most of it. My question is, how to apply the Maximum Entropy principle to do the classification on a continuous data set? I am asking for some sources where I can learn the formalism, not the software packages.


title: Few main variables in the best model still show multicollineairty. Should I remove them despite their importance?

text: I have a binary logistic regression with 5 IVs and all of their 2-way interactions. I have reduced/removed multicollinearity between main variables and their interactions by centering the main variables. I have found one of the best models according to LRT. The model also accords reasonably with theory. Everything is fine. However, one of the odds ratios is very large and its 95% CI is strangely super-broad. I think this is caused by a new case of multicollinearity. This OR is for a main variable which is highly correlated (> 0.8) with a couple of other main variables, according to the correlation matrix. What do you suggest me to do with this new source of multicollinearity which happens to be one of the main variables in this research?   1. Should I remove it from the model and analyze it in isolation?   2. Should I report it as is, but discuss that the very large OR is a case of multicollinearity?   3. Should I report it as is, but add its isolated analysis as well? Do you have any other suggestions that is routinely done in the case of a resisting multicollinearity? Many thanks in advance.


title: Predict binary occupancy vector from history of vectors

text: I have a set of binary vectors where each vector represents one day of occupancy in a house and consists of 48 elements (each element for 30 minutes of the day). Each element can be 1 meaning that house was occupied and 0 for non occupied house. My task is to predict the next day based on the history of the same days (Monday from history of Mondays etc.). So far I am using hamming distance to find 5 most similar days in the history and from them I calculate the probabilities of the occupancy as a mean of those 5 numbers. When the probability is higher than some X, in my case 0.4, I predict it to be occupied. But there is definitely some more efficient way to do this, any algorithms that would capture the trend in the history?


title: Calculate a monthly uncertainty from an annual uncertainty

text: Given a normal distribution of annual mean values (sunny days per year), I know the P50 and the annual uncertainty. Is there any method to calculate the monthly uncertainty? Thank you very much.


title: Reconstructing a vector after projection

text: Suppose one has a matrix of data $X$, which is $n$ observations by $p$ dimensions. Let $P_\perp$ be a projection onto some $k<p$ dimensional subspace. Suppose one computes the principal direction of $P_\perp(X)$. Is there a principled way to reconstruct this principal direction in the full $p$-space?


title: In a random selection of addresses, determine if users are unable to access one of the addresses

text: Our client software randomly selects a server to connect to. There are four separate addresses the software randomly chooses from. Two addresses are "old" and two are "new". As many firms need to explicitly allow outbound access to IP addresses, there will be some firms (users of some firms) that cannot connect to the "new" addresses. We record connections over time. Based on when the new addresses were introduced (June 2012), I'm interested in finding a means of indicating which firms (users of those firms), since June 2012 to now, are unable to connect to the new addresses. (We have no way of explicitly determining this. I am wondering if we can infer it with a probability model of some sort.) So far I have done the following:   * Retrieved connection history of all individual users since June 2012   * Done a count of the number of times each user has connected to one of the "old" addresses   * Discarded any users from the list who have connected to a new address at least once (since this would suggest there is no problem connecting to new addresses) In the end, I have a list of users who have only connected to the old addresses, and how many times they have done so.


title: Choosing a better data-set

text: I have two data-sets for same samples. But they are produced using two different instruments. I want to choose one data-set for further analysis. How can I find/prove which data-set is better?   To make it more simple - let's say we took photos of 50 males and 50 females using two different types of camera, and produce two different data-sets. Both data-sets will have same dimension but the data-values will be different. We want to use one of these data-sets for machine-learning (classification). But the first step will be to choose the better data-set. How to do that? Note 1- As some people suggested simplest way will be to apply machine- learning method on both data-sets and do cross-validation to see which data- set gives better accuracy. This is exactly what I thought in the beginning. Unfortunately this is not a feasible path in case of my data-sets primarily because of the fact that I have several (>1000) data-sets and every data-set have very high dimension with many classes. So I am looking for some kind of statistics which can give some indication about quality of data which will allow to pick few data-sets for further analysis. Note 2 - In this case quality of data means 'In a data-set how well separated are classes compare to other data-set'.


title: testing for mediation of interactions

text: I am somewhat familiar with various ways of testing mediation for factors in different types of regression analysis. (I'm using R and currently working with a multilevel binary logistic regression.) But now I have a situation in which I'd like to test whether one interaction between factors mediates another, and I'm not sure how this could be done properly. To give a simplified example of what I am interested in doing: I have a multilevel binary model using student characteristics to predict pass/fail for students who are in a control versus experimental group. Let's say that the "intervention" appears to affect women more strongly than men, because the interaction gender*intervention is significant. But then adding in a number of co-variates (and their interactions with the intervention), results in a decrease in the magnitude of the coefficient and the significance of the fenale*intervention interaction, suggesting that once we control for these co-variates and their interactions with the intervention, differences between how the intervention "affects" men and women are no longer significantly different. I would like to be able to say something about mediation, and I understand how to test the individual factors for mediation of the gender*intervention interaction, but what if there is another interaction, such as (hours spent on childcare)*intervention, which I think may mediate the gender*intervention interaction? Is there a way to test whether the first interaction mediates the second one? EDIT: As requested, here is a simple example equation which I think explains what I want to do. For the purposes of simplicity, I am specifying this as a simple binary logistic regression model instead of a multilevel binary logistic regression model. Let's suppose there are three IV being used as factors: intervention = whether student was in control or experimental group gender GPA And the DV is whether the student passed or failed. And let's suppose we consider the following models (I'm just listing the factors here, without coefficients or error terms, for ease of readability): M1: OUTCOME = INTERVENTION + GENDER + INTERVENTION*GENDER M2: OUTCOME = INTERVENTION + GENDER + GPA + INTERVENTION*GENDER + INTERVENTION*GPA Suppose in M1 that INTERVENTION*GENDER was significant, so that women benefited significantly more than men from the intervention. Then supposed that INTERVENTION*GENDER was not significant (and the female*experimental coefficient had a smaller magnitude) in M2, and we suspect that this is because INTERVENTION*GPA mediates INTERVENTION*GENDER. What I would like to know is if there is a way to test whether or not INTERVENTION*GPA mediates INTERVENTION*GENDER for these two models....


In an experiment with 5 technicians from different companies, we investigated effects of new equipments. The interested response variable is categorical. One thing I am considering is that the technician effect is a random effect. But it's hard to include a random effect in the multinomial logistic model. If there is not much difference, I would enter the effect as a fixed effect. So my question: How bad is considering a random effect as a fixed effect? From my old statistics class, I can see variability of response variable will be underestimated because of the true random effect and thus might affect statistical tests in the model. But I am not sure how bad it will be.


title: representative baseline for time series-like data

text: I have data for several individuals that takes the following form over time. Many individuals are flat over the course of a year, meaning their measurements stay roughly the same. Others have a peak or two and are then flat. And others have many peaks. For each individual I must find an average (representative) value that is not influenced by the peaks. For example, the overall average will not be a good measure because if an individual has several peaks, this will increase the average. Overall, I need an estimate of the baseline for each individual that is not influenced by the number of peaks an individual has.


title: Which test should be used when paired data (pre/post treatment) have different distributions?

text: My question is if data are in paired groups (before and after treatment), where the pre-treatment group is normally distributed and the post-treatment group is not normally distributed (i.e. it is skewed), then which test should be applied? Is the Wilcoxon sign rank test appropriate?


title: Non-independence of IVs in a random forest model

text: How is a random forest model affected if some of the variables are not independent?


title: How to write the error term in repeated measures ANOVA in R

text: My question is very closely related to a previous post Error term in repeated measures ANOVA in R However, I would like to get more insight into how to define the error term. Suppose I have a two-way repeated ANOVA, The factor for between group effect is the Treatment (control vs. placebo), while Time is the within group effect measured repeatedly over 4 times (T1~T4). Patients ID are recorded as Subject. Here I borrowed the data from an example from the tutorial in http://gjkerns.github.io/R/2012/01/20/power-sample-size.html so the data looks like this                Time Subject Method      NDI      0min    1     Treat 51.01078      15min   1     Treat 47.12314      48hrs   1     Treat 26.63542      96hrs   1     Treat 20.78196      0min    2     Treat 42.61345      15min   2     Treat 32.77171      To apply ANOVA:               aovComp <- aov(NDI ~ Time*Method + Error(Subject/Time), theData)     summary(aovComp)     Error: Subject               Df Sum Sq Mean Sq F value Pr(>F)     Method     1    113   112.7   0.481  0.491     Residuals 58  13579   234.1                        Error: Subject:Time                         Df Sum Sq Mean Sq F value  Pr(>F)         Time          3  13963    4654 103.789 < 2e-16 ***     Time:Method   3   1221     407   9.074 1.3e-05 ***     Residuals   174   7803      45       I have also tried the other error term:               aovComp1 <- aov(NDI ~ Time*Method + Error(Subject), theData)     summary(aovComp1)          Error: Subject                     Df Sum Sq Mean Sq F value Pr(>F)     Method     1    113   112.7   0.481  0.491     Residuals 58  13579   234.1                         Error: Within                  Df Sum Sq Mean Sq F value  Pr(>F)         Time          3  13963    4654 103.789 < 2e-16 ***     Time:Method   3   1221     407   9.074 1.3e-05 ***     Residuals   174   7803      45      Can someone help me explaining the differences between these two error terms? If the first term is the correct one, what does the results from the second error term mean? Thanks


title: Logistic regression model manipulation

text: I would like to understand what the following code is doing. The person who wrote the code no longer works here and it is almost completely undocumented. I was asked to investigate it by someone who thinks " _it's a bayesian logistic regression model_ "               bglm <- function(Y,X) {         # Y is a vector of binary responses         # X is a design matrix              fit <- glm.fit(X,Y, family = binomial(link = logit))         beta <- coef(fit)         fs <- summary.glm(fit)         M <- t(chol(fs$cov.unscaled))         betastar <- beta + M %*% rnorm(ncol(M))         p <- 1/(1 + exp(-(X %*% betastar)))         return(runif(length(p)) <= p)     }      I can see that it fits a logistic model, takes the transpose of the Cholseky factorisation of the estimated covariance matrix, post-multiplies this by a vector of draws from $N(0,1)$ and is then added to the model estimates. This is then premultiplied by the design matrix, the inverse logit of this is taken, compared with a vector of draws from $U(0,1)$ and the resulting binary vector returned. But what does all this **_mean_** statistically ?


title: Visualizing Likert responses using R or SPSS

text: I have 82 respondents in 2 groups (43 in Group A and 39 in Group B) that completed a survey of 65 Likert questions each ranging from 1 – 5 (strongly agree - strongly disagree). I therefore have a dataframe with 66 columns (1 for each question + 1 indicating group allocation) and 82 rows (1 for each respondent). Using R or SPSS does anyone know a nice way to visualize this data. I need something like this: ![enter image description here](http://i.stack.imgur.com/xjoyI.png)   (from Jason Bryer) But I cannot get the initial section of code to work. Alternatively, I found really good examples of how to visualize Likert data from a previous Cross Validated post: Visualizing Likert Item Response Data but there are no guides or instructions on how to create these centered count graphs or stacked barcharts using R or SPSS. I would be very grateful for a step by step guide if possible using R or SPSS. Appreciate your time, Adam


title: Variance as a function of parameters

text: I am currently building a statistical model. According to data, the variance is non-constant, and it is likely based on a factor. Is there any statistical model with variance as a function of other parameters? What is its name, and can that model be fitted in R?


title: How is $\chi^2$ value converted to p-value?

text: I'm brand new to statistics and am studying the math behind split testing (A/B and multivariate). I've learned how to calculate $\chi^2$ with given test data, and I understand how to translate this into a probability via a table, but I'd like to be able to calculate the probability myself. I've read through a couple of explanations online, but I'm not getting it. Does anyone know of a resource or book that breaks this down?


title: How can I correct for residual autocorrelation in a fixed effect panel model?

text: The residuals have an AR(2) structure. Is it appropriate to add AR terms to a fixed-effects panel model?


title: Penalized methods for categorical data: combining levels in a factor

text: Penalized models can be used to estimate models where the number of parameters is equal to or even greater than the sample size. This situation can arise in log-linear models of large sparse tables of categorical or count data. In these settings, it is often also desirable or helpful to collapse tables by combining levels of a factor where those levels are not distinguishable in terms of how they interact with other factors. Two questions:   1. Is there a way to use penalized models such as LASSO or elastic net to test for the collapsibility of levels within each factor?   2. If the answer to the first question is yes, can, and should, this be set up in such a way that the collapse of levels and the estimation of model coefficients occurs in a single step?


title: Post-hoc after GLM

text: I am running a GLM, using the function glm.nb (pscl package) trying to figure how what could influence a particular trait in several locations and years. The output as follow (with slight modification and removing the things beyond this question)                         Estimate Std. Error z value Pr(>|z|)                   state2        -0.163190   0.807307  -0.202  0.00081                   state3        -0.530588   0.783023  -0.678  0.00302                state4         0.942953   0.737328   1.279  0.00094                   year2         -1.759102   0.733214  -2.399  0.01643                year3         -0.633870   0.662903  -0.956  0.33897          I understand the output as such: compared to year1, year 2 is statistically different, but year3 is not different from year1. but what about year2 and year3? For the state factor, state2, state3 and state4 are all statistically different than state1. so I need a post-hoc test between state2 state3 and state4, am I right? thanks


title: Normalized mutual information using cardinalities

text: I have the cardinalities of sets $N_i, \forall i \in 1,2..n$, and the cardinalities $|N_i \cap T|, |N_i|, |T|, |N_i \cup T|, \forall i$, are known. Here, the set $T$ and sets $N_i$'s are all subsets of a mother :) set M and $|M|$ is also known. Using these quantities, I would like to compute the normalized mutual information $NMI(N_i,T)$, for each set $N_i$. Do I have all the necessary info, to compute the $NMI(.)$ using the cardinalities-and what are the assumptions that go with it? In specific, the probabilistic assumptions required over these sets, and the normalization of the mutual information is a grey area to me.


title: How to model sum of two values in JAGS?

text: Imagine, that we measured two values and we know, that in reality one measurement directly corresponds to the latent variable "s1", and the other measurement is in fact sum of two values: "s1" and unknown "s2". Given the measurements input1 for "s1" and input2 for "s1+s2" what is our belief about the distribution of "s1" and "s2" separately? We assume, that s1 and s2 are from normal distribution with the same variance. We measure input[i] together with its SE $=0.2$, that acts as a weight in this "fixed effects model". I did the following JAGS model:               var        input[2], #Input        tau,        s1,         s2     model{     input[1] ~ dnorm(s1,25)     input[2] ~ dnorm(s1+s2,25)     s1 ~ dnorm(0,tau)     s2 ~ dnorm(0,tau)     tau~dgamma(0.001,0.001)     }      As I understand it, it states something like this:   1. Model the variable s1 as mean of normal distribution so the input1 is a probable outcome from it.    2. Favor s1 and s2 such, that `input[2]` is a probable outcome from normal distribution with mean s1+s2 and SD = 0.2.  Because SD for both normal distributions is the same, model should choose s1 and s2 such that `input[1]` and `input[2]` are (on average) equally plausible random samples from N(s1,0.2) and N(s1+s2,0.2). Condition 2. should make the sampler favoring adjusting s2 over s1, since it is the only place that puts any stochastic constraints on s2 (except from a vague prior), in contrast to s1, which is already constrained to make `input[1]` plausible result from N(s1,0.2) So I would expect the model to yield the `mean(s1)` $\approx$ `input[1]`, and `mean(s2)` $\approx$ `input[2] - mean(s1)`... Now what I get is something completely different. What is the error in my understanding? For example, setting the `input[1]<-0.3` and `input[2]<-0.9` produces `s1` $=0.485 \pm 0,254$ and `s2` $=0,321 \pm 0,174$ (mean $\pm$ SD), where I would **I expect`s1` $\approx 0.3$ and `s2` $\approx 0.6$**! I used 120000 adaptation and burn-in steps. I sampled 4000000 (4 million) samples that were thinned with factor 400. Maybe I should use the `dsum` probability distribution (a unique feature of JAGS)? But then, this model does not run:               var        input[2], #Input        tau,        i1,        i2     model{     input[1] ~ dnorm(i1,25)     input[2] ~ dsum(i1,i2)     i1 ~ dnorm(0,tau)     i2 ~ dnorm(0,tau)     tau~dgamma(0.001,0.001)     }      I've got the following error: > >     RUNTIME ERROR: >     Inconsistent arguments for logDensity in distribution dsum >   And I am really not sure how to fix it, besides the fact that this model doesn't let me specify SE of the `input[2]`. * * * Ultimately I am trying to model Bayesian Network Meta-Analysis (NMA) for Hazard Rate (HR) data (survival analysis, proportional hazard rates model). This type of NMA is different to the one in e.g. Gelman A.,Carlin J.B.,Stern H.S.,Rubin D.B. - Bayesian Data Analysis(2003) in that all data I have is the estimate of the _difference_ of log Hazard Rate between groups (as opposed to estimates of each treatment group _separately_ ). So what data tell me is the measured sum (actually the difference) between latent effects of treatment A vs B, (and B vs G, G vs X, A vs X and so on) together with reported SE of the estimate. The model above is a minimal model that illustrate the principle of the NMA with fixed effects for HR, that can be stated in simpler words: > Knowing bunch of values with standard errors of differences between pair of > latent (unknown) parameters describing different combinations of treatments, > what is the best estimate of the difference between treatment X and Y? Is > the difference greater than zero? * * * *I cross posted this question to the J. K. Kruschke's blog and M. Plummer's JAGS forum in hope of reaching a broader audience. I'll update all the posts with any information I gather on this problem. *


title: The meaning of translation completeness w.r.t a random variable

text: I stumbled upon this term in _McFadden - Analysis of qualitative choice behavior_ (page 111). It is said that > "A random Variable $X$ is _translation complete_ if for a function _h_ of > bounded absolute variation with $h(\pm\infty)=0$, the condition > $E(h(X+a))=0$ for all real $a$ implies $h\equiv 0$ (except possibly on a set > of measure zero)." Further it is said that most common distributions have this porperty; in particular, the Gumble distribution is translation complete. Another qoute says that > "A distribution whose characteristic function is nonzero for real arguments > is translation complete [apply Feller, 1966, An Introduction to Probability > Theory and Its Applications, Vol. 2, Wiley, New York, p. 479]..." if I check there I cannot find any reference. **My question are:** If I google _translation complete_ in connection with _random variable_ I get no feasible results. Is _translation completeness_ connected to sufficient statistics: What are complete sufficient statistics? Though I understand the mathematical meaning of _translation completeness_ I dont get the intention of this definition. Why does _translation completeness_ imply $G(X - logK) ) = G(X)^K$ if $G$ is a cdf and $X$ a random variable?


title: When generating a corpus for latent Dirichlet allocation, does word frequency within documents matter?

text: I have my set of documents and have extracted the unique words from them, including a count of the number of times each word appears in the document. But it would seem from the documentation on the Python library I'm using, that word count within documents is immaterial, and that co-occurrence is the main thing. Am I correct in that understanding?


title: Can I use independent-measures ANOVA instead of repeated-measures when I have only one subject?

text: I have some data where I've measured the response times of the same participants in two sessions and within 36 conditions (each participant took part in all 36 conditions). So the data are repeated measures across both IVs, however I've been informed that I can run it as though it was an independent- measures ANOVA because all the data comes from one participant. Is that true? If so, why does having only one participant make it OK?


title: Continuous response distribution bounded by zero

text: My lab has data from an experiment with a treatment and a control, and pre and post tests. The data will be analyzed by standard 2X2 repeated measures ANOVA. However, the response variable is of cell counts in tissue expressed as a percentage of total cells. There are many subjects with 0%'s and the rest somewhere between 1 and 10%. There are many fractions of a percent, so it is not count (i.e. suitable for Poisson). I'm worried that assuming it is a continuous normal distribution is making it more likely that an increase in %cells is being declared significant (when in reality there is no difference; Type I error) because the model assumes variation above and below the mean, when in reality the variation will primarily be upwards due to the tight bound with zero.   1. Am I thinking about this data in the right way and, if so, what would be a suitable distribution for this type of data?   2. Can this easily be incorporated into PROC GLM in SAS or using a package in R?


title: Estimating effect of binary event on multiple time series

text: I have multiple time series where at any point in the time series, an event can occur that I believe has an effect on the time series. This event can happen at different times for each of the different time series. How would I go about estimating the effect of this event on a time series where the event has not yet occurred? Ideally this would be some function of the time series itself. Thanks!


title: Pitfalls of using k-sample Anderson-Darling test to perform two sample Anderson-Darling test

text: Could I please check if there are any theoretical or practical pitfalls if one uses implementations of the $k$-sample Anderson-Darling test [3] to perform the two sample Anderson-Darling test [1,2]? For example, there is the R package adk, which is an implementation of the $k$-sample Anderson-Darling test. * * *   1. Anderson, T. W. and Darling, D. A. (1952) 'Asymptotic Theory of Certain "Goodness of Fit" Criteria Based on Stochastic Processes'. The Annals of Mathematical Statistics, **23** , 193-212.   2. Anderson, T. W. and Darling, D. A. (1954) 'A Test of Goodness of Fit'. Journal of the American Statistical Association, **49** , 765-769.   3. Scholz, F. W. and Stephens, M. A. (1987) '$K$-Sample Anderson-Darling Tests'. Journal of the American Statistical Association, **82** , 918-924.


title: Can we do a regression for statistically insignificant data?

text: I'm testing the sensitivity of different factors on the benefits of precision irrigation. The change in the inputs is changing the output but not significantly. How is it possible to do a regression for these factors?


title: Can Random Forest Methodology be Applied to Linear Regressions?

text: Random Forests work by creating an ensemble of decision trees where each tree is created using a bootstrap sample of the original training data (sample of both input variables and observations). Can a similar process be applied for linear regression? Create k linear regression models using a random bootstrap sample for each of the k regressions What are the reasons NOT to create a "random regression" like model? Thanks. If there's something I'm just fundamentally misunderstanding then please let me know.


title: How do you interpret the coefficients in PROC GLIMMIX and PROC GENMOD?

text: For clustered data, how would you interpret the coefficients? Would you just exponentiate the coefficient to get an odds ratio?


title: Tool to draw normal QQ-plot with 45 degree reference line

text: I have some data (a set of numbers) and I want to compare them with the normal distribution using QQ-plot. The only statistical tool that I am aware of is Octave but Octave doesn't draw the reference line. ![enter image description here](http://i.stack.imgur.com/RMTrF.png) As can be seen in the above picture the x-axis has a range from - 3 to 3 but the y-axis has a range from 20 to 140. I think drawing a 45 degree slope line would be wrong. How can I draw such a reference line in Octave? Or does there exist any easy tool that would draw a QQ-plot with that reference line?


title: Can somebody offer an example of a unimodal distribution which has a skewness of zero but which is not symmetrical?

text: In May 2010 Wikipedia user Mcorazao added a sentence to the skewness article that "A zero value indicates that the values are relatively evenly distributed on both sides of the mean, typically but not necessarily implying a symmetric distribution." However, the wiki page has no actual examples of distributions which break this rule. Googling "example asymmetrical distributions with zero skewness" also gives no real examples, at least in the first 20 results. Using the definition that the skew is calculated by $ \operatorname{E}\Big[\big(\tfrac{X-\mu}{\sigma}\big)^{\!3}\, \Big]$, and the R formula sum((x-mean(x))^3)/(length(x) * sd(x)^3), I can construct a small, arbitrary distribution to make the skewness low. For example, the distribution x = c(1, 3.122, 5, 4, 1.1) yields a skew of $-5.64947*10^{-5}$. But this is a small sample and moreover the deviation from symmetry is not large. So, is it possible to construct a larger distribution with one peak that is highly asymmetrical but still has a skewness of nearly zero?


title: SEM with binary dependent variable

text: Much like with regression, handling binary dependent variables in SEM requires special considerations. In particular, some of these are noted on Dave Garson's Structural Equation Modeling and include: >   1. Polychoric correlation. LISREL/PRELIS uses polyserial, tetrachoric, and > polychoric correlations to create the input correlation matrix, combined > with ADF estimation (see below), for variables which cannot be assumed to > have a bivariate normal distribution. > >     * Sample size issue. ADF [Asymptotically distribution-free] estimation > in turn requires a very large sample size. Yuan and Bentler (1994) found > satisfactory estimates only with a sample size of at least 2,000 and > preferably 5,000. Violating this requirement may introduce problems greater > than treating ordinal data as interval and using ML estimation. This is also > a reason cited for preferring the Bayesian estimation approach to ordinal > data taken by Amos since Bayesian estimation can handle smaller samples than > ML or ADF. > I'm currently trying to use the package sem in R to test my model, and the author of the model suggests using polychoric correlations on R-help. The problems are:   1. I don't know what estimation method is being used with these correlations (i.e., ADF or ML).    2. My sample size is small (N = 173).    3. I'm not familiar with how to interpret polychoric associations (in the case that it is appropriate for me to use them). All the other variables in my model are continuous in nature.  Any help and/or links would be greatly appreciated. I'm also considering using other software like OpenMX, but I'm still reading about how it handles binary data. Help with what other software I might want to use would also be appreciated.


title: Measuring Document Similarity

text: To cluster (text) documents you need a way of measuring similarity between pairs of documents. Two alternatives are:   1. Compare documents as term vectors using **Cosine Similarity** - and TF/IDF as the weightings for terms.   2. Compare each documents probability distribution using f-divergence e.g. **Kullback-Leibler divergence** Is there any intuitive reason to prefer one method to the other (assuming average document sizes of 100 terms)?


title: How to make stochastic gradient descent algorithm converge to the optimum?

text: (Background info taken from my blog) In logistic regression, the hypothesis function, which models the relationshiop between the dependent variable $P(y = 1)$ and the independent variable $X$, is : \begin{align*} H_i = h(X_i) &= \frac{1}{1 + e^{-X_i \cdot \beta}} \end{align*} where $X_i$ is the $i$th row of the design matrix $\underset{m \times n}{X}$, or in matrix form: \begin{align*} H &= \frac{1}{1 + e^{-X \beta}} \end{align*} H is a $m\times 1$ matrix. Except for $X\beta$ all operations are element- wise. The cost function $J$ is a measure of deviance of the modeled dependent variable from the observed $y$ \begin{align*} J &= (1/m)\sum_{i = 1}^m [-y_i\log H_i - (1-y_i)\log (1-H_i)] \\ &= (1/m)\sum_{i = 1}^m \left[ -y_i \log \frac{1}{1+e^{- X_i \cdot \beta}} - (1 - y_i) \log \left( 1 - \frac{1}{1+e^{- X_i \cdot \beta}} \right) \right] \\ &= (1/m)\sum_{i = 1}^m \left[ y_i \log (1+e^{-X_i \cdot \beta}) + (1 - y_i) \log \left( 1+e^{X_i \cdot \beta} \right) \right] \\ \end{align*} \begin{align*} \frac{\partial J}{\partial \beta_j} &= \dfrac {1} {m} \sum_{i=1}^m \left[ y_{i}H_{i}e^{-X_{i}\cdot \beta }\left( -X_{ij}\right ) + \left( 1-y_{i}\right) \dfrac {1} {1+e^{X_{i}\cdot\beta }}e^{X_{i}\cdot\beta }X_{ij} \right] \\ &= \dfrac {1} {m} \sum_{i=1}^m \left[ y_{i}H_{i}e^{-X_{i}\cdot \beta }\left( -X_{ij}\right ) + \left( 1-y_{i}\right) H_i X_{ij} \right] \\ &= \sum_{i=0}^m H_{i}X_{ij}\left( -y_{i}e^{-X_{i}\cdot \beta }+1-y_{i}\right) \\ &= \sum_{i=0}^m H_{i}X_{ij}\left( 1-y_{i}\left( 1+e^{X_i\cdot\beta }\right) \right) \\ &= \sum_{i=0}^m H_{i}X_{ij}\left( 1-y_{i} / H_i\right) \\ &= \sum_{i=0}^m \left( H_{i}-y_{i}\right) X_{ij} \\ &= (H - y) \cdot X_j \\\\ \frac{\partial J}{\partial \beta} &= X^T (H-y) \end{align*} Let $f_i' = \left( H_{i}-y_{i}\right) X_{ij}$, then according to this video: ![enter image description here](http://i.stack.imgur.com/6VIAl.png) batch gradient descent can be described as: Until convergence: for all $j$: $$\theta_j := \theta_j - \alpha \sum f_i'$$ and stochastic gradient descent can be described as: Shuffle the rows of data, and until convergence: for all $i$ in $1\cdots m$: for all $j$ in $0\cdots n$: \begin{align*} \theta_j := \theta_j - \alpha f_i' \end{align*} This looks straight-forward, but when I implement stochastic gradient descent in R, it's unable to converge anywhere close to the optimum, here is the code:               logreg = function(y, x) {         alpha = 1.15         x = as.matrix(x)         x = cbind(1, x)         m = nrow(x)         m1 = sample(m)         n = ncol(x)              b = matrix(rep(1, n))         newb = b + .1         h = 1 / (1 + exp(-x %*% b))         J = -(t(y) %*% log(h) + t(1-y) %*% log(1 -h))         newJ = J+.5              while(1) {             cat("outer while...\n")             for(i in m1) {                 Vi = exp(-as.numeric(x[i, ]%*%b))                 Hi = 1 / (1 + Vi)                 Ei = (Hi - y[i])                 sDerivJ = matrix(Ei * x[i, ])                 newb = b - alpha * sDerivJ             }             h = 1 / (1 + exp(-x %*% newb))             newJ = -(t(y) %*% log(h) + t(1-y) %*% log(1 -h))             if((newJ - J)/J > .15) {                 alpha = alpha/2                 next             }             print(b)             print(newb)             b = newb             J = newJ             if(max(abs(b - newb)) < 0.001)             {                 break             }         }         b     }          nr = 5000     nc = 20     set.seed(17)     x = matrix(rnorm(nr*nc), nr)     y = matrix(sample(0:1, nr, repl=T), nr)     testglm = function() {         res = summary(glm(y~x, family=binomial))         print(res)     }     testlogreg = function() {         res = logreg(y, x)         print(res)     }     print(system.time(testlogreg()))     print(system.time(testglm()))      I am wondering what went wrong.


title: Discrete-Time Event History (Survival) Model in R

text: I'm trying to fit a discrete-time model in R, but I'm not sure how to do it. I've read that you can organize the dependent variable in different rows, one for each time-observation, and the use the `glm` function with a logit or cloglog link. In this sense, I have three columns: `ID`, `Event` (1 or 0, in each time-obs) and `Time Elapsed` (since the beginning of the observation), plus the other covariates. How do I write the code to fit the model? Which is the dependent variable? I guess I could use `Event` as the dependent variable, and include the `Time Elapsed` in the covariates. But what happens with the `ID`? Do I need it? Thanks.


title: SAS and Principal Components Analysis

text: I'm new to SAS coding from SPSS point and click. I don't currently have SPSS, so SAS is my only option. My question is: how do I get the percent variance explained and cumulative variance explained following varimax rotation?? I get the total that I would get under 'Rotation Sums of Square Loadings' but not the % variance or cumulative %. Example SAS code below:               PROC FACTOR DATA=survey     METHOD=PRIN     PRIORS=ONE     NFACT=2     MINEIGEN=1     SCREE     ROTATE=VARIMAX     FLAG=.32     OUT=ex4comp4;     VAR pn1 pn2 pn3 pn4 pn5 pn6 pn7 pn8 pn9 pn10 pn11 pn12 pn13 pn14 pn15 pn16 pn17 pn18 pn19 pn20;     RUN;      Essentially I want to know what % variance each factor explains but all I get is the Total. So frustrated and confused! Thank you in advance for any suggestions or solutions KWAK


title: How to calculate multi-step prediction intervals for time series data?

text: I need to manually calculate multi-step prediction intervals for time series data. I know packages like 'forecast' in R provide these, but I cannot use these packages as the production infrastructure does not support R. Could someone point me to an article that explains how to generate multi-step prediction intervals?


title: Can I perform a stratified test againts a true value?

text: I have a study where they randomized patients into 2 groups (A and B). The primary aim is to compare group A against a fixed/ true value of 0.75. The binomial test is not stratified and cochran mantel haenszel needs 2 groups. How can I perform a stratified test for risk differences against a true value?


title: How to extend PDF of normalized sample to original sample?

text: To calculate the PDF function using Shannon entropy I have scaled my original sample by simply doing $x'=(x-a)/(b-a)$; where $b=\text{max}(x)$, and $a=\text{min}(x)$ and then I found the $\text{PDF}(x')= exp(-\lambda_0-\lambda_1x'-\lambda_2x'^2)$ using optimization with subjection to moments constraints. However it's customary to plot PDF against the real sample (x). Therefore I wonder how I can map the PDF(x') for real sample (x). Thank you.


title: The effect of oversampling on the positive predictive value

text: I need to calculate the positive predictive value for a validation set for a rare event. The problem is that the validation set was oversampled for the rare event. The event occurs in 5 percent of the population, however the oversampling has adjusted it to be in 50 percent of the sample. How does the oversampling effect the calculation of the ppv?


title: Is it possible to model BOTH censoring and truncation in BUGS?

text: Survival times are often right censored and left truncated. From my experience, it does not seem like OpenBUGS allows for both. Truncation is denoted as T( , ) and censoring as C( ,). For instance, a censored observation known to lie above 10 days, but cannot have been observed below 2 days (left truncation) seem like it should be coded as thus:               y[i]~dexp(lambda)T(2 , )C(10, )      where the absent second values denote no upper limits. But it seems OpenBUGS can only do one or the other; using both in the code invariably results in a compilation error. Since this is a typical and reasonable modeling scenario, is there any way to do this in OpenBUGS?


title: Is eigen-gene the principal components of gene dataset just like eigen-faces in face recognition?

text: I have read about eigen faces and eigen voices etc in signal processing which involves PCA. They used to refer to the principal components selected to process for image recognition etc. But now I heard from somewhere about eigen gene in genetics and computational biology. OK if they are the principal components, how do we assign numbers to nucleotides A, G , T and C etc which are string data. Or isn't the data format like strings?


title: Analysis using ANOVA

text: I want to compare 6 designs of spoons (D1,D2,D3,D4,D5,D6) in 20 children (blocks).I also want to see whether holding the spoon in right or left hand affects food-pinching response(number of M&M's picked and placed in cup). I have dataset with response % noted for each design type on all of the 20 children, and the information about which hand is used is represented as zero or one. I am using Randomized Complete Block Design (with children as blocks) to perform the test. I can use ANOVA in R to check the effect:               summary(aov(response~blocks+designs))      Should children and hand data be included in the model. Would the information about hand impact the outcome. **EDIT** : I am confused about how to compare whether holding spoon in right or left hand affects the number of M&M's picked. Should this information change the above mentioned ANOVA analysis. Below is the sample of data set which includes data for only 10 children for two design types.               Response  Design Children Hand     20.11   D1  1   0     30.26   D1  2   1     28.56   D1  3   1     23.19   D1  4   1     23.29   D1  5   0     27.62   D1  6   0     29.0    D1  7   1     30.5    D1  8   1     28.23   D1  9   0     29.98   D1  10  1     25.7    D2  1   0     26.2    D2  2   1     19.2    D2  3   1     24.32   D2  4   1     25.3    D2  5   0     30.45   D2  6   0     28.98   D2  7   1     18.85   D2  8   1     32.7    D2  9   0     29.8    D2  10  1      Any pointers will be helpful.


title: Estimate the proportion of the variance of $Y$ taken into account by the linear regression $Y=aX+b$ and regression coefficient r

text: It is a homework question from coursera course Core Concepts in Data Analysis. I think we need to know at least variance of $X$ to calculate the value. The staff gives a hint we can calculate it with the coefficient of determination. But I have no idea what the coefficient is.


title: High correlation will help detect spurious regression over cointegration?

text: I'm analyzing two financial time series with Johansen method. A high Correlation coefficient using the Pearson method will help me to detect spurious cointegration models to avoid? If this is not the case, which is the best method would provide clue about it? Thank you


title: Is multilevel modelling simpler, more practical, or more convenient using Bayesian methods or frequentist methods?

text: In this community wiki page a twice-upvoted comment asserted by @probabilityislogic asserted that "Multi-level modelling is definitely easier for bayesian, especially conceptually." Is that true, and if so/not why?


title: Pre Processing data for prediction or classification

text: I have a time series which I can preprocess by using centering and scaling. In this respect I have the following questions   1. Should I preprocess the response as well as the variables or only the predictive variables?    2. If I preprocess the dependant variable as well as the independant variables should not I reverse the preprocessing of the dependant variable for the predictions as to return it to the mean and scale of the actual response? Many thanks


title: ANCOVA with complicated research design

text: I am having a few problems with my analysis and hope that somebody can help me out. About my design: I have 120 respondents, who I showed each 7 facebook posts out of a pool of 16 (randomly assigned). Each facebook posts was "ranked" by them on 1 DV (shares) and 3 IV (entertainment, information and emotion). Due to the repeated measures design I decided to use ANCOVA and include "participant" as a factor, controlling for general differences among participants (so 'groups'), and my three IV's as covariates. I am not interested in group differences though. Up to now makes sense, doesnt it? I could also perform a linear regression but it would be kind of crazy to create a dummy for 120 participants? Now I have the problem that I dont meet the assumptions. If I check for interactions between the factor (participants) and IV's, they are all significant. And in a model with participant as factor and the three covariates, the Levene's test is significant as well...Nevertheless, overall there is a linear relationship between DV and IV's. But each group is only 7 responses, this may be too small? Anybody knows why and can help? Since I use ANCOVA for another purpose than most of the people anyway, I am not sure if I can ignore the violations or how I can help it. Highly appreciate any help!! :) Wiebke


title: Possible Problem in Naive Bayes Spam Filters

text: I realized that, never been observed words in training sets, decrease test probaility of documents being spam and not spam after i mentioned here and this link "... you need to account for it as a highly unprobable one." was helpful. I think it is an important problem. Spammers might use this by adding their spam mails, words probably not in mail corporations dictionary (for example using rarely used languages words, meaningless words etc. ) Isn't this a big problem?


title: Residuals in poisson regression

text: Zuur 2013 Beginners Guide to GLM & GLMM suggests validating a Poisson regression by plotting Pearsons residuals against fitted values. Zuur states we shouldn't see the residuals fanning out as fitted values increase, like attached (hand drawn) plot. But I thought a key characteristic of the Poisson distribution is that variance increases as mean increases. So shouldn't we expect to see increasing variation in the residuals as fitted values increase? ![enter image description here](http://i.stack.imgur.com/YQTRP.jpg)


title: Should I use ANOVA or MANOVA?

text: I have data on two groups...I have one independent variable (times 1, 2, and 3) and a dependent variable fitness score (which includes score of a person's heart rate in five stages of the test, minutes they achieved, distance they achieved and rated perceived exertion). I'm leaning toward one-way repeated measures ANOVA, but I'm also starting to wonder if I should be using one-way repeated measures MANOVA because my dependent variable has 4 levels. Or does this not matter? Could I still use repeated measures ANOVA? Thank you.


title: Understanding Maximum Entropy Bootstrap (meboot) algorithm

text: The Maximum Entropy Bootstrap (meboot) let's you create random realizations of time series. I have tried to make sense of the description of the algorithm in H. D. Vinod's paper **'Maximum entropy ensembles for time series inference in economics'** where it is divided into 7 steps.   1. First the data $x_t$ is sorted to create order statistics $x_(t)$, then intervals $I_t$ around these sorted points are constructed and an according set of uniform densities.    2. Next the mean $m_t$ in the interval $I_t$ is computed as a weighted sum of the order statistics $x_(t)$. I am having problems understanding the 2nd step and consequently the following steps in the paper. What does recomputing the mean of a uniform distribution in a given interval signify? Do I have to change the interval to satisfy the new mean? Thank you in advance to anyway who can provide me with an explanation!


title: Parameters of $ y_i = \beta_0^2 + \beta_0 \beta_1 x_i$

text: I have this model, nonlinear in the parameters $ y = \beta_0^2 + \beta_0 \beta_1 x_i $ exist a known estimation of parameters ?


title: How should I interpret these strange density and mixing plots when fitting a generalised pareto distribution using MCMC with JAGS?

text: I'm trying to fit a generalised pareto distribution to a simulated dataset using JAGS and runjags. When doing so, I get very strange density and mixing plots for the mu parameter. The sigma and xi parameters are estimated as expected. **How should I interpret these plots, and proceed to fix my model?** The plot and a reproducible model are below. ![Diagnostics for mu](http://i.stack.imgur.com/lcwLt.png)               # Simulate gpd     require(evir)     N <- 1000     set.seed(1984)     x <- rgpd(N, xi = 0.5, mu = 5, beta = 2000)          # fit distribution using JAGS     require(runjags)     require(rjags)          library(mcmcplots)          model <- "model {       for (i in 1:N) { #data# N             x[i] ~ dgenpar(sigma, mu, xi) #data# x         }         mu ~ dunif(0.0001, 1e6)       xi ~ dgamma(0.001, 0.00001)         sigma ~ dunif(0.000001, 10000)     #monitor# sigma, mu, xi     }"          results <- run.jags(model, method="rjags", n.chains = 3, modules="runjags", inits = list(list(sigma=1, mu=5, xi=1),                                                                                list(sigma=5, mu=1, xi=10),                                                                                list(sigma=100, mu=0.1, xi=50)))     load.runjagsmodule()          mcmcplot(results, dir=getwd())      Note that increasing the chain length doesn't solve the problem. After increasing the sample from 4,000 to 400,000, the density for mu isn't anything like a unimodal distribution centred around 5: ![enter image description here](http://i.stack.imgur.com/yogpJ.png)


title: Mixed-design with split-plot and mixed effect

text: From Wikipedia > a **mixed-design** analysis of variance model (also known as a **split- > plot** ANOVA) is used to test for differences between two or more > independent groups whilst subjecting participants to repeated measures. > Thus, in a mixed-design ANOVA model, one factor (a fixed effects factor) is > a between-subjects variable and the other (a random effects factor) is a > within-subjects variable. Thus, overall, the model is a type of **mixed > effect model**.   1. what relations are betwen mixed design and split-plot design?  Is a split-plot design a within-subject design, or a between-subject design, or both?   2. Why is the model of a mixed design a mixed effect model? Which input variable has random effect and which input variable has fixed effect? Thanks!


title: Difference between BIRCH and Two-Step Clustering

text: I don't have SPSS. But judging from the documentation you can find on the internet, SPSS "Two-Step Clustering" seems to be closely related to BIRCH clustering: > Zhang, T., Ramakrishnan, R., & Livny, M. (1996, June).   >  **BIRCH: an efficient data clustering method for very large databases.**   >  In ACM SIGMOD Record (Vol. 25, No. 2, pp. 103-114). ACM. In particular, the first step of Two-Step clustering apparently is to build a **CF-Tree** , a structure which is introduced in the BIRCH article. Does anyone have an insight into what the _differences_ of BIRCH and Two-Step clustering are, or where additional details can be found?


title: Inclusion of significant interaction term in logistic regression table versus stratification for data presentation

text: This is a general question on logistic regression result reporting for a publication. We have an example where two well correlated ($r=0.4, p=0.001$) blood parameters (`blood parameter1` and `blood parameter2`) are associated with blood pressure in mixed sample of men and women, however, the effect is much stronger in women for blood `parameter1. UNIANOVA analysis includes `age`, `gender`, and one of the two blood parameters at a time as covariates, and `blood pressure` as independent variable. Each analysis shows the blood parameters as significant, however, when including `gender × blood parameter1`, both the interaction term and `blood parameter1` are significant. In women but not men `blood parameter1` is highly correlated with `blood pressure`. The same interaction term using `blood parameter2` (`gender × blood parameter2`) is not significant. We are wondering whether for the publication of the data, we should include the table showing the interaction term where it is significant and not further include it/show it in the table where it isn't? I.e. show `age`, `gender`, `blood parameter1`, `gender × blood parameter1` in one table and `age`, `gender`, `blood parameter2`, in the other? I would be grateful for an opinion. Alternatively, we could stratify the analysis for `blood parameter1` by gender. The information seems to be the same, but it might be perceived in a different way.


title: Minimum number of instances to create a decision tree

text: I'm kind of new doing data mining, so sorry if my question is not very clear. I'm working in a project that is aiming to do data mining over the interactions of the students with a e-learning platform. So, I'm trying to generate decision trees with some data that I collected (number of times that student use a resource, number of activities done, number quizzes taken, etc). The trees that I'm getting are interesting for me and according with the cross-validation tests, it's kind of accurate. So my doubt is the following, I generated the tree with the data collected over two months and there were 44 students in the course. Is this enough data to trust in the tree? I have only the 44 instances that gather the two months of interactions... Thanks in advanced for your help.


title: Simulating Gaussian Mixture signals

text: I am trying to simulate a signal produced from a GMM of 3 mixtures. Here's MATLAB code               mix1 = -0.4 + 0.01*randn(10000,1);      mix2 = -5 + 0.02*randn(10000,1);     mix3 = 10 + 0.1*randn(10000,1);          gmm = mix1 + mix2 + mix3;          [f,x] = ksdensity(gmm);     plot(x,f)      I expected to see 3 mixtures (i.e. peaks), instead I see just one gaussian. Why is that?


title: What are some examples of public datasets that have randomized instruments?

text: Sometimes they ask questions in different orders, or use different prompts. Or datasets with instruments (with at least one variable randomized)? I would like to use at least one of them for my causal modelling course (Stat 566), whose syllabus is at https://www.stat.washington.edu/tsr/s566/syllabus566.pdf


title: Steps in making a Global Regression Model

text: I saw a journal article [1] saying that he constructed the following S-curve model : $$y=\exp\left({\beta_1 + \frac{\beta_2}{x}}\right) + \mathrm{residual}$$ The topic was about a global regression model. **I just cant understand how did he make it. Do you guys know how?** I emailed him already but I guess he would not have time to answer. By the way this what he said in the paper: > Since flood damage increases with flood depth, the following S-curve model > was constructed: > > $$y=\exp\left({\beta_1 + \frac{\beta_2}{x}}\right) + \mathrm{residual}$$ [1] L.F. Chang & M.D. Su (2007),   Using the Geographically Weighted Regression to Modify the Residential Flood Damage Function   _World Environmental and Water Resources Congress 2007: Restoring Our Natural Habitat_


title: Clustering can be plotted only with more units than variables?

text: I am using R software (R commander) to cluster my data. I have a smaller subset of my data containing 200 rows and about 800 columns. I am getting the following error when trying kmeans cluster and plot on a graph: > 'princomp' can only be used with more units than variables I then created a test doc of 10 row and 10 columns whch plots fine but when I add an extra column I get te error again. Why is this? I need to be able to plot my cluster. When I view my data set after performing kmeans on it I can see the extra results column which shows which clusters they belong to. Is there anything I am doing wrong, can I ger rid of this error and plot my larger sample?


title: Machine learning with trinomial features

text: I have 100,000 students who have each answered some multiple choice questions. Given their performance I want to work out what the chances are of a particular student answering the next question correctly.   * not each student has answered every question   * the probability of getting a question correct by chance changes with each question Can I turn this into a machine learning problem using past questions as features? If so, how should I code them up? From what I know of other ML systems (e.g. bag of words) people generally use 1 to show presence of a word and 0 to denote its absence. But in this case I would need 1 to show they got the question right, -1 if they got it wrong and 0 if they didn't answer it. Or can I just use 1 for correct, 0 for incorrect and then omit the feature if the student never answered that question? Which algorithms are best suited for this kind of problem? I've been reading a bit about Latent Trait/ Latent Class analysis and Item Response Theory but unsure how to implement them. Using Python/Octave/Matlab. The representation of the dataset in Zhubarb's answer below is correct btw. To make my question clearer, imagine i have been given a new question today, to which nobody knows the answer. I can see that 20 people have answered it and I can look through their exam history where i know which questions they got right and which they got wrong. I want to find out who the smartest person is so that I know who to copy but everyone has answered different questions and different numbers of questions (e.g. person 1 may have answered 20 and got 18 correct whereas person 2 has answered 12 and gotten them all correct. Person 3 has 5/5 correct but they were all quite easy questions, whereas person 4 only got 3/4 but the ones they got right were tricky).


title: Granger Causality coefficient comparison

text: I am working with a time series data and found that $a_t$ Granger causes $b_{t+1}$ and $b_t$ granger causes $a_{t+1}$. The results were obtained through Stata.   1. With the coefficients obtained, is it possible to determine which Granger causal direction is stronger? If it can't be determined by just looking at the coefficients, is there a method to do so?  Also, Granger causality was tested on 2 conditions. In one subset, a policy was adopted. In the other subset, the policy was not adopted. The results was that a granger causes $b_{t+1}$ in both cases.   2. Is there a way to test whether the Granger causality coefficient in condition one was significantly different in condition two? One idea is to introduce an interaction term, but is this possible with Granger causality tests?


title: How to stationarize profit and loss data with an increasing variance and large negative values for time series analysis?

text: PnL can take large negative values, and its variance increases over time as the firm grows. If we do differencing, an increasing variance remains. If we take log, negative values cannot be defined. How do you stationarize this type of data for time series analysis? Please assume that we cannot decompose PnL into underlying factors such as profits and losses, which can be experessed as positive and modeled separately.


title: What is the effect of having skewed dependent variable on scatterplot result

text: The histogram of my dependent variable is as following: ![enter image description here](http://i.stack.imgur.com/H6XlW.jpg) I draw the scatter plot of my dependent variable and independent variable, and the result is as following picture? I am wondering if the skewness of the dependent variable may affect the result. It seems that the relationship is not linear and we have a lot of outliers! ![enter image description here](http://i.stack.imgur.com/lEypk.jpg)


title: C-Index for cluster analysis in Stata

text: I'm wondering how to calculate the C-Index for determining a 'good' number of groups in a cluster analysis in Stata? It's mentioned in this post (What is an acceptable value of the Calinski & Harabasz (CH) criterion?) for R, but it seems Stata does not provide a built-in solution. Thank you!


title: manipulating arrays in SAS

text: English is not my first language. Here are my questions. I am using proportions ((0.1, 0.25, 0.25, 0.4) to generate 4 classes of data (1 variable 4 levels) in SAS based on the code below, * * *               data NativeBMI (drop = i) ;                          array prob [4] (0.1,0.25,0.25,0.4);           call streaminit(1234);                     do i = 1 to 100;                     BMI = rand("Table", of prob[*]);          output;                     end;            run;      * * * I have no problem with this code above. However, if I change the line from                   array prob [4] (0.1, 0.25, 0.25, 0.4);                                          call streaminit(1234);                         do i = 1 to 100;                         BMI = rand("Table", of prob[*]);      to                   array prob [4] z;         call streaminit(1234);                         do i = 1 to 100;                         BMI = rand("Table", of prob[*]);      it doesn't work. **My goal is to use a set of proportions saved in an array elsewhere (in this example z) and pass that array (z) into the rand function**                    rand("Table", of ...);      Thank you in advance.


title: F1 score for biased binomial data

text: I am applying a Bayesian classifier and would like to find out the f1 score. I determined the TP, TN, FP, TP. Unfortunately I had to find out that in my cross-validation almost in all test scenarios only TP and FP are non-zero (so TN=0, FN=0). Probably this is due to the fact that the ratio of the amount of data of the two classes is 80% to 20%. Of course I can calculate f1 score but as soon as I just swap the classes, I am not able to do so anymore because it would imply a division by zero for the $precision=\frac{TP}{TP+FP}$. Does anyone know a better way to calculate the f1 score or can anyone recommend me a better measure, preferably one that considers the bias among the two classes. Thanks for your help!


title: In R and INLA, how do I graph results of 'rw2' model?

text: I'm working with model 2 of INLA's Zambia example in R:               formula.mod2 = hazstd ~ f(agc, model = "rw2") +       f(district, model="besag", graph = g, param = c(1,0.01)) +       f(distr.unstruct, model="iid", param = c(1,0.01)) +       edu1 + edu2 + tpr + sex + bmi          mod2  =  inla(formula.mod2, data = Zambia, control.family = list(initial = 1),       control.inla = list(h = 1e-4),       control.compute = list(dic = TRUE,cpo=TRUE),       verbose = TRUE)      Everything works as advertised. In the Martino and Rue (2010) materials (section 5.2) I noticed nice graph of the 'Fig. 2a Estimated effect of age (posterior mean together with 2.5% and 97.5% quantiles) using MOD2'. ![enter image description here](http://i.stack.imgur.com/Bn13k.png) PS I asked that question on INLA user group; cross posting here to get more attention. My question is where do I find the data and how do I plot a graph like that?


title: Kalman Filter to correct model simulation bias

text: I am working with a large scale deterministic model, which attempts to simulate CO2 emissions in different regions. When compared to historic data, the model output suffers from systematic biases. Here are plots for two regions (model output in blue, historic data in red) : ![Plot 1](http://imgur.com/McCB6if.png) ![Plot 2](http://i.imgur.com/v0eCj0j.png) I am searching for a method to estimate the systematic bias in the model output and correct the forecasts of the model. In my view, the (scalar) Kalman Filter is an appropriate choice: $x_t = A x_{t-1} + B u_t + w_t$ $z_t = H x_t + v_t$ where $w_t$ and $v_t$ are gaussian white noise. The first equation represents the state equation, which is assumed to be unknown. The second one is the measurement equation. Here are my questions:   1. Can the kalman filter be applied to this kind of bias correction? If so, I am not sure about the interpretation of the two equations: Can I say, the historic time series is the state equation and the simulation is the observed measurement $z_t$? How is the bias estimated then?   2. I am searching for an implementation of this bias correction (preferable in R or MATLAB, pseudo code would be also very nice). Important is, that no special package is used (to my knowledge, one can implement the scalar Kalman Filter with a for loop). Thank you in advance! Here are the data for the second plot:               time    historic    model     1   0   0     2   0   0     3   0   0.323     4   0   0.323     5   0   0.323     6   0   0.299     7   0   0.299     8   0   0.299     9   0   0.303     10  0   0.299     11  0   0.252     12  0   0.372     13  0   0.298     14  0   0.292     15  0   0.398     16  0   0.432     17  0   0.435     18  0   0.424     19  0   0.407     20  0   0.404     21  0   0.406     22  0   0.378     23  0   0.327     24  0   0.353     25  0.239   0.375     26  0   0.356     27  0.229   0.408     28  0.28    0.394     29  0.251   0.409     30  0.268   0.411     31  0.281   0.41     32  0.264   0.401     33  0.286   0.402     34  0.257   0.412     35  0.284   0.4     36  0.276   0.417     37  0.273   0.367     38  0.28    0.39     39  0.278   0.431     40  0.305   0.437     41  0.334   0.445     42  0.309   0.461     43  0.306   0.458     44  0.289   0.458     45  0.288   0.454     46  0.264   0.427     47  0.283   0.473     48  0.28    0.459     49  0.308   0.426     50  0.393   0.432     51  0.326   0.418     52  0.292   0.406     53  0.257   0.412     54  0.232   0.409


title: beta differs from overall subgroup

text: My plain regression looks as follows: tendaystockreturn = alpha + beta * dummy variable that takes a value of 1 in case a director purchases shares Now I want to test whether a subgroup of directors will have a higher beta. In other words, I want to know whether the beta of a subgroup differs from the beta of the overall group. The overall director dummy variable has 52 out of 267 obseverations where it takes a value of 1. The subgroup director dummy variable has 17 out of 267 observation where it takes a value of 1. How should I put this in the regression? Should I make an interaction of both dummy variables or not? Or should I run both regression separately and then compare?


title: Linear regression of 0/1 response (Fig. 2.1 of The elements of statistical learning)

text: In chapter 2 ESL book authors write: Let's look at example of linear model in a classification context![enter image description here](http://i.stack.imgur.com/ANaTc.png) They fit a simple linear model $g = 0.3290614 -0.0226360\cdot x_1 + 0.2495983 \cdot x_2 + e$, where $g$ is given with values $0$ or $1$. Then they made a decision boundary where $\hat{y}$, if $\hat{y}>0.5$ then yellow. Question: There is a separation line on the picture. Where did intercept and slope for this line come from? In the `ElemStatLearn` R package, they simply put as `abline( (0.5-coef(x.mod)[1])/coef(x.mod)[3], -coef(x.mod)[2]/coef(x.mod)[3])`, where first term is the intercept, and second term is slope for this line


title: Statistical method to find capacity limits?

text: Im analyzing time-series to detect when the y-value is so flat that one can assume there is an underlying factor limiting y from being higher. Is there a methodology or statistical discipline that do this kind of analyis automatically? Im after a name I can put in Google that points me to books, articles etc. In the first example below, the eyes can easily detect that there is a cutoff at y=8. It looks like the curve should go higher but something is blocking it from doing that. The second curve, although having the same max of y=8, does not seem to have a bottleneck or limitation. I suppose our eyes and brain do a lot of post-processing, like looking at the variance before and after the plateau etc.               library(ggplot2)     the.data.1 <- data.frame(x.vals=seq(1:40),y.vals=c(1,2,3,2,3,2,3,4,5,4,5,6,6,5,6,7,8,8,8,8,7,8,8,8,7,8,8,8,7,6,7,6,6,5,4,5,4,3,2,1))     p <- ggplot(the.data.1, aes(x=x.vals, y=y.vals))     p <- p + geom_line()     p     the.data.2 <- data.frame(x.vals=seq(1:40),y.vals=c(1,2,3,2,3,2,3,4,5,4,5,6,6,5,6,7,8,7,6,5,6,6,7,8,7,8,7,6,5,6,7,6,6,5,4,5,4,3,2,1))     p <- ggplot(the.data.2, aes(x=x.vals, y=y.vals))     p <- p + geom_line()     p      Any suggestions for methods are welcome.


title: What's the resulting distribution when you take the N times self convolution of normal(mu, sigma) where N~Poisson or binomial

text: So if you have $N \sim \mathrm{Poisson}(\lambda)$ or $N \sim \mathrm{Binomial}(n,p)$ I'm curious about the convolution $$X = \sum_{i=1}^N X_i$$ (where $X_i \sim \mathrm{Normal}(\mu, \sigma)$ and $X=0$ if $N=0$) I guess alternatively this is the compound distribution $$\mathrm{Normal}(N\cdot\mu, \sqrt{N}\cdot\sigma)$$ but that's not much more help. Another way to think about this, if you know the number of customers entering a bank in a day is determined by a binomial distribution, and the amount they each deposit/withdraw is a normal distribution, what is the distribution for the amount deposited/withdrawn in a day (not just the EV)?


title: Obtaining easy-to-interpret decision tree

text: I am trying to create a decision-tree out of a number of attributes, where there are only two final classes and the classes are highly unbalanced (Class 1: 95.5%; Class 2: 4.5%). The idea is to profile the members of Class 2. For example, Class 2 members have attribute 1 >= 8, attribute 2 < 6, attribute 3 between 1/1/2013 and 12/31/2013. A profile would be successful if:   1. it is easily understandable by the business    2. it shows an improvement in predicting Class 2 over the trivial classifier; i.e.: those classified using the profile should have a > 4.5% probability of actually being in Class 2. Using WEKA's J48 algorithm along with its CostSensitiveClassifier on one attribute alone (treading carefully here) gives me the following tree: ![enter image description here](http://i.stack.imgur.com/J3b1w.png) Accuracy is 36.615%, which I can live with. 75% of Class 2 records were correctly classified, which is great. Most importantly, those classified as Class 2 had a 5.4% chance of actually being in Class 2, which meets goal (2) above. Here's the problem: that tree is _awfully_ complicated for just one attribute. I am ideally looking for a single rule of the type: attr >= x: Y; attr < x: N. I have tried many different weights in the cost matrix. My questions are:   * What techniques (different algorithms; different options in J48) can I try to obtain a simpler tree?   * Perhaps the data just doesn't admit a simpler tree. Can people speak of their experiences interpreting complex trees for non-technical audiences? I am really new to data mining, so if I've missed relevant information or asked the wrong question, please do let me know and I will correct. Thanks in advance!


title: power for estimating treatment heterogeneity in RCT

text: Could someone point me to a good resource for estimating power for detecting treatment heterogeneity via interactions in MR? Let's assume that individuals will be randomized to two study arms. To estimate heterogeneity (moderation), we would regress the continuous outcome on an indicator of assignment to treatment (0/1), a pre-treatment covariate (dichotomous or continuous), an interaction of the two, and possibly a vector of other baseline covariates, study strata, etc. I've read 'rules of thumb' about sample sizes for moderation needing to be 4x as large as the sample required to detect a certain main effect. Looking to estimate more precisely.


title: small sample size, large number of variables (most categorical) - how to proceed?

text: I would be grateful for general guidance/advice about data analysis with some data that is problematic for me because of the small sample size, and the large number of categorical data. I realize this question is a bit vague, but that's because I'm not sure what to do. I'd like to come up with some good descripive characterization of the data, and see if ther's any way to make any meaningful inferences. Sample size = 13, number of variables 93 (81 categorical, only 12 numerical). In the past I've worked with mostly numeric data with large sample sizes so I'm not sure how to proceed. Given the small sample, I don't feel I can make any assumptions about normality. Other than generating basic descriptive statistics (mean, std dev for the numeric data, and tables for the categorical data) what else can I do to meaningfully summarize this data? Is it possible to generate some reliable confidence intervals for the means via some nonparametric tests with such limited data? In addition to descriptive statistics, I am also wondering about simple linear regression. From what I have read, automatic selection of relevant variables, say via stepwise regression, is questionable to start with, and especially with such a mall sample won't be reliable. Also, I'm not so much interested in prediction, as to exploration of the relationships between a numeric response variable and the rest of the data. I also worry about detecting collinearity. With numeric data I could generate a correlation matrix, not sure this makes sense with this small sample, nor am I certain how to do the equivalent for the large set of categorical data. So not sure if there is a any automatic or semi-automatic way for stats to guide me to the relevant variables as a start and then take it from there. I.e., other than manually considering various combinations of independent variables, is there another way? I'm using R.