title: Derivation of sample autocovariance

text: The autocovariance is defined as $$\gamma(t,s) = Cov(X_{t}, X_{s})=E[(X_{t}-\mu_{t})(X_{s}-\mu_{s})]$$ When we have a stationary process the only thing that matters is the lag between the variables: $$\gamma_{k} = Cov(X_{t}, X_{t-k})=E[(X_{t}-\mu)(X_{t-k}-\mu)]$$ However, the expectation means that we are summing over all possible values of the random variable $X$. For example: $$\gamma_{k}=\int\cdot\cdot\cdot\int(X_{t}-\mu)(X_{t-k}-\mu)f(X_{t},...,X_{t-k})dX_{t},...,dX_{t-k}$$ Given this definition of $\gamma_{k}$, how can we derive the sample autocovariance: $$\gamma_{k} = \frac{1}{N}\sum_{t=0}^{N}(X_{t}-\mu)(X_{t-k}-\mu)$$ in which we are summing over $t$ instead of $X$. I think this is related to an ergodicity assuption, but I haven't found a very detailed explanation about it. Another question: Is this definition of sample autocovariance correct only when the process is stationary?


title: Reducing sample size to change from a paired to unpaired design

text: I have two groups of data, one consisting of 9 subjects and another consisting of 5 subjects. In the 9 subjects group, 5 of those could be paired with the 5 subjects from the other group. I could run a paired samples t-test and cut out the data for the 4 subjects that cannot be paired or I could run an independent samples t test and completely ignore the paired nature of the data. I think similar situations have been discussed on this site, but what I want to know is, if I go with the unpaired test instead of the paired test, would the only difference be a loss/change in power? Or in other words, are my results somehow invalid/ suspect if I don't take into account the paired nature of some of the data? Also I'm wondering if anyone has a good reference for what to do in this kind of situation. Thanks.


title: What software (paid or free) exists for learning large datasets?

text: Is there available software (or even just relevant papers) that can perform multiclass learning on datasets of 200m+ samples with 50+ classes and 1000+ features? What are the limits on dataset sizes for neural networks? Decision tree ensembles? SVM? As an example: Microsoft has developed code that can build decision trees on a 1kilonode cluster of 600m samples per tree over 32 classes with 2000 features. It takes a day to do the training of 3 trees. Are there publicly available programs that can do this for ANY of the above learning algorithms?


title: Modelling rainfall to size storage tanks

text: I have daily rainfall data for a given site going back about 30 years. I have a building with an average daily demand for water of $L$ litres and a catchment of $A$ m$^2$ with a runoff coefficient of $c$. How can I construct a model so that I can determine for a tank size T litres, The mean and variance of the   1. days per year the tank will be empty    2. Shortfall of water   3. days per year the tank will overflow   4. Amount of water that overflows I have googled this extensively and found plenty of rules of thumb but noting that uses a probabilistic model.


title: Heteroscedastic test does not solve the problem

text: As my previous questions I'm trying to solve a problem with my stocks tests. I tried Breusch-Pagan test for heteroscedasticity but some residuals still pass these tests. My procedure is:   * Get two stocks prices (I have a matrix with two columns that represent the price lists)   * I do a linear regression, like: `lm(prices[,1] ~ prices[,2])`   * Then I test the residuals of the linear regression with Unit Root tests (PP and KPSS)   * After these tests I do Breusch-Pagan test because I only need to work with stocks that have constant variance during all the period. I do: `bptest(prices[,1] ~ prices[,2])` Ok, now I still get strange results, with strange i mean that when i plot the residuals i see that the variance is not constant (take a look at the chart below). So now, I need to understand how better test is the variance is constant. I read someone use GARCH (1,1) but I never used it, could someone exmplain it or maybe give me other tests to try? ![enter image description here](http://i.stack.imgur.com/OYqKq.png)


title: How is Hyndman's explanation of proper Time Series Cross Validation different from Leave-One-Out?

text: Hyndman's great explanation of proper time series CV is at the bottom of the page in the following link: http://robjhyndman.com/hyndsight/crossvalidation/ Leave-One-Out illustration in the following link: http://i.imgur.com/qrQI4LY.png In the 'leave one out' illustration if the dataset in the illustration was time series data and was sorted from past to present from left to right, wouldn't it be identical to Hyndman's explanation on Time Series CV? If not, how so?


title: vcovHC, vcovHAC, NeweyWest – which function to use?

text: I am trying to update my lm() based model to get correct standard errors and tests. I am really confused which VC matrix to use. The `sandwich` package offers `vcovHC` , `vcovHAC` and `NeweyWest`. While the former only accounts for heteroskedasticity the latter two account for both serial correlation and heteroskedasticity. Yet, the documentation does not tell much about the difference between the latter two (at least I don't get it). Looking to the function itself I realized that NeweyWest actually calls vcovHAC. Empirically the results of `coeftest(mymodel, vcov. = vcovHAC)` and `coeftest(mymodel, vcov. = NeweyWest)`are mad different. While `vcovHAC`is somewhat close to the naive lm results, using NeweyWest all coefficients turn insignificant (tests even close to 1).


title: How to combine likelihoods from multiple binomial trials. Can we also determine the overall FDR?

text: A king has rounded up 1000 people suspected of counterfeiting coins, which look and feel the same as the official coin. However, only the official coin is truly fair (Pr(heads)=0.5), while all fake coins produce strongly skewed coin flip results (biased towards heads or tails). The king decides to identify the counterfeiters by flipping a coin taken from each of the suspects. He leaves the actual flipping to 6 of his guards, handing them a suspected fake coin **AND** one official coin. The latter is to assure the king that his guards carried out the coin flips faithfully. **Unfortunately** , the king forgets to tell the guards how often to flip each coin, only that both coins (fake and official) must be flipped an **equal number** of times. All 1000 suspected fake coins tested in this manner go through the hands of all 6 guards. Because heads of proven counterfeiters will roll, the king cannot afford to be wrong and commit too many innocent people to die. But he's comfortable executing up to 50 innocent people. How should the king analyze the coin flip data from his 6 guards to identify the counterfeiters among the 1000 suspects? Can he also find a way to keep the number of innocently convicted under 50? To get started, here's a sample from the first suspect (heads/tails):               Guard 1: suspect coin (88/11), official coin (42/47)     Guard 2: suspect coin (38/5), official coin (22/21)     Guard 3: suspect coin (115/15), official coin (70/60)     Guard 4: suspect coin (39/33), official coin (35/37)     Guard 5: suspect coin (70/13), official coin (43/40)     Guard 6: suspect coin (22/18), official coin (19/21)      Thank you!


title: How to figure out what numbers often appear together in a dataset?

text: I have a lottery style dataset we produce internally (example below). I am trying to figure out which numbers appear most frequently together. Example questions: What are the top 10 pair of numbers that appear most frequently together? What are the top 10 three numbers that appear most frequently together? What methods/techniques would I need to use to answer these questions? I was thinking clustering at first but a test run (ward) I performed, didn't return the results I was expecting. ![enter image description here](http://i.stack.imgur.com/Cc06x.png)


title: Z score standardization of 3 patient measurement scales

text: I am involved with a survey where 3 different patient measurement scoring systems were used. They all are used as surrogate outcome measured for patient self- management. All scoring systems are on a different scale - that is one ranges from 0-100, another 1-8 and other 1-10. The research wants to z-score each measurement and then perform ICC calculations between the measures (2 at a time) and generate bland Altman plots. The standardization results in the same correlation results but is it appropriate to perform Bland Altman plots and ICC statistics? The measurements are only taking once per subject - so each individual as 3 variables , one for each question - only measured once. Thank you for any help/suggestions KJ


title: Anomaly detection with one-class SVM

text: As I understood from my experimentation on some multivariate data-sets, anomaly detection (AD) heavily depends on the underlying distribution of data. Like, for ex., you can devise a method to detect anomalies supposing a Gaussian distribution for your data, and then end up with a failure if it is actually a power-law... Then a visual investigation followed by a "fit a model" approach seems to be the best way for AD. And I started to think that ML algorithms like one-class SVM is useless in most situations for AD. So here is my silly and novice question: "Suppose you have a data-set with many variables (IDs) distributed in many ways (some gaussian, some power-law). How one can use a one-class SVM approach in this situation without taking care of distributions"? Or do I need to study SVMs in more detail?


title: Statistics of a simple coin toss

text: Original simple question: "Given a fair coin that has been tossed 100 times, each time landing heads, would it be more likely that that the next coin flip be tails or heads" The proportion of head to tails is 1 to 1 as number of trails approach infinity (this is a fact correct?). Therefore, if this must be the case, mustn't there be an enacting "force" per say that causes this state of being to (admitted unreachable, but _technically eventual_ ) case? Hence if there exist a force pushing the number of heads and tails towards one another, shouldn't there exist a higher chance of tails than (regardless of how significant) heads? Or is this thought process just illegitimate simple because the state exist on at the infinity case?


title: Heavy-tailed distribution with closed-form ML fit from data

text: Which (if any) heavy-tailed distributions can we compute the maximum likelihood parameters of, given some data to fit the distribution to?


title: How to apply the random statement in my model?

text: i am trying to model the following data in SAS using proc glimmix and i would like feedback if i have modelled my data correctly. My data consist of individual chickens (ID) who are grouped by treatment in the enclosure they occupy (chickens receiving different treatments are kept in the same enclosure). I want to know if treatment influences if they appear at popholes present in their enclosure (registered, in which '1' describes an event). Furthermore if two weather variables (Factor1, Factor2) may influence this appearance. Repeated measurements (maximum of once a day) have been taken of the same individual chickens over a course of days (Date). I have attempted the following models but i am unsure which to use specifically since i am fairly new to using both proc glimmix and sas. Model 1               proc glimmix data=chicken ;     class date id treatment;     model registered (event='1') = treatment|factor1|factor2 / link=logit dist=binomial solution;     random intercept/subject=id(treatment);     Lsmeans treatment/pdiff;      run ;      or Model 2               proc glimmix data=chickens ;     class date id treatment;     model registered (event='1') =  treatment|factor1|factor2   / link=logit dist=binomial solution;     random treatment/subject= id;     Lsmeans treatment/pdiff;      run ;      Furthermore, in order to model repeated measurements on the same individual i have attempted to enter in both models the following code:               random date/ residual subject=id type=ar(1);         though then either model would only converge if pconv=1e-4. So my two main concerns are which model would be correct to use (with a focus on both random statements), and should i include the extra random statement. Furthermore, for statistics involving either weather factor the DF are increadibly high. observations have been made for about 500 treated chickens measured once a day over a period of nearly a month and while treatment DF=479 for factor1 and factor2 the DF= >12000. I hope my description is clear enough and thank you for your time in advance. Marielle


title: Chi-square test for one group

text: I'm struggling with this problem. For the analysis of my test, I was instructed to use chi-square test for one group. I use Statistica and everywhere I try to do this I need the expected values, which I don't have (or I don't know how to prepare it). The analysis is crucial for my PhD. In the test I have a question:               Which economic indicator and how much influences on the choice of the heating system in:                                |A.Single family buildings | B. Multifamily  | C. Public buildings      1. NPV               |                          |                 |     2. IRR               |                          |                 |     3. SPBT              |                          |                 |     4. Investment costs  |                          |                 |          Replies: 1 (no effect) to 5 (significant effect).       The respondents had to enter answers for each type of building, for a total of 3x4 = 12 numbers. As a result I need to know which economic indicator in each group (type of building) and how strong influence on the choice of heating system. The analysis can be performed separately for each type of building... For any idea how to solve this, I would be very grateful !


title: ARIMA - is it normal to be unable to recover DGP of simulated data with 200,000 datapoints?

text: When learning how to use a statistical model, I often create a simulation of the DGP, and run simulated data through the software to make sure I know what I'm looking at. I've done this for Proc ARIMA in SAS, and for the most part the results are expected. One thing I'm finding surprising is that if I overparamatise the model, then the estimated coefficients do not necessarily correspond to the DGP. I'd like to know whether this is to be expected. An example will clarify. Suppose I simulate data according to the following model (say 100,000 datapoints): $y_t = e_t + 0.9e_{t-1}$ In SAS I then estimate a model, allowing it to estimate 3 MA and 3 AR terms:               proc arima data=ts;     identify var=y (0) nlag=8;     estimate q=3 p=3;      run;      When I do this, often there are statistically significant AR terms, and the MA terms at lags 2 and 3 are also significant: ![ARIMA outputs](http://i.stack.imgur.com/JQWny.png) Is this because there are many different ways of writing the same data generating process? (e.g. any AR process has an equivalent MA representation). Or is something else going wrong? Thanks very much for any help! I note that from the ACF and PACF it's clear it's an MA(1), and if I estimate the model with q=1 and p=0 (i.e. a single MA term), it gives me the right results.


title: Interpreting the marginal effects in terms of standard deviations

text: I have a very simple logistic regression in which the binary variable Y is regressed on three continues variables, X1 and logX2 and X1*LogX2. X1 is a proportion (is between 0 and 1), logX2 is the natural logarithm of the continues variable X2 and X1*LogX2 is the interaction of these two variables. I have calculated the marginal effects of X1, LogX2 and X1*LogX2 at their mean. How would you interpret the marginal effects in terms of the standard deviation of the actual variables? lets say the following is the results of the marginal effects estimations and the standard deviation of the three variables. Variable: X1 , LogX2 , X1*LogX2 , X2 Mar.Eff.: 9.34744 , 3.72533 , -1.07857 , --- Std. Dev.: 0.05845 , 0.17723, 0.80012, 0.52 Every unit of increase in X1 increases the probability of Y by 9.34744 (930% ?) . Since X1 is proportion I can say that every percent of X1 increases the probability by 0.0934 (or 9.3%). Does this mean that if X1 increases by 1 standard deviation (0.05845) then the probability increases by 0.546 (9.34744*0.05845)? For LogX2, I can say that every 1 percent of increase in X2 results in 0.037 (3.72533/100) increase in probability. To interpret the results in terms of the standard deviation, can I just multiply the marginal effects of LogX2 by the standard deviation of X2 and say that every standard deviation of increase in X2 increases the probability by 0.019 (0.52*3.72533/100)? It does not seem right, because the changes in X2 should be studied in percentage, not the mere units! What is your suggestion? also How would you interpret the interaction term?


title: Multivariate hypergeometric confidence intervals

text: Suppose you have an urn filled with money. You know there are $N$ total bills, and that each is either a \$1, \$10, or \$100. You draw $n$ bills without replacement from the urn, and wish to construct, say, a confidence lower bound for the total amount of money in the urn. The obvious solution is if $n$ is sufficiently large, then you can use the normal approximation - but let's assume that $n$ is small and/or the distribution of bills is highly skewed, so that we don't expect anything resembling decent coverage from CLT-based confidence intervals. Are there any effective exact methods for computing confidence bounds for these sorts of problems? If you have only two types of bills, then it's not difficult (e.g. Buonaccorsi). If you know that each bill is either a \$1 or \$10, and your sample returns $k$ \$10s and ($n-k$) \$1s, then an $\alpha$-confidence lower bound for the total number of \$10s is the smallest $K$ satisfying: $P$(draw at least $k$ out of $n$ \$10s | $K$ out of $N$ are \$10s) > $\alpha$. Of course, in this case, having a lower bound for the total number of $10s is equivalent to having a lower bound for the total amount of money - but this no longer holds when you have more than two types of bills.


title: Estimating the running time of gbm grid tuning

text: I am trying to estimate the time it will take me to run a tuning grid on gbm (I am using the R Caret package but this is irrelevant as I am interested in the relative processing time). I can see many parameters:   * number of rows and columns of the data set,   * number of resampling (I am using method = "repeatedcv", number = 4, repeats = 2 for example)   * the tuning grid dimensions: interaction.depth, n.trees, shrinkage but I have a feeling the interaction.depth might have a non-linear impact... In any case is there a way to estimate the time it will take depending on all parameters? (of course it depends on hardware, but I will compare on small data set and grid on same hardware).


title: Feature learning with a deep learning aproach?

text: How to create a feature vector from text with a deep learning aproach?. Im new at this topic, could anybody advice me where to start and how to aproach this task?.


title: Resources for developing prognostic index-scores

text: I am looking for online resources or books explaining how to develop prognostic indexes-scores step by step. I am mainly wondering about model validation and transformation of regression coefficients into scoring systems.


title: Ancillary statistic: $X_i \sim N(\theta, \theta^2)$

text: Let $X_1, X_2, ... , X_n$ i.i.d random variables with probability density function $N(\theta, \theta^2)$. Show that $$T(X) = \frac{X_{(1)}-X_{(n)}}{X_{(2)}-X_{(n)}}$$ is ancillary to $ \theta$. My attempt: Since $\frac{X_i}{\theta} \sim N(1,1)$, which does not depends on theta, we have that if $$Y_i=\frac{X_i}{\theta}$$ the $ Y_i $ are iid $N(1,1)$ and, for example, $$Y_{(n)} = 1/\theta X_{(N)}$$ The same occurs to the other order statistics. Then, $$\frac{X_{(1)}-X_{(n)}}{X_{(2)}-X_{(n)}} =\frac{X_{(1)}/\theta-X_{(n)}/\theta}{X_{(2)}/\theta-X_{(n)}/\theta}$$ Does not depends on $\theta$. Is that correct? Thanks!


title: At what VIF level should you switch from OLS to ridge-regression?

text: Regarding multicollinearity, is it recommended to use ridge-regression if you have some covariates with VIF values around 10 in the OLS model? What would be the best VIF level to use to decide whether or not to do ridge-regression?


title: Double exponential smoothing in multivariate multilevel panel regression

text: I would like to use double exponential smoothing to predict prevalence rates of care dependency in Austrian federal states. My data is very detailed, thus I would like to make use of that in order to refine my predictions. I have the percentage of people in care dependency levels 1–7 aged 50–99 in 9 Austrian federal states.                str(daten[1:12][daten$jahr>1996,])     'data.frame':   39600 obs. of  12 variables:      $ age       : num  50 51 52 53 54 55 56 57 58 59 ...      $ gender    : Factor w/ 2 levels "male","female": 1 1 1 1 1 1 1 1 1 1 ...      $ bundesland: Factor w/ 9 levels "Bgld","Ktn","Noe",..: 1 1 1 1 1 1 1 1 1 1 ...      $ jahr      : num  1997 1997 1997 1997 1997 ...      $ PfSt0     : num  0.992 0.989 0.985 0.985 0.985 ...      $ PfSt1     : num  0.001458 0.000967 0.001459 0 0.002199 ...      $ PfSt2     : num  0.00437 0.00193 0.00802 0.00793 0.00587 ...      $ PfSt3     : num  0.00146 0.0058 0.00073 0.00433 0.0044 ...      $ PfSt4     : num  0.000729 0 0.002188 0.002163 0.000733 ...      $ PfSt5     : num  0 0.000967 0.002188 0.000721 0.002199 ...      $ PfSt6     : num  0 0.000967 0 0 0 ...      $ PfSt7     : num  0 0 0.00073 0 0 ...      DES is a time series analysis method. Time series analysis explains a data series by its past values only. While it is true that I use only past data of care dependency, one could regard age, gender and federal state as explanatory variables. Instead of computing individual double exponential smoothing forecasts for each age, gender, federal state combination, I could assume structural uniformity within these time series. Thus, my data might be regarded a multilevel panel dataset, with 50 observations per year (age groups) nested in 9 federal states each. (I plan to do separate regressions for males and females.) I would like to use federal state, age and age squared as explanatory variables apart from previous value and previous trend, as done in double exponential smoothing. However, in panel data analysis, time trends are typically covered by including the year variable in the regression, and rarely ever by including lags. **How could I realize a forcasting method similar to double exponential smoothing in a panel dataset, i.e. including also other explanatory variables?** (Preferably in R) (Matters are complicated further by the fact that I have 7 instead of 1 dependent variables.)


title: What are possible reasons for differences in numerical output from the same procedures performed in different statistical packages?

text: My colleague and I do the exactly same thing with the exactly same data and get different results. He uses Stata, I use R. We both use 32-bit Windows machines. Are there any possible explanations apart from human failure for this? The concrete problem is calculating c-indices for Cox models. His are up to 0.005 higher.


title: Concepts behind fixed/random effects models

text: 1. Can someone help me to understand fixed/random effect models? You may either explain in your own way if you have digested these concepts or direct me to the resource (book, notes, website) with specific address (page number, chapter etc) so that I can learn them without any confusion.    2. Is this true: "We have fixed effects in general and random effects are specific cases"? I would especially be grateful to get help where the description goes from general models to specific ones with fixed and random effects


title: How to build a predictive model using sound (frequency and amplitude) as the input

text: I'm working on a medical research project where we're looking at the changing resonances in femur bones during hip replacement operations. Some surgeons use the increasing pitch of the bone resonance to tell them when to stop hammering/reaming the shaft of the feamur. If they apply too much pressure when creating the cavity for the replacement hip, they risk fracturing the femur. Basically, I've recorded a number of these changing resonances using a bone conduction mic and now I have a number spreadsheets with one column listing frequencies in 5Hz bands and 3 columns showing the amplitude of that frequency when the bone is loose, medium to tight (referring to the tension of the reamer/chisel creating a cavity for the replacement joint). I'm coming at this more from a sound perspective so I would really appreciate some advice from anyone with a handle on stats and predictive modelling, I'm trying to find a pattern to the change in frequency which is consistent in all cases. Any advice would be very much appreciated.


title: What's the difference between dlmSmooth and dlmFilter in R's dlm package?

text: Could someone please explain what the difference is between the two, and perhaps avoid the worst statistical jargon? I am currently using the `dlm` package to model dynamic regressions as can be seen on p.122-5 in Dynamic Linear Models with R, and I don't really grasp the theoretical difference between using `dlmFilter` and `dlmSmooth`. More specifically, why is the `dlmFilter` example on p. 123 not considered dynamic, while the `dlmSmooth` example on p.124-5 is? And also, why are the two results fairly different?


title: How do I model repeated measure + random effects with SPSS mixed

text: I'm having trouble understanding what a correct model for my data should be, and also if it would even be a good idea to make it a single model. At the basis it is a rather simple regression model with two predictors. The predictors are 'perceived ease of use' (ease) and 'Perceived Usefulness' (useful), and the dependent variable is 'Intention to Use' (intention). Now the thing is: we measured these variables for 10 different technologies. So for each person we have a multiple of these variables. On top of that we have two measurements for each person (repeated measures) If at all possible, I would like to fit this all into one model (mixed model in SPSS). Ofcourse an 'easy' solution would be to do a 'simple' repeated measure model for each technology separately, but then I think you lose a lot of information, cause you cannot partial out the effect of the person, and difference between the models would be hard to explain. So I guess my first question is, is it valid to treat all dependent variables (intention) for the different technologies as the same outcome variable? or do we _have_ to do 10 separate analyses? As I said the reasons to make a single model is that more information is preserved (the effect of the person), thus hopefully resulting in a more general model that is based on 'more' data. And if it is valid to do this, I'm not sure how to model this in SPSS. Right now I have this:               MIXED intention BY personID time technology WITH  ease usefulness     /FIXED = usefulness ease | SSTYPE(3)     /RANDOM=technology time|  SUBJECT(personID) COVTYPE(ID)     /REPEATED = time*technology | SUBJECT(personID)     /PRINT = SOLUTION TESTCOV.      However I'm not sure this models the dependencies in the data correctly. I also don't know what kind of covariance matrices to choose. Can anyone tell me if this model fits my description well, or what I can do to improve it? Thank you very much!


title: Kalman smoothing of returns vs. prices with dlmSmooth in R's dlm package?

text: So I am using the R code behind Fig. 3.14 in Dynamic Linear Models With R (p. 124-5) to make a dynamic version of a simple pair trading model: $$ Y = \alpha + \beta X. $$ If I use log returns (`diff(log(P))` in R, P being a timeseries of equity prices) I get results that makes sense, they are pretty similar to the static regression of the whole dataset (comparable to the R code on p. 123). So far so good. But if I instead use log prices (`log(P)` in R, P being a timeseries of equity prices), the results don't make sense at all, $\alpha$ has the wrong sign and $\beta$ is almost constant throughout the timeseries at half of what I would expect compared to the regular static regression. So what am I doing wrong here? How can I change the code to produce better results for log prices?


title: Find trendline for minimum (not mean) values in distribution

text: I would like to perform something like a linear regression on my distribution of data, but I'm interested in a trendline that estimates the **minimum** , _not mean_ , value for each time bin. I'd like to do this in R. The image below shows a scatterplot of the minimum value for each time bin. The black line is a typical linear regression, which estimates the mean. What I'd like is something like what I painted in red - an estimation of the minimum. ![enter image description here](http://i.stack.imgur.com/D0O7Q.gif) My data look like this: ![enter image description here](http://i.stack.imgur.com/qGAiS.png) Those are just the first few lines but you get the idea. Thank you!


title: R sparse smoothing splines

text: I would like to fit a model of the form `weight ~ bs(height) + id`, where id is a factor (so I have several observations per every id). There are perhaps 100k values of `id`. Does anyone know the best way to go about fitting this? I have tried the mgcv package, but I do not think it was designed for problems this large and sparse. I can define my own knots and figure out the best smoothing parameter, if necessary. Also, if anyone has a solution in Python, that's great too. Or just the right phrase to google would be helpful, also.


title: Target motion analysis with support vector machines

text: Does anyone have experience or knows a good source where SVM are used to perform target motion analysis? > A general maneuver tactic against submarines is a zig-zag. A submarine > usually relies on passive detection, not risking active sonar or a periscope > observation. So to determine where a unit is heading the submarine needs > target motion analysis (TMA). This requires several minutes of passive > contact and if the contact starts to zig-zag this process must restart What I a interested most how the submarines over passed the zig zag maneuver and effectively applied TMA. I searched over the net but no resources could be found except the classic TMA.


title: Analysing a target based trial

text: I have a trial where the measurement is the deviation about a target value. There are 8 subjects - 4 treatment, 4 control. Each subject has undertaken the trial 4 times. Thus, I have 32 values, nested within groups. Obviously this is a rather small number Essentially, my treatment is effective if it minimises the absolute values of the deviations (ie 10 above the target value is the same as 10 below). Can I ask for any suggestions on how this data should be analysed, particularly how I should account for the nesting?


title: Unexpected intercept in linear regression with dummy variables

text: I am modeling a binary DV with a factor variable consisting of 8 clasess. I know I could just analyse the distribution of the DV across all factor classes but I will add more variables and wanted to start the predictive model with a regression approach. The training set is :               > head(learn_set[,c(1,12)])       learn_IV factor     1        0      1     2        1      5     3        0      2     4        0      2     5        0      3      Just to make sure that the DV is binary :               >table(learn_IV)      learn_IV      0      1       102970  44323      And the distribution of the IV factor variable:               >table(factor)     factor         1     2     3     4     5     6     7     8      47048 33975 24982 30023 19022   668  7655   286       Now let's see the call of the `lm()` function               >summary(LM)          Call:     lm(formula = as.numeric(learn_IV) ~ factor, data = learn_set)          Residuals:      Min      1Q  Median      3Q     Max      -0.6905 -0.1200 -0.1198  0.4518  0.8917           Coefficients:                  Estimate Std. Error  t value Pr(>|t|)         (Intercept)  1.548154   0.001997  775.299   <2e-16 ***     factor2     -0.428366   0.003083 -138.962   <2e-16 ***     factor3     -0.439816   0.003393 -129.607   <2e-16 ***     factor4     -0.428185   0.003202 -133.718   <2e-16 ***     factor5     -0.083549   0.003721  -22.454   <2e-16 ***     factor6      0.142362   0.016884    8.432   <2e-16 ***     factor7     -0.061965   0.005331  -11.624   <2e-16 ***     factor8     -0.433530   0.025916  -16.728   <2e-16 ***      How come the intercept is `1.54` when all variables are binary? All the IV's are mutually exclusive so this model would always predict a value greater than `1` even for in-sample testing. I have seen in another thread that a constant `-1` was added in the formula but this doesn't make much sense to me. Has this something to do with the binary variable trap?


title: Using different variables for each dependent value in multinominal logistic regression

text: I am not sure if this is possible, but my question is - Is it possible to use different sets of independent variables in the same model for different values in a multinomial logistic regression model? In other words, I have a set of variables, which are applicable to only one dependent value and not to other values of my DV. So, I don't want to use this variable for other observations. Is this possible?


title: Calculating $R^2$ when one variable can only take integer values

text: I am a programmer, not a statistician, so pardon my botched use of the terms. My basic problem is this: I am wanting to calculate $R^2$ between a known concentration (which can be any non-negative value) and a discrete measurement (where the values are all integers). There are 92 possible observations, and each of these 92 has a known actual concentration. We wish to measure the accuracy of the measurements by looking at the $R^2$. The different cases (each a different molecule) have concentrations which vary across orders of magnitude, so we use the log (base 2) of the values when calculating the $R^2$. This is an industry standard convention for this case. In most cases, this works well. However, in cases where the measurements are relatively low, I am thinking that it may cause an artificially low $R^2$. For example, if the known concentration should result in a measurement of 0.1, since the only measurements possible are 0 or 1, then it will interpret this as error. Even if there are around 10 cases with concentrations that should give a measurement around 0.1, and we get one detection of one of them and 0 for the other 9, this will be interpreted as error. For higher concentrations, this is obviously less of an issue (if the detections should have been 118.5 and we got either 118 or 119, it will correctly interpret this as not much error). However, I don't want to just make up my own correction for this. My guess is there is some standard way of handling the calculation of $R^2$ between a continuous and a discrete variable. Can you point me at it? I'm doing my calculations in Python using the scipy.stats module, but if you just know the name of the proper calculation method and don't know the python code that's perfectly fine. p.s. To be more clear, there are 92 molecules of known concentration, and we are measuring their concentration using a technique. We want to know if a given measurement run went ok, and so the $R^2$ of the measurement run (which is discrete, i.e. how many counts do you have for that molecule) vs. the known true concentration (which is continuous) is being used to determine how this run's accuracy compares to other runs (for the exact same set of molecules). Hopefully the fact that it is always the same 92 x axis values (where y is the measurement), and we are comparing only one measurement run of this type to another, makes $R^2$ not too bad a metric to use here.


title: Can someone please explain the back-propagation algorithm?

text: What is the back-propagation algorithm and how does it work?


title: Significance of most frequent answer

text: I have the following data: I have participants who answered on their favorite out of 6 positions of various information types on a display (e.g. 1 = bottom left, 2 = top left etc.). I can now say for example the information “CPU status” was preferred to be shown on the top left by most of the participants, i.e. I can report the frequencies of the answers. But how can I test if these frequencies are not due to sampling or chance? In other words, how can i say that position a was significantly more often chosen than position b and c? My only idea is to do a chi squared test with the two most frequent positions per information, but this seems cumbersome to me.


title: Why do US and UK Schools Teach Different methods of Calculating the Standard Deviation?

text: As I understand UK Schools teach that the Standard Deviation is found using: ![alt text](http://upload.wikimedia.org/math/e/e/4/ee485814ab9e19908f2b39d5d70406d5.png) whereas US Schools teach: ![alt text](http://upload.wikimedia.org/math/8/3/a/83a7338b851dcaafaf5b64353af56596.png) (at a basic level anyway). This has caused a number of my students problems in the past as they have searched on the Internet, but found the wrong explanation. Why the difference? With simple datasets say 10 values, what degree of error will there be if the wrong method is applied (eg in an exam)?


title: Measuring the effectiveness of a pattern recognition software

text: Are there any objective methods of assessment or standardized tests available to measure the effectiveness of a software that does pattern recognition?


title: What is meant by a "random variable"?

text: What do they mean when they say "random variable"?


title: Hidden extrapolations in linear regression

text: On page 261 of Kutner, the author cautions readers about hidden extrapolations in linear regression. As shown in the figure, the allowed region for prediction is not a rectangle delimited by the ranges of $x_1$ and $x_2$, but a ellipse. Why is this the case and how is this ellipse generated? ![enter image description here](http://i.stack.imgur.com/LCZTP.png)


title: How to compute the free energy of a RBM given its energy?

text: The free energy is defined as such: $F(v)=-ln(\sum_he^{-E(v,h)})$ What I don't understand is how to compute this $\sum_h$ ? If I know $E(v,h)$, how can I compute $F(v)$. I have seen several definitions of $F(v)$ for a binary-binary RBM, but I have not been able to work it out from $E(v,h)$ of the RBM. And can I do the same with gaussian visible units (changing $E(v,h)$) of course) ? And what about ReLU units ?


title: How would I calculate statistical significance for a range?

text: I am running an A/B test, and the following are samples of how many questions a user has answered. I am trying to figure out which test is better (A or B), and how confident we are that it is better. For example, in Sample A, first user answered 4 questions, next two users answered 5 questions. I know how to calculate statistical significance for A/B tests, but am not sure how to do it for a range of numbers. I am trying to get users to answers the maximum questions possible. Sorry for my incorrect lingo but i'm not a stat guy, so i'm not sure what you refer to these terms. Essentially, I would like to figure out it if my A/B test, A sample or B sample yields a higher number.               Sample A :     4, 5, 5, 9, 11, 14, 15, 15, 16, 19, 27, 30, 31, 32, 58, 65, 67, 79,     98, 99, 100, 103, 106, 204, 232, 341, 354, 359, 360     Med: 58     Avg: 98          Sample B:     1, 3, 4, 4, 4, 5, 9, 11, 12, 15, 19, 28, 37, 48, 50, 54, 59, 72, 74,     78, 80, 81, 89, 91, 99, 101, 103, 103, 104, 120, 121, 174, 203     Med: 59     Avg: 62


title: Bayesian updating of continuous variables given mutual information

text: I have very little stats training, so this may be a very obvious and boring question, in which case I apologise. Given two real-valued continuous random variables A and B, and given prior probability distributions $f_A$ and $f_B$ for each, how do I perform a Bayesian update when given the extra evidence that b>a with probability p? Thanks in advance.


title: check if the numbers belong to the same cluster

text: I would like to find if the particular two numbers in a vector belong to the same or different groups, For instance, I have a set of vector data: $$d_1=\{100, 120, 104, 110, 20, 25, 30, 35\}$$ I want to know whether $100$ and $35$ belong to the same group or not. By just looking at this, it is found that this set of data can generally be clustered into 2 groups $g_1=\{100,120,104,110\}$ and $g_2=\{20,25,30,35\}$, thus the probability that $100$ and $35$ belong to the same group is lower as compared to another example: $$d_2=\{100 ,90, 80, 65, 60 ,50 ,40 ,35\}$$ This set doesn't have a clear cut between the numbers (they decrease gradually). Is there a formula to calculate the 'score' (e.g.: $s$ for $d_2$ is higher (near to 1), as $100$ and $35$ belong to the same group and vice-versa) based on the examples given? I thought of utilizing the mean and standard deviation, such as $$a(d)= \left( \left(\frac{100+35}{2}\right)- \frac{\mbox{mean}(d)}{\mbox{std}(d)} \right)$$ and $$s(d)=\exp(-c\times a(d))$$ where $c$ is a constant. But I am not sure if this is a valid way of computing it, does it make sense? Please advise. Or is there any better way to model the function?


title: AR(1) Covariance

text: So I'm trying to derive the covariance between $z_t$ and $z_{t-1}$ in the $AR(1)$ model: $z_t =\phi z_{t-1} + a_t$. Can anyone give me some advice on where to start?


title: Regression - significance of a variable

text: I have been running a linear regression, and one of the variables was a discrete variable (3 categories), decoded into three binary variables. I ran the regression in Stata. The first of those three variables was omitted, the other two were tested for significance (a ttest I assume); the resulting p-values were strikingly different: $0.682$ and $0.049$. Before removing entirely that discrete variable from regression, I decided to test those two binary variables for joint significance; the result was a p-value of $0.0746$, which would indicate that for a significance level of, say $\alpha=0.1$ (still reasonable, I guess), we reject the null asserting joint insignificance. I know that in general, there is no "recipe" which would directly tell us which variables are significant, and which should not be kept. Nonetheless, in this case, is there any sound rationale for keeping the variable?


title: Use squared correlation in regression without intercept

text: If I want to compare the goodness-of-fit of two regression models, with and without intercept, is it valid to compare the squared correlation coefficient between the fitted values and the data? Since the squared correlation would get back the $R^2$ for the model with intercept, it seems to make sense to compute the squared correlation for the model without intercept and use that to make comparison. I am not entirely sure whether this is legitimate. Am I missing anything? If this approach is not valid, what methods can I use to make such comparisons?


title: How do you calculate the fastest rate at which a disease increases globally?

text: I have a q and for with the answer, but cannot get the same. The q is based in table form with a list of major diseases around the world and the date collected for three years total deaths x1000. for each year. i.e: TUBERCULOSIS IN 1999 WAS 1669 DEATHS(x1000), In 2000 was 1660 (x1000) and in 2001 was 1644(x1000). How would I calculate the four diseases increasing globally at the fastest rate? Here the answer is 150%


title: predict.glmnet test data

text: I used `glmnet` to build a predictive model with ~200 predictors and 100 samples for a binomial regression/classification problem. This was my training data. I selected the best model (16 predictors) that gave me the maximum AUC. I have an independent test data set with only those variables (16 predictors) which made it into the final model from the training data. I was wondering if there is any way to use the `predict.glmnet` based on the optimal model from the training data with the new test data set, which has data for only those variables that made it into the final model from the training data. Any suggestions would be appreciated.


title: How to cross-classify two categorical variables?

text: I have a data frame with two categorical variables (representing different methods) that "classify" objects into three categories: `cat1`, `cat2`, `cat3`. The objects are classified somewhat differently by the two variables and I would like to perform a test described as "Measuring Agreement Between Observers" in Agresti's CDA, Section 10.5 (after Thompson 2009, p. 192) to check if both methods perform differently. Wanting to follow the example in package exactLoglinTest, I struggle to create a cross-classification data frame similar to `library(exactLoglinTest); data(pathologist.dat)` My variables are in the following form:               classification.1    classification.2     cat1                cat1     cat1                cat2     cat3                cat3     cat2                cat2     cat2                cat3      ...etc I have more than thousand classified observations and (as far as I understand) I would like to count all the combinations - i.e. `cat1` & `cat1`, `cat1` & `cat2`, `cat1` & `cat3`, etc. How can I achieve this? References: Thompson, L.A. 2009. R (and S-PLUS) Manual to Accompany Agresti’s Categorical Data Analysis (2002) 2nd edition. available online


title: Pooled logistic regression with irregular intervals

text: I'm planning an analysis that will require pooled logistic regression to analyze a longitudinal dataset with time-dependent explanatory variables (see here - D'Agostino, et al. 1990). However, the standard analysis uses regularly spaced intervals. My study has followup points at 0, 2, 6, 12, 24 months. What is a good way to model the effects of time given this structure of unevenly spaced intervals? Would a restricted cubic splines for the effects of time address the irregularity issue? Thanks!


title: Compare Ordered Lists that Include Frequncies

text: We have a lot data generated from a genetic algorithm. We run this algorithm multiple times with different 'initial populations', meaning different items chosen for the initial group of items. We wish to show that the results from different runs are the same (statistically very similar) independent of the initial population. The results are an ordered list of items from most frequently chosen to least frequently chosen. The data also includes the number of times that each item was chosen. We have many such runs and therefore many such lists to compare. (Example: Item_1 499, Item_2 327, Item_3 101, Item_4 99) We considered performing a spearman rank correlation on these lists, but we would like to preserve the frequencies. Many of the frequencies after the first few items are similar, making it seem that the ordering for results with similar frequencies in unlikely to be exactly the same. Can someone recommend a different method, which we may not be familiar with, which would be a better way to analyze this data?


title: Results of bootstrap reliable (hyperbolic distribution)?

text: I already asked a question about the interpretation of the results of the bootstrap algorithm in case of a normal mixture. This time, I fit a hyperbolic distribution. According to this thread: > the use of standard errors in this family is completely unreliable since the > distributions of the estimators are quite asymmetric as shown by the > histogram of the bootstrap samples To keep things short I will just give you the histogram plots. I am also using the alpha,beta,delta,mu parameterization. **Alpha** ![alpha](http://i.stack.imgur.com/EPbDR.png) **Beta** ![beta](http://i.stack.imgur.com/1wVEl.png) **Delta** ![delta](http://i.stack.imgur.com/Qj96P.png) **Mu** ![mu](http://i.stack.imgur.com/enwmL.png) From these samples of estimators I calculate the standard deviation and these give the standard errors of the alpha, beta, delta and mu estimates of the `hyperbFit` command. Are the distributions of the estimators asymmetric? Are my result reliable? I fit the hyperbolic distribution to financial return data.


title: Calculate pseudo-$R^2$ from R's zero-inflated negative binomial regression

text: I'm looking into calculating a Pseudo $R^2$ used McFadden's method for a zero- inflated negative binomial regression. I'm unclear how to go about evaluating $\hat L(M_{intercept})$ in R. Any suggestions for how this might be easily done? R Code thus far:               > require(pscl)     > require(MASS)     > ZerInflNegBinRegress<-zeroinfl(y~.|x+z, data=DATASET, dist="negbin", EM=TRUE)      Which returns the Log Likelihood for the model using the summary function. It's finding the Log Likelihood for an intercept only function that I'm unsure of.


title: Basic epidemiological models in R

text: To support an education project around neglected infectious diseases in Sub- Saharan Africa, I aim to estimate and report some very basic epidemiological models in R to build intuition and understanding with the students. I have been searching for a little bit (to no avail) - does anyone know of basic R-based modules / models (compartmental or agent-based) that I could use? For example - SIR models similar to these Python models (http://wiki.deductivethinking.com/wiki/Program_3.1)


title: Confused about relevance of Bayes' theorem to my problem

text: I've spent the morning teaching myself about Bayes's theorem, because I assumed that it was required to help me solve my problem. However, the answer I've ended up with is the same as I would have got if I had used my naive approach prior to understanding Bayes's rule. I'll explain my problem: Suppose there are four possible methods of purchasing some item, $A_1, A_2, A_3, A_4$. These methods have different levels of attractiveness to customers. Over the previous month, for example, 14,456 people have used one of the methods. Thus, knowing how many people used each of the methods, one can clearly see the probabilities that customer has used one of the methods:   * $A_1 = 1500 / 14456 = 0.103$   * $A_2 = 11000 / 14456 = 0.76$   * $A_3 = 56 / 14456 = 0.0038$   * $A_4 = 1900 / 14456 = 0.131$ Customers are pretty fickle though, and most of the time they don't go through with the purchase. In fact, the number of customers who actually _convert_ are pretty low. It's dependent also on the different method that they are using to make the purchase (imagine some methods are easier than others). In Bayesian language, I'm talking about the probability of a conversion given the method as $P(B | A_k)$. Using our example past month, we also know these probabilities, because we can see the number of conversions for each of the methods. There were varying conversion numbers for customers using each method.   * $P(B|A_1) = (76/1500) * A_1 = 0.0051$   * $P(B|A_2) = (24/11000) * A_2 = 0.0015$   * $P(B|A_3) = (10/56) * A_3 = 0.00065$   * $P(B|A_4) = (34/1900) * A_4 = 0.0022$ Now, let's suppose that sometimes the method we use to track the methods of purchase, and the number of conversions therein is a bit clunky and doesn't always work. Suppose a 145th conversion is in the system, but we don't know which method was used to make the purchase; we want to estimate this using the available data. Bayes' theorem tells me I need to calculate:   * $P(A_1 | B) = 0.5277$   * $P(A_2 | B) = 0.1666$   * $P(A_3 | B) = 0.0694$   * $P(A_4 | B) = 0.2361$ Great! So I now know that there was just over a 52% chance that the 145th customer purchased using method 1. However, this involved quite a bit of maths in making the $P(A_k | B)$ calculations. When I looked at the numbers that I got out I realised that if I had just divided the number of conversions for each channel by the total number of conversions, I got exactly the same answers.   * `76 / 144 = 0.527`   * `24 / 144 = 0.166`   * `10 / 144 = 0.069`   * `34 / 144 = 0.236` I guess I'm a bit surprised that I didn't need to go to all the effort with the Bayes's Rule stuff, so I don't quite understand situations where I'll need to use the rule, and not just do other summary calculations.


title: Definition of autocorrelation time (for effective sample size)

text: I've found two definitions in the literature for the autocorrelation time of a weakly stationary time series: $$ \tau_a = 1+2\sum_{k=1}^\infty \rho_k \quad \text{versus} \quad \tau_b = 1+2\sum_{k=1}^\infty \left|\rho_k\right| $$ where $\rho_k = \frac{\text{Cov}[X_t,X_{t+h}]}{\text{Var}[X_t]}$ is the autocorrelation at lag $k$. One application of the autocorrelation time is to find the "effective sample size": if you have $n$ observations of a time series, and you know its autocorrelation time $\tau$, then you can pretend that you have $$ n_\text{eff} = \frac{n}{\tau} $$ independent samples instead of $n$ correlated ones for the purposes of finding the mean. Estimating $\tau$ from data is non-trivial, but there are a few ways of doing it (see Thompson 2010). The definition without absolute values, $\tau_a$, seems more common in the literature; but it admits the possibility of $\tau_a<1$. Using R and the "coda" package:               require(coda)     ts.uncorr <- arima.sim(model=list(),n=10000)         # white noise      ts.corr <- arima.sim(model=list(ar=-0.5),n=10000)    # AR(1)     effectiveSize(ts.uncorr)                             # Sanity check         # result should be close to 10000     effectiveSize(ts.corr)         # result is in the neighborhood of 30000... ???      The "effectiveSize" function in "coda" uses a definition of the autocorrelation time equivalent to $\tau_a$, above. There are some other R packages out there that compute effective sample size or autocorrelation time, and all the ones I've tried give results consistent with this: that an AR(1) process with a negative AR coefficient has _more_ effective samples than the correlated time series. This seems strange. Obviously, this can never happen in the $\tau_b$ definition of autocorrelation time. What is the correct definition of autocorrelation time? Is there something wrong with my understanding of effective sample sizes? The $n_\text{eff} > n$ result shown above seems like it must be wrong... what's going on?


title: Reporting bootstrapped output for an independent samples t test

text: When reporting results of a bootstrapped independent samples t-test, can I just report the t statistic alongside the bootstrapped p value? Or do I also need to report confidence intervals?


title: Confidence interval for the biggest difference between means

text: I have a number of samples $X_1,\ldots,X_n $ from Gaussian distributions centered at $\mu_1,\ldots,\mu_n $ respectively. I need to construct a confidence interval for the biggest difference between means $\max_{1\leq i<j\leq n}|\mu_i-\mu_j|.$ How can this be done? Probably that interval should be wider than the ordinary t-test confidence interval for the difference between two means (due to selection bias), but tighter than, for example, confidence interval based on Tukey's range test (since only two of all the means are needed to be compared). This seems like a textbook question, would you please suggest where should I look for the answer? Thanks! Update. The task which I think might be solved with this confidence interval is the following. In a framework of equivalence testing the ordinary hypotheses $H_0\colon\mu_1=\mu_2 $ and the alternative $H_1\colon\mu_1\neq\mu_2 $ are usually swapped, resulting in null like $H_0\colon|\mu_1-\mu_2|\geq \varepsilon $ against the alternative $H_1\colon|\mu_1-\mu_2|<\varepsilon $ (see, for example, "Testing Statistical Hypotheses of Equivalence"). The problem with this approach is that $\varepsilon $ must be specified beforehand, while sometimes it would be more convenient to report something like the smallest $\varepsilon $ at which the null hypotheses of nonequivalence is rejected. For two-sample problem confidence interval for the difference $|\mu_1-\mu_2| $ provides a reliable alternative to equivalence testing: as far as I understood, $\varepsilon $ might be chosen as the smallest meaningful $\varepsilon $ such as $[-\varepsilon, \varepsilon] $ contains that interval. For multiple samples the situation is not clear to me. Sometimes it is suggested to construct CIs for all pairwise differences, but if we want all of them to be simultaneously valid, they must be enlarged with, say for simplicity, Bonferroni correction. But we actually don't need all of CIs to hold, because to claim the equivalence of all the means only $\max_{1\leq i<j\leq n}|\mu_i-\mu_j|<\varepsilon $ is needed. And, of course, the assumption of normality is not really desired, so it would be perfect to get rid of it in the end, but it might be kept for a start.


title: encrypted data on CART, ID3

text: Some data are confidential such as patient data. Therefore sometimes companies does not want to give original patient data instead they first encrypt it(for instance with SHA1) and then give. If we are given some encrypted data and run the following algorithms;   * ID3   * CART to generate decision trees, and then give those trees back to the company. After that company can decrypt it. However my question is; `Will the decrypted decision trees still be meaningful to company?`


title: What will happen if we sample the most probable value in the Gibbs sampling?

text: I am now working with the Gibbs sampling. One problem that puzzled me is that when we use the Gibbs sampling, we always sample randomly from the conditional probability. What will happen if we sample the most probable value instead? Thank you in advance.


title: Transform a non-monotonic value before decision tree (concrete example)?

text: Newbie question here. I am building a toy decision tree to differentiate personal names from business, government, or organizational names, like: AAA ENTERPRISES LLC DBA AAA BBB SERVICE   SMITH BARBARA   EXAMPLE FLOOR COVERINGS   BROWN JOSEPH A   2013 HOLDINGS L L C I'm mainly using continuous values, such as the # of tokens, the length of the string, and the % of tokens which are "dictionary words." However, I've noticed by eyeballing the data that if there are one or two trailing singleton chars (as in "BROWN JOSEPH A" or "HUTTON E F"), it's often a name, whereas if there are more than two, it's often a business (like "2013 HOLDINGS L L C"). The value of this signal is non monotonic and probably goes like this:               Num Singletons:    Relative size of effect:     0                  0     1                  +Large     2                  +Smaller but positive     3                  -Moderate     4                  -Large      I recall reading about "recentering" although this isn't quite exactly what I need, I fear. But I'm fairly confident there is a well-known technique for handling these non-monotonic suspected correlations early in the processing / transformation pipeline. Can someone point me in the right direction? Heavy lifting not required, just looking for a Wiki or CV link, or a couple terms of art I can Google.


title: ML Bank transactions assignement to invoices

text: In a effort to reduce human intervention, I'm trying to optimize the process of assigning bank transactions to invoices. This task should be done once every year, so we can assume our dataset won't take more much more than 200 entities. About our invoices, they are already categorized, for example as "Water expenses","Electricity", etc... Each invoice has a value (€), a date and a description. About the information the bank transactions provide, we have a value (€), a date and a description too. We can't naively compare each attribute in both dataset's because there might be outliers. What I mean is:   * Difference between the invoice date, and transaction date are going to happen.   * A transaction might not have the exact value of an invoice, it's possible a transaction was made to pay for different invoices and there's also the possibility an invoice required multiple payments at different times.   * About the bank descriptive, it has no similarity with the invoice description, and there might be no pattern. My idea, preprocess the data through clustering. Apply K-means on the bank transaction feed, has there will be similarity between groups (water expenses hopefully won't vary much in value, descriptives will be more or less similar), and after clustering try to assign each cluster to each category. The problem is: the features all have different types, the description similarity we can compute through jaro-winkler's distance, the price tag: a simple euclidean distance will suffice, and the date: we aren't really looking for similarity, instead what I want is periodicity between dates (any metric I can use here btw?), how can i apply k-means with different feature types?! This is my first real world problem with machine learning, so I'm looking for advice, is there any better way to solve the problem? Are there any flaw's in my plan? Would you approach the problem in a different way?! Thanks


title: Testing change in categories?

text: I hope that someone will be willing to help me with the following question. I have data on the numbers of hours participants spent in 6 programs over 2 years. I have already tested whether the number of hours increased over time, but I also want to know whether the proportion of time spent in each of the programs changed over time. Eg. Program1/Time1 = 2920 hrs/10000hrs; Program1/Time 2 =3200/12000; Program2/Time1 = 4200/10000; Program2/Time2 = 3800/12000, etc. Would a 2 sample z test for proportion change on each program (perhaps with Bonferoni correction) be appropriate? And/or is there another test(s) That I should be using? Any help would be much appreciated! I am comfortable with SPSS and am starting to learn R, if that makes a difference...


title: Multi-class classification via all pairwise classifications with LDA

text: I have trained linear discriminant analysis (LDA) classifiers for three classes of the IRIS data and struggling with how to make the classification. Here is the procedure: For the Iris data, I have 3 combinations i.e. (0,1), (0,2) and (1,2). So, I trained a simple binary LDA classifier for each combination, and ended up with three classifiers:               Classifier(0,1)     Classifier(0,2)     Classifier(1,2)      Now, say I need to classify an input, say k = [1.2, 2.3, 5.0]. What I am doing is passing this input through all the classifiers individually, which are giving me their respective scores, like:               Classifier(0,1)[k] = {0: some score, 1: some score}     Classifier(0,2)[k] = {0: some score, 2: some score}     Classifier(1,2)[k] = {1: some score, 2: some score}      In a simple binary case of two classes, what we are taught to do is to take the class with maximum score as the result. My question is, what to do in such a scenario, where I have three results from three different classifiers, and I want to classify the output. Please note that I am not using a multiclass LDA. I am just using a binary LDA for all the possible combinations, a technique which is stated here: http://en.wikipedia.org/wiki/Linear_discriminant_analysis#Multiclass_LDA > quoting the last paragraph of this section: " _Another common method is > pairwise classification, where a new classifier is created for each pair of > classes (giving C(C − 1)/2 classifiers in total), with the individual > classifiers combined to produce a final classification._ " Can somebody please enlighten me about what needs to be done in such a case for classification? Thank you.


title: How to do the prediction of user action?

text: I have a dataset, which is based of a log.The data is as following:                  bid LOGIN UPDATE_BILLING CREDITCARD_APPROVED CREDITCARD_DECLINED UPDATE_MONTHLY_BUDGET CREATE_AD UPDATE_AD_BID UPDATE_AD_BUDGET UPDATE_AD PAUSE_AD     1   -1     0              0                   0                   0                     0         0             1                2         0        0     2    1     0              0                   0                   0                     0         0             1                2         0        0     3   -1     0              0                   0                   0                     0         0             2                3         0        0     4    1     0              0                   0                   0                     0         0             3                3         0        0      bid is the label. 1 means increase bid, -1 means decrease. the other column is the user actions. For each time a user change the bid, I count how many different actions this user did during the past one day. If the action happens once, it will be 1. AS you can see, most of the table is 0,because this action does not happen. I was asked to predict the bid will increase or decrease. I feel it is not predictable. I train a simple naiveBayes model to see what happens. the accurency is only 0.48. Do you think the increase/decrease is predictable? do you have any suggestion about how to do the prediction?


title: Best method to predict binary outcome with multiple records per subject

text: I am interested in building a model to predict the binary outcome, retention (1 - retained; 0 - not retained) with various potential predictor variables (either continuous or categorical). With that being said, I have a dataset containing multiple records (magazine subscriptions) for some subjects. For example, I have 4 records at the magazine subscription level for Joe Smith (3 of which he is retained, while the other record he is not retained) and 7 magazine subscription records for John Doe (4 of which he is retained, while the other 3 records he is not retained). My initial thoughts were to use logistic regression using a single randomly selected magazine subscription for each subject. For example, randomly selecting 1 of the 4 magazine subscriptions for Joe Smith and 1 of the 7 magazine subscriptions for John Doe. Obviously, I would be losing a great deal of my data, which makes me think that there has to be a better way of modeling this data. What method would you be best to predict retention with data such as this? Thank you!


title: Pearson's coefficient and p-value. How do I work out proportions from this?

text: My results of a Pearson's Correlation are shown below:               Pearson’s Correlation – Test for Association       Pearson's product-moment correlation:     Test statistic: r = -0.1146235     Degrees of freedom = 363     Statistical significance: P = 0.02855      Reject null hypothesis, there is a small negative correlation between TV hours and the number of books read. Based on the negative $r^2$ value. The p-value is ≤ 0.05. What proportion of the variation in number of Books read has been explained by the variation in hours of TV watched? My answer is r2= -0.01313854675 x 100% = 1.3% in favour of TV watching. Is this correct?


title: The distribution of STD/MAD for a Student

text: Where X ~ a symmetric Student T Distribution $t_\alpha$, with power law tail $\alpha>2$, looking for the distribution of $\frac{(\sum_{i=1}^n x_i ^2 )^{½}}{\sum_{i=1}^n |x_i|}$, in order to estimate the properties of the statistic STD/MAD for a sample size $n$. The problem is that both measures, MAD and STD will be correlated for a given draw, so we cannot work with two different distributions.


title: Weighted least squares for energy data

text: I am trying to do a weighted multiple least squares regression on utility data. Basically I have utility bills where I use : Days billed, Consumption. For the same billing period I calculate Heating Degree Days and Cooling Degree Days. $ HDD = (18-Tout)^{+} $ and $ CDD = (Tout - 21)^{+}$ Where the + sign here means it equals zero if negative. Tout = Outside Temperature, given in Celsius. So, I used to do it like this : $$ Y = X\beta $$ Where : Y = Column vector, size (12,1), with each "monthly" consumption X = Matrix of size (12,3). Column 1 is the number of Days in the billing period (i). Column 2 is HDD(i) and Column 3 is CDD(i) $\beta $= The coefficients I'm looking for. Size (3,1) $$ \beta = (X^{T}X)^{-1}X^{T}Y $$ That works almost perfectly in most cases. But the fact is that each billing period is NOT equal in length, even if billed monthly (take February...). And in most cases the reading date varies. So that method will give me a non-zero Net Bias Error due to data length variation. That's where I'm struggling. I tried doing the following:   * redefining X as having the first column being full of 1.   * Defining a matrix W that is diagonal. Each $ Wii = Days(i)/DaysTot $ And solving for $$ \beta = (X^{T}WX)^{-1}X^{T}WY $$ But I cannot seem to make it work. What am I missing? Can you point me to the right direction please? Sample Data below, each line corresponding to one billing period (almost a month, but not exactly and that's the problem). You have Y (the kWh reading), the days in that billing period, and the corresponding sum of daily HDD and CDD.               Y      Days  HDD    CDD     30,343  31  0.00    0.00     25,193  31  0.00    29.67     24,719  30  0.00    70.67     66,993  31  41.61   0.00     202,271 30  166.61  0.00     297,333 31  243.44  0.00     393,891 28  341.33  0.00     505,941 31  433.83  0.00     485,715 31  431.06  0.00     412,160 30  350.67  0.00     346,495 30  298.44  0.00     138,180 36  86.22   0.00      Here is a spreasheet where I explain the problem more, show one example of what I'm trying to do in a case where there is only heating involved, and then one case where there is heating and cooling. For heating and cooling, I show what I've been doing so far, and then three attempts at correcting my methodology. The last attempt could be the solution but I am not sure at all. Download the Excel file here (.xlsx, 56ko): http://www.filedropper.com/weightedregressionforenergy ASHRAE Guideline 14 describes the problem and its solution for a case where only heating is involved (one variable): http://gaia.lbl.gov/people/ryin/public/Ashrae_guideline14-2002_Measurement%20of%20Energy%20and%20Demand%20Saving%20.pdf Read pages 139 (141 in the PDF) and the couple next ones. I am just trying to find a solution for a case where heating AND cooling are involved (two or more variables)


title: Re-ranking a list

text: I have a ranked list say A,B,C,D,E,F,G and a subset of this list ranked with another strategy lets say D,B,C,F. And my assumption is that the strategy used for ranking the sublist is better because its a feedback strategy. Now I want to re-rank my initial ranked list according to feedback ranking of sublist and initial ranking available. Please, suggest some good algorithm in this direction


title: Question about Tukey post-hoc ANOVA test results

text: I have used a Tukey post-hoc test after seeing that the ANOVA test indicated some significant differences between my groups (A,B,C,D). I am somewhat confused about the results. ![Plot and results](http://i.stack.imgur.com/cERWq.png) Why, for example, is A significantly different from B (although they are very similar); yet A is NOT significantly different from C and D (even though they are clearly VERY dissimilar)? Perhaps I do not fully understand the Tukey post-hoc test? If someone can direct me to a decent text on statistical analysis (graduate level), that would be greatly appreciated. Any comments are greatly appreciated.


title: Unexpected residuals plot of mixed linear model using lmer (lme4 package) in R

text: I have conducted an experiment with multiple (categorical) conditions per subject, and multiple subject measurements. My data-frame in short: A subject has one property, `is_frisian` which is either 0 or 1 depending on the subject. And it is tested for two conditions, `person` and `condition`. The measurement variable is `error`, which is either 0 or 1. My mixed linear model in R is:               > model <- lmer(error~is_frisian*condition*person+(1|subject_id), data=output)      However, the residuals plot of this model gives an unexpected (?) result. ![Residuals lmer model](http://i.stack.imgur.com/nz2KY.png) I was taught that this plot should show randomly scattered points, and they should be normal distributed. When plotting the density of the fitted and the residuals, it shows a reasonable normal distribution. The lines you can see in the graph, however, how is this to be explained? And is this okay? The only thing I could come up with is that the graph has two lines due to the categorical variables. The output variable `error` is either 0 or 1. But I do not have that much knowledge of the underlying system to confirm this. And then again, the lines also seem to have a low negative slope, is this then perhaps a problem? **UPDATE:**               > model <- glmer(error~is_frisian*condition*person + (1|subject_id), data=output, family='binomial')     > binnedplot(fitted(model),resid(model))      Gives the following result: ![binned residual plot](http://i.stack.imgur.com/XMXFx.png) **FINAL EDIT:** The density-plots have been omitted, they have nothing to do with satisfaction of assumptions in this case. For a list of assumptions on logistic regression (when using family=binomial), see here at statisticssolutions.com


title: Evaluating a Survival Analysis model

text: I have just begun studying survival analysis and I am using it for an industrial problem where I am trying to build survival curves for customer liftimes (from when they sign-up to when they cancel). I am using the Aalen additive model using explanatory variables that I think affect the customer lifetime. What I am having trouble understanding is the following. What I am really interested in are the survival curves of active customers (i.e. censored). This way I can get some insight into how likely they are to cancel given how long they have been a customer. Coming from a ML background, I find it troubling to train the Aalen model with data points that I am interested in testing (i.e. overfitting). So would it be a better idea to only train the model on cancelled accounts, then test it on the active ones? I guess a generalized version of my question would be how do you know if your SA model is overfitting? This particular model uses regression, so wouldn't overfitting be a possible issue? Perhaps I am completely missing the point of SA and maybe it is not meant for prediction in this way


title: Which statistical analyses should I use?

text: I want to compare data from a treatment group vs a control group over time (baseline, 1 yr, 2 yr). The measurements are repeated for all subjects across all time points (i.e., all participants complete the same set of questionnaires at all time points). Some participants have missing data (some did not complete questionnaires at certain time points). Control group has more subjects than treatment group (unequal sample sizes). I am not sure whether I should use mixed models analyses to handle missing data because I am not sure whether mixed models can be used for unequal sample sizes. Or should I use Scheffe's method (non-parametric variant of anova) to handle the unequal sample sizes. * * * I am interested to investigate whether scores on a questionnaire (that measures severity of symptoms) differ between the treatment group (given medication) and the control group (not on medication) across the different time points. I was thinking of using repeated measures mixed model analysis (specified as 'Linear Mixed Models analysis' in SPSS which involves selecting the appropriate covariance structure [e.g.compound symmetry and 1st order auto regressive]) because all subjects undergo the same questionnaire across all time points and also due to incomplete data from some participants. The only issue is that I have 300+ participants in the control group and 30+ participants in the treatment group. (very unequal group sizes). In this case, will the repeated measures mixed model analysis still be appropriate? I agree with Peter that performing multiple t-tests will not be advisable in this case as it also increases family-wise error (hence increases tendency of committing a type 1 error). Please kindly clarify. Thanks!


title: Am I using Fisher method correctly?

text: Suppose I want to test the hypothesis that speed of light is always below $c$ (note: regardless of the medium). [Also, this is just a fake example, so pretend we don't know any physics theory] Now, I did measurements of speed of light in 30 different mediums. Suppose the readings come with error (Gaussian with known s.d.). Suppose one of my readings is $1.2c \pm 0.1c$ - I should be somewhat worried for my hypothesis, since the p-value is about 0.025 (half of $2\sigma$ distance). But say my other 29 readings are something like $0.6c \pm 0.0001c$, $0.47c \pm 0.0001c$,... (basically confidently below $c$). Then the p-value for each is something like 1. Then if I use the Fisher method, I obtain a total p-value of 1. It seems that all worries are gone! But in reality I dont really care about the "confidently-below $c$" points. Is it wrong of me to be combining all the p-values? After all, my hypothesis is about all mediums... I'm supposed to devise the statistic before performing the experiments, so I am not supposed to "cheat" and only pick the medium for which $1.2c \pm 0.1c$. I mean, if my data instead has $1.2c \pm 0.1c$, $1.3c \pm 0.2c$, $0.9c \pm 0.3c$,... I won't know how to "cheat" so I think we have to combine results from all the mediums?


title: How to analyze data which might come from a few normal distribution concatenate together in order?

text: For example, I have a series of values for example like the following:               data <- c(rnorm(10000,40,1500),rnorm(9000,-35,1400),rnorm(11000,30,1300))      I don't know the mean, don't know the sd of my data, and am particularly interested in the "sign of the mean" of my data (I can do a t-test but most often it is not significant, and even if it is, it only tells half of the story), and would be better to know if it changes over time. (how ever, I believe the data is not a time series, i.e. consecutive data points are independent of each other). But, I also doubt that my data might have different baseline distributions (assumed normal) at different time stages. So, what is a standard way of analysing such data? You can see the sd of my data is quite big. Do I have to pretend it is a time series and decompose it to see the trend?               df <- data.frame(x=1:30000, data=data)     ggplot(df, aes(x=x)) + geom_bar(aes(y=cumsum(data)), width=1, stat="identity")      This is the graph of cumulative values of the data. Can we tell if it is from one normal distribution, or several normal distributions with different means and variances? If it is possible, what assumptions do I have to make? ![enter image description here](http://i.stack.imgur.com/SOsHq.png) Compared to this:               data <- rnorm(30000, 3, 1800)     df <- data.frame(x=1:30000, data=data)     ggplot(df, aes(x=x)) + geom_bar(aes(y=cumsum(data)), width=1, stat="identity")      ![enter image description here](http://i.stack.imgur.com/KLnYR.png) # Update: I guess that `changepoint` package in R is good enough               > a <- cpt.meanvar(data=c(rnorm(10000,40,1500),rnorm(9000,-35,1400),rnorm(11000,30,1300)), method="BinSeg")     > summary(a)     Changepoint type      : Change in mean and variance      Method of analysis    : BinSeg      Test Statistic  : Normal      Type of penalty       : SIC with value, 20.61791      Maximum no. of cpts   : 5      Changepoint Locations : 10264 19013           > b <- cpt.meanvar(data=rnorm(30000, 3, 1800), method="BinSeg")     > summary(b)     Changepoint type      : Change in mean and variance      Method of analysis    : BinSeg      Test Statistic  : Normal      Type of penalty       : SIC with value, 20.61791      Maximum no. of cpts   : 5      Changepoint Locations :        PELT would give some small intervals thought I had penalty to include the number of change points.               > a <- cpt.meanvar(penalty="SIC1",data=c(rnorm(10000,40,1500),rnorm(9000,-35,1400),rnorm(11000,30,1300)), method="PELT")     > summary(a)     Changepoint type      : Change in mean and variance      Method of analysis    : PELT      Test Statistic  : Normal      Type of penalty       : SIC1 with value, 30.92686      Maximum no. of cpts   : Inf      Changepoint Locations : 10049 18984 18987      >      > a <- cpt.meanvar(penalty="SIC1",data=c(rnorm(10000,40,1500),rnorm(9000,-35,1400),rnorm(11000,30,1300)), method="PELT")     > summary(a)     Changepoint type      : Change in mean and variance      Method of analysis    : PELT      Test Statistic  : Normal      Type of penalty       : SIC1 with value, 30.92686      Maximum no. of cpts   : Inf      Changepoint Locations : 8810 8814 18962


title: Problem with big Data in TraMineR

text: I would like to analyse the trajectories of patients hospitalized for CVA (cerebrovascular accident) with the package TraMineR. I have a large dataset with around 25,000 trajectories observed during 365 days. There are 5 states in my data. When I try to use the seqdist function which computes the distances between sequences, I have problems with computer's memory. Some people suggested to aggregate identical sequences, compute the distance and cluster the sequences using weights. But this solution is not sufficient since I have 20000 unique sequences in my data. An other solution proposed in the web is to analyze a sample of the data. So, I have analyzed a sample of 5000 trajectories in my data. Is there another solution by using all the trajectories ? Thank you in advance for your answer. Christelle FERREIRA


Can unimodal prior and unimodal sampling distributions lead to a multimodal posterior distribution? The Bayes rule tells us that $$ f(y|x) = f(x|y) f(y) / f(x). $$ which, I think, implies a unimodal prior and unimodal sampling distributions can only lead to a unimodal posterior distribution. If $y$ is multidimensional, can the marginal posterior distribution of some component of $y$ be multimodal? (Related question: can a unimodal multivariate distribution has some multimodal marginal distribution? If all marginal distributions are unimodal, can the multivariate distribution be multimodal?)


title: Cronbach's alpha for one scale but different types of items

text: I developed a questionnaire on art expertise. There are 3 rather subjective questions("I know much about art","I know more than the average citizen","I am able to identify an artwork by a popular artist") evaluated on 7-point likert scales (not at all - absolutely). Furthermore, there are multiple rather objective questions like ("How many hours do you deal with art publications per month") etc. with differing numbers of answering options, scaled ordinally and coded with some number from 1-n. Now I want to use Cronbach's alpha to assess the internal consistency of my questionnaire. I used it on the subjective scales and got a alpha=0.93 which is not too bad. Now I wonder: can I use it on the objective scales as well, although they have different numbers of options, or can I even use it on the whole questionnaire? Answer option alternatives range from 2-6. Thanks and best, Christian


title: Example of the impact of statistical inference in society

text: When giving a talk or explaining my latest results, several colleagues of mine downsize the importance of the statistical analysis part as being "only statistics". In my opinion, the impact of statistical inference is often underestimated, even among so-called educated people. Could you provide a clear example of an important and high impact societal/political decision we based on statistical inference?


title: Statistically prove the effectiveness of a treatment using GLM repeated measurements

text: I have two lots of samples: one is the control lot and the other undergoes some treatment. I did three measurements for the samples: one at the initial time (`T1`) and the other two later.                          Descriptive Statistics                        Lot         Mean  Std. Deviation      N     T1  1.00    124.3043         3.21127     23         2.00    124.1333         1.94286     30         Total   124.2075         2.54467     53     T2  1.00    112.8261         5.81262     23         2.00    107.7000         7.42387     30         Total   109.9245         7.18398     53     T3  1.00    90.3478          8.47783     23         2.00    114.7000         4.43458     30         Total   104.1321        13.77852     53      The plotted data look like this: ![enter image description here](http://i.stack.imgur.com/5NjiA.jpg) **What I did:** Use IBM SPSS (ver 20) and run `General Linear Model/Repeated Measures`. Established the factor to `MeasureTime` with three levels and the measurement. Between subjects factor was set to the lot number (`1` or `2`), `Model full factorial`, `Contrasts` to `MeasureTime`, `Simple`, `Reference category` set to first. SPSS output a lot of statistics. Almost all are relevant (sig < 0.000) except Mauchly's Test of Sphericity. **What I want to know:**   1. How can I see if the mean differences in `T1`/`T2`/`T3` are significant? (I expect that `T1` is not and `T3` is.)   2. What statistic (from SPSS output) tells me that the treatment works?  P.S. I read other related questions suggested by Stack Exchange, but no luck with my questions. I think that ANOVA would work for a factor with multiple levels, but as I read on other questions, in case of repeated measurements, ANOVA assumptions fail. The data are not independent.


title: Regression specification

text: Hello statistic friends, do I make any econometric mistakes if I regress: $D/A = a + b_1 * A + b_2 * X/A + b_3 * Z + U$ D = Debt, A= Assets, X = Cash Flows, Z = Controls, U = error terms I am a little bit confused since I have A in the denominator for my dependent variable and for one independent variable (but not in my controls). Simultaneously I use A as another independent variable. Any help is highly appreciated. Thanks!


title: Best introduction to longitudinal data analysis

text: I have an undergraduate in statistics and I like some mathematical explanation and exposition. I would like to learn the basics of longitudinal analysis for my leisure and edification. What, if any, is the best mezzanine treatment of this subject? EDIT: to give gung's clarifications.


I do not know any multimodal distribution. Why all known distributions are unimodal? Is there any "famous" distribution that has more than one mode? Of course, mixture of distributions is often multimodal, but I would like to know whether there exist any "non-mixture" distribution that has more than one mode.


title: Deduce variance from boxplot

text: I was wondering how to deduce the variance of a variable using a boxplot. Is it at least possible to deduce if two variables have the same variance observing their boxplot?


how can i use regression model selection technique and to carryout ANOVA for calculated data I have problem regarding to know , what kind of relationship is excised between parameters to produce a good response and how to eleminate the other factors. And a relationship between response and factors of it. so, will any body tell me how to do it, and difference between response surface methodology and regression analysis and even I wana know is there any other techniques for it. what is anovo and how it has to be done


title: Any Statistical Benefit of Using Composite Scores?

text: I know composites gives an overall assessment of multiple tests. But is there any statistical benefits to using composites? I have noticed in my composite scores, the composite turned out to be significant, but the individual components are not. I used the z scores to calculate composites of 3 variables. Why is this? Thank you so much


title: Creating a model for prices including supply

text: I'm working on modeling secondary market ticket prices for sporting events, but the issue I'm running into is that the model (a linear regression) assumes that more season ticket holders and more tickets on places like StubHub mean that the game is more desirable, and thus the tickets should be worth more. The problem is that one year has a significantly higher number of both due to increased team performance, so it thinks prices should be going up when the increased supply of tickets on the secondary market is actually driving prices down. Is there anything I can do to address this problem?


title: Linear Mixed Models with variable time points in SPSS

text: I am analyzing a clinical study measuring patients' symptoms and brain structure sizes over three time points. The "visits" should be at baseline, 3 months and 12 months but they vary considerably and I would prefer to run the measurement as "days from baseline". Just not sure how to set this up in SPSS. So I have set: Symptom score as a DV and "days from baseline" as a covariate. "days from baseline" as fixed effect. "days from baseline" as random effect with intercept, variance component. Subject grouping by ID. When I create a line plot of the predicted values vs days each ID has a reasonable trend but each has three growth curves that do not connect into one line. I am expecting this is wrong but unsure how to fix this. Advice?


title: Where does the bootstrap beat the "classical" estimates?

text: Although this question is kind of similar to Please provide an example of when bootstrap has less bias than classically approximated estimates? I would like to look at the topic from a more general point of view. As far as I understand it, the major advantages of the bootstrap are:   1. The bootstrap can deal with any distribution of the underlying sample, no matter if parametric or non-parametric. Maybe you'd have to do some smart trick, e.g. smoothing the edf but in general you can tackle any given sample.   2. The bootstrap can deal with statistics, that are rather uncanny to handle with "classical" approximation, e.g. constructing CI's for the median. But are there, aside from these two benefits, any real advantages in performance? Can anybody give me an example of a distribution and a statistic, where the bootstrap estimate is actually more accurate (has less bias/ less variance/ a better covering probability) than estimating by doing calculations based on CLT or equivalent methods that do not involve resampling? Our professor said, the accuracy of the bootstrap compared to the "classical" method depends on the skewness of the underlying distribution, i.e. the less skewness, the better the bootstrap. Unfortunately he didn't explain it in more detail.


title: Comparing correlations for significance

text: I am wondering how I would do the following in SPSS. I have a question (question 1) that was asked once before treatment, then the group was split into two different treatment groups (let's say `A` and `B`), then both groups were asked question 1 again. I would like to find out whether the answer to question 1 changed significantly from time 1, before treatment, to time 2, after treatment. Then I would like to see if there is a significant difference between groups `A` and `B` on the question 1 correlations mentioned in the previous sentence.


title: Understanding the Rank Probability Score

text: The ranked probability score (RPS) is a measure of how good forecasts that are expressed as probability distributions are in matching observed outcomes. Both the location and spread of the forecast distribution are taken into account in judging how close the distribution is to the observed value. \begin{equation} \mathrm{RPS}=\dfrac{1}{r-1}\sum\limits_{i=1}^{r}\left(\sum\limits_{j=1}^{i}p_j-\sum\limits_{j=1}^{i}e_j\right)^2, \end{equation} where $r$ is the number of outcomes, $p_j$ is the forecasted probability of outcome $j$ and $e_j$ is the actual probability of outcome $j$. The special case $r=2$ gives the brier score. The score lies between 0 and 1, with lower scores being better, but does a value of 0.8, say, imply your forecasts are not very good? Or is this similar to $R^2$, for example, where a low value does not necessarily imply you have a bad model and it is better suited for comparing between models. At first I thought it's unlikely that it would say anything about your forecasts (not relative to another model, at least), but then consider two outcomes ($r=2$) and say your model predicted each outcome to occur with probability 0.5. Then the RPS gives 0.25 regardless of the outcome. So generally, a value greater than 0.25 for the RPS says your model predicted the wrong outcome and a lower value implies your model predicted the correct outcome. Unfortunately, this idea does not work for $r=3$. Taking the probability of the three outcomes to be 1/3 gives RPS values of 5/9, 2/9 and 5/9 when outcome 1, 2 and 3 actually occurred, respectively. This does makes sense though, predicting outcome 1 when outcome 2 happened is better than predicting outcome 1 when outcome 3 occurred and this is exactly the motivation for the predictive performance evaluator. Nonetheless, I am back to thinking it is only suitable for comparing models. A second question: $R^2$ is the proportion of response variation "explained" by the regressors in the model. Is there a similar interpretation for the RPS? Finally, if a model makes a number of forecasts - say it predicts the next outcome, the outcome is observed, the model is updated and then the model predicts the next outcome - would simply averaging the individual ranked probability scores for each predicted outcome be appropriate, for comparing between two models? I think this is reasonable, but some formal reasoning would obviously be preferred.


title: Simulate missing data for multivariate distribution?

text: (Following the private request from a more senior CV member, I am editing this question to make it more readable and comply with CV standards). I am looking for a method to simulate multivariate, non-normal data with a pre-specified pattern of missingness while making sure that the missingness on does not invalidate my non-normality. I thought about doing this as a two-step process: (1) Generate multivariate, non-normal data following the method described in Vale & Maurelli 1983 (the multivariate extension of Fleishman 1978). The Vale & Maurelli process allows me to specify the population values of skewness/kurtosis in all 1-dimensional marginals. (2) Once a sample is taken from said distribution with given population values of skewness/kurtosis on the one-dimensional marginals , generate some missing data pattern on the sample data. I immediatley recognized that by inducing missingness, my non-normality could disappear (e.g. if I sample from a skewed distribution and then have missing data on the higher values, the skewness could be severely diminished or disappear all togehter). My attempt at a solution was to introduce a 3rd step where (3) Once the data were sampled and the missingness introduced, I would calculate the sample estimates of skewness/kurtosis and only retain those datasets for further analysis that fell within certain "acceptable" (<--- this would be set by me) range of skewness/kurtosis. A professor commented on the fact that my step (3) would be invalid because #1 the definitions of skewness/kurtosis we have only exist for complete datasets and #2 if I only kept certain datasets and threw away others I could no longer claim I was sampling randomly from my distribution and I could only do this under MCAR (Missing Compeletely at Random) but not under MAR (Missing At Random) and NMAR (Not Missing At Random). So I am stuck now and would very much appreciate any insights on anyone who could have some idea about how to implement said algorithm. I believe maybe going about with a different sampling scheme than the one I suggested?


title: Synthetic Model Object

text: I want to test a multinomial logit model that has been published by another author against one that I have developed. I have his coefficients, but neither the standard errors nor the model fit statistics. Is there a way to create an `mlogit` object without estimating it? Can I create a blank object and fill it with the coefficient estimates? Could this be used for prediction? Are there pitfalls to doing so?


title: How to build the best regression model using stepwise regression?

text: I have a data set that contains rental cost (Dependent variable) and I have another data sets that are potential predictors like city, area, No. of rooms etc. (Independent variables). I have performed a stepwise regression to find the best model. I stopped when all other independent variables showed high p-values and I ended up with three variables (all of them have below 0.05 p-value). When I tried to validate the model by entering only these three variables as independent variables, one of the variables had a high p-value different than it showed in the stepwise steps. In short, why even after I have performed the stepwise regression my model included an above 0.05 p-value variables?