title: Causal models in pymc3

text: I'm trying to fit a causal model. Participants in a task are trained on the model then asked for their belief in all the joints over the variables (e.g., what are the chances of observing an item with `cause_1` and `effect` present and `cause_2` absent?) I want to infer people's beliefs in the strengths of two causes, as well as various other beliefs (background causes of the effect, likelihood of the causes, etc.). This is a very simplistic model to get my toes wet. My specification is as follows:               # Common effect model with 2 causes (c1 and c2) causing an effect e.     # Causes combine using a 'noisy or' rule: two causes independently bring about the effect with probability m.     # r is the likelihood of the causes.     # bg is the probability of the effect coming about without a cause.          ncauses = 2     nvars = ncauses + 1          joints = np.array([[c1, c2, e] for c1 in [0, 1] for c2 in [0, 1] for e in [0, 1]])     resps = (np.array([10.0, 1.0, 3.0, 10.0, 4.0, 10.0, 4.0, 11.0])-1)/10 # converting from their 11-point scale.               def noisy_or(r, m, bg):         pc = switch(joints[:,:-1], r, 1-r)         pe_absent = (1-bg) * prod(1-r*m)         pe = switch(joints[:,-1], 1-pe_absent, pe_absent)         return prod(pc, axis=1) * pe          def beta_recoded(name, ratio, confidence, **kwargs):         beta = confidence / (ratio+1)         alpha = confidence - beta         return Beta(name, alpha=alpha, beta=beta, **kwargs)          with Model() as model:         r = Beta('r', alpha=1, beta=1, shape=(1, ncauses))         m = Beta('m', alpha=1, beta=1, shape=(1, ncauses))         bg = Beta('bg', alpha=1, beta=1)         pjoints = noisy_or(r, m, bg)         ratings = beta_recoded('ratings', pjoints, 5, observed=resps)      With this setup, find_MAP chokes. Did I do things reasonably? One of the weird things is that when calculating the probability of each of the joints, there's a sort of if-then structure (if var is absent, we want 1-p, else p) which I turned into abs(is_absent-p), I wasn't sure what was the most elegant way to do that. EDIT: this has been fixed using switch. EDIT: `pymc3` appears to require the development version of Theano in git.


title: Calculating the std dev of a 30 team league with each team having a 50% chance of winning

text: I just want to preface this with the fact it is not a home work question. I am also not a stats person so I am not sure even how to start with calculating this, what it is called or where to look. When I google I get not much that looks like it helps. So my problem, assume you have a 30 team league. Any game between two teams has a 50% chance of winning or losing, essentially a coin flip. I know that at the end of the runs all the teams wont have equal amount of wins, you'll get a normal distribution or bell curve. How do I calculate what the standard deviations of the win% for each team would be? Let's assume 500 games. How would I solve this? How would I even start? Thanks


title: How is a statistic represented as a mapping with varying sample size?

text: Suppose there is a sample of varying size $n\in \mathbb N$, each sample point taking values in $\mathbb R$, and a statistic $T$. If I am correct, a statistic can accept arbitrary sample size. What are some standard ways to represent/write the statistic $T$ with varying sample size $n$? For example,   1. We may write the statistic $T$ as a mapping from $\mathbb R^n$ to $\mathbb R$ for each different $n$ value. How can we make sure that, the different mapings for different $n$ values are from the same statistic $T$ not from another different one?   2. Can we write $T$ as a single mapping for all $n$ values? How is the statistic $T$ represented as a mapping with varying number of arguments (here arguments are sample points in the sample) then? Thanks and regards!


title: How to specify a path model with exogenous variables jointly affected by an unmeasured variable

text: I’ve simulated some data consisting of one response variable (‘y’) and two collinear predictor variables (‘Amount’ and ‘MPS’), where collinearity arises from one of two causes: (1) Amount causes MPS, or (2) Amount and MPS are jointly affected by an unmeasured variable. What I'm trying to do is figure out whether path analysis can discriminate between these two causes of collinearity. But I'm having trouble specifying a path model for collinearity scenario (2). My question: **Is it possible to specify a path model that implies that two exogenous variables are jointly influenced by an unmeasured variable?** I'm working in lavaan, but answers for how to do this conceptually would also be appreciated (if you aren't familiar with lavaan). Here are my data, simulated in R:               # collinearity cause (1)     Amount <- rnorm(n=350, mean=0, sd=1)        MPS <- rnorm(n=350, mean=0.76*Amount, sd=0.653)      y <- rnorm(n=350, mean=0.367*Amount + 0.367*MPS, sd=0.72)           # collinearity cause (2)     Lurking <- rnorm(n=350, mean=0, sd=1)      Amount <- rnorm(n=350, mean=0.872*Lurking, sd=0.486)       MPS <- rnorm(n=350, mean=0.872*Lurking, sd=0.486)       y <- rnorm(n=350, mean=0.367*Amount + 0.367*MPS, sd=0.72)               And this is my path model for (1), specified in lavaan:               model1 <- '       #regressions       y ~ Amount       y ~ MPS       MPS ~ Amount       '      And this is a path model I tried for (2):               model2<- '       #regressions       y ~ Amount       y ~ MPS       #residual correlations       MPS ~~ Amount       '      so for path model (2) my approach was to specify a residual correlation between MPS and Amount. I'm uncertain if this is the correct approach. but even if it is, it doesn't work – to make it work I have to specify that exogenous variables are not fixed, and this uses up my degrees of freedom so I can't test the model. If anyone has any suggestions for how I can do this – or if it is possible at all - I'd really appreciate it.


title: use Markov model to analyze communication systems?

text: I want the basic steps to analyze a communication system using Discrete Time Markov Chain? Can you help me with a book name that concern with this issue?


title: Sandwich covariance for robust regression using M estimators for data exhibiting heteroskedasticity

text: Following the answer and comments on Python Statsmodels Testing Coefficients from Robust Linear Model based on M-Estimators: I'm wondering how the sandwich form of an m estimator's covariance matrix would be calculated. The RLM regression will devalue outliers, but my understanding is that it does not control for heteroscedasticity. I have both systematic outliers and heteroscedasticity in my data and am interested in how to come up with a heteroscedastically controlled covariance matrix for the RLM. Additional information: using the Huber's t norm at tuning parameter = 1.345 and MAD scaling.


title: Optimal probability distribution with no knowledge about the base rates

text: If we have absolutely no knowledge about the base rates (maybe except for the number of possible outcomes), what is the most accurate probability distribution to assume, is it uniform? Here is a practical example. I don't know whether the Universe is real or I am just living in a simulation. In either case, my observations would be exactly the same. Furthermore, I haven't got the slightest clue (let's assume for the sake of the question that I don't) about the probability of the existence of a simulation, or the probability of the existence of a universe like ours, or the number of each. Also, a simulation does not necessarily have to occur in a universe with the same physical laws - as far as I know, in the underlying universe even the laws of logic might be so different that 1 + 1 equals 3. Note that this question is not about philosophy, it is about the probabilities.


title: How can I simulate runtime data containing workflows in R?

text: I am a PhD student, and I am trying to construct a data set to use for proof of concept for my research. However, I am not sure how to construct the model to simulate the data in R. I want to construct data that mimics software usage logs. In other words, each record would represent a user running an application that is part of a 50 application distribution, so the data fields I need are:   * Application name   * User name   * Session ID (one or more applications are run by a user in each session)   * Start time (time the application was launched)   * End time (time the application was closed)   * Runtime (End time - Start time) Within this set of runs, I would like to embed the following characteristics:   * Clusters of applications that are frequently used together (same session and same user)   * Frequent workflows (during same session, users frequently use application A followed by C then E)   * Applications that can be removed from the distribution dues to infrequent use   * Applications that should be considered for aggregation (applications D and F are almost always used together, so further testing should be done to see if they should be combined into a single application) Given these features, I can not seem to come up with any sort of model to fit the data to. I am at a loss for where to even start. Any recommendations you could give would be greatly appreciated with regards to an R script that would generate data such as this. Thanks in advance.


Does there exist a conjugate prior for the Laplace distribution? If not, is there a known closed form expression that approximates the posterior for the parameters of the Laplace distribution? I've googled around quite a lot with no success so my current guess is "no" on the questions above...


title: What are the stationarity requirements of using regression with ARIMA errors for inference?

text: What are the stationarity requirements of using regression with ARIMA errors (dynamic regression) for inference? Specifically, I have a non-stationary continuous outcome variable $y$, a non- stationary continuous predictor variable $x_a$ and a dummy variable treatment series $x_b$. I would like to know if the treatment was correlated with a change in the outcome variable that is more than two-standard errors away from zero change. I am unsure if I need to difference these series before performing the regression with ARIMA errors modelling. In an answer to another question, IrishStat states that `while the original series exhibit non-stationarity this does not necessarily imply that differencing is needed in a causal model.` He then goes on to add that `unwarranted usage [of differencing] can create statistical/econometric nonsense`. The SAS User Guide suggests that it is fine to fit regression models with ARIMA errors to non-stationary series without differencing so long as the residuals are non-stationary: > Note that the requirement of stationarity applies to the noise series. If > there are no input variables, the response series (after differencing and > minus the mean term) and the noise series are the same. However, if there > are inputs, the noise series is the residual after the effect of the inputs > is removed. > > There is no requirement that the input series be stationary. If the inputs > are nonstationary, the response series will be nonstationary, even though > the noise process might be stationary. > > When nonstationary input series are used, you can fit the input variables > first with no ARMA model for the errors and then consider the stationarity > of the residuals before identifying an ARMA model for the noise part. On the other hand, Rob Hyndman & George Athanasopoulos assert: > An important consideration in estimating a regression with ARMA errors is > that all variables in the model must first be stationary. So we first have > to check that yt and all the predictors $(x_{1,t},\dots,x_{k,t})$ appear to > be stationary. If we estimate the model while any of these are non- > stationary, the estimated coefficients can be incorrect. > > One exception to this is the case where non-stationary variables are co- > integrated. If there exists a linear combination between the non-stationary > $y_t$ and predictors that is stationary, then the estimated coefficients are > correct. Are these pieces of advice mutually exclusive? How is the applied analyst to proceed?


title: Cluster Analysis in R

text: I have a table of employees' information and each employee has the following attributes. I wanted to do an analysis to find out what similarities they share among themselves and possibly cut them into 3 or 4 groups. I started with cluster analysis using `fit = kmeans(mydata,5)` and `fit = dist(mydata,method="euclidean")` but got an error message `NA/NaN/Inf in foreign function call (arg 1) In addition: Warning message: In kmeans(sales1, 5) : NAs introduced by coercion` Do all my variables have to be numbers? Or I am using completely a wrong approach? List of the variables:               $ GENDER: int  2 2 2 2 1 2 1 2 2 1 ...     $ HISP_LATINO_DESC: Factor w/ 2 levels "Hispanic/Latino",..: 2 2 2 2 2 2 2 2 2 2 ...     $ APPRAISAL_CRITERION: Factor w/ 5 levels "","GD","GN","TN",..: 3 5 3 3 3 5 3 3 3 5 ...     $ P2_RT: Factor w/ 7 levels "","HH","HL","HM",..: 4 6 4 3 7 7 7 3 3 6 ...     $ GRADE: Factor w/ 8 levels "21","22","24",..: 6 4 5 2 8 4 6 5 2 4 ...     $ ACTION_RSN_DESC: Factor w/ 9 levels "IVT-Misconduct",..: 4 4 4 7 4 3 4 4 4 3 ...     $ Generation: Factor w/ 4 levels "Baby Boomer",..: 2 1 1 2 1 2 2 2 2 2 ...`     $ Date.of.Birth: Factor w/ 38 levels "1/4/1978","10/23/1968",..: 28 6 29 15 11 1 18 3 32 8 ...     $ Age: int  43 55 49 37 48 34 35 45 43 33 ...`     $ Date.of.Hire       : Factor w/ 33 levels "1/1/2012","1/2/1998",..: 13 32 7 10 21 16 6 5 15 23 ...     $ Term.Date          : Factor w/ 30 levels "1/30/2012","1/31/2012",..: 18 20 4 29 2 1 27 25 22 10 ...     $ Tenure             : num  13.4 3.9 5 5.9 4.7 0.8 1.9 7.6 7.4 1.8 ...     $ Tenure.Category    : Factor w/ 5 levels "1 to less than 3 years",..: 2 3 3 3 3 5 1 4 4 1 ...     $ Gender_Text        : Factor w/ 2 levels "F","M": 1 1 1 1 2 1 2 1 1 2 ...     $ Ethnicity          : Factor w/ 6 levels "Asian","Black or African American",..: 6 6 6 6 6 6 6 6 6 6 ...     $ Reg..vs..Non.Reg.  : Factor w/ 2 levels "Non-Regrettable",..: 2 1 2 2 2 2 2 2 2 1 ...     $ lateral            : Factor w/ 6 levels "0","1","2","3",..: 2 2 1 2 1 1 1 2 1 1 ...     $ promotion          : Factor w/ 3 levels "0","1","2": 2 1 3 1 2 1 1 2 2 1 ...     $ Geo                : Factor w/ 17 levels "Bentonville, AR",..: 1 1 1 1 1 1 2 2 2 3 ...     $ Avg..Time.for.Promo: num  9 0 2.4 0 2.4 0 0 5.5 2.3 0 ...


title: What are the 2nd derivatives of the log multivariate normal density?

text: I develop open-source statistical software (http://openmx.psyc.virginia.edu/), but matrix calculus is not my strong point. I need the 1st and 2nd derivatives of the log multivariate normal density. I was happy to find the 1st derivatives here on CrossValidated, How to take derivative of multivariate normal density? However, the 2nd derivatives are left as an exercise to the reader. I am sure the 2nd derivatives have been independently derived many times. However, I cannot find them exhibited anywhere. Is there a calculus expert out there who can detail the 2nd derivatives? Many thanks.


title: confusion regarding confidence interval of normal distribution

text: I am a bit confused when it comes to 68-95-99.7 rule of confidence interval of normal distribution. Normally I could see confidence interval of 95% for a sample statistic (s) with margin of error lets say e. So the confidence interval is s +- e. So 95% means that if we do the sampling 100 times and calculate the sample statistic s 100 times then for each we will have different s+-e. However 95 times the true population parameter p lies within s+-e. This is what I have understood about confidence interval. The confidence interval is s+-e and the confidence level is 95%. Now when it comes to normal distribution, how the above idea fits in. I mean here confidence interval is mui+-sigma and the confidence level is 68%?. According to wikipedia it is something like this About 68% of values drawn from a normal distribution are within one standard deviation σ away from the mean; about 95% of the values lie within two standard deviations; and about 99.7% are within three standard deviations I am confused where is the sample statistic here, it says the values drawn(which is directly the random variable).


title: Reusing predicted values as independent variables for new linear regression

text: **Edit** I've rephrased this question severely. Suppose I have a fitness center, and instead of a monthly fee, the people are paying for courses they are taking. I've got the monthly total revenue of all participants who started in January, for seven months (up to July). I've also collected the monthly revenue generated by participants who started in February and from those who started in March, also up to July. These cohorts don't need to have the same size. The data looks like this:                       Jan     Feb     Mar     Apr     May     Jun     Jul      Coh1    6970    9890    8610    7200    6270    4680    4860     Coh2    NA      5970    7850    6720    5720    5020    3680     Coh3    NA      NA      4350    5420    4570    4010    3460      Note that the data above is not tab-separated - the head and the NA have excessive tabs for a nicer display here. I'd like to predict the revenue generated by the second cohort in August, and for the third cohort in August and September, assuming that their training habits over the time are the same as these of the first cohort. Plotting the first row ($x_{1, {\rm Jan}}, \dots, x_{1, {\rm Jun}}$) over the second (shifted by 1 to the left, $x_{2, {\rm Feb}}, \dots, x_{2, {\rm Jul}}$), and the third (shifted by 2 to the left) looks pretty linear, so linear regression seems to be suitable here. Employing this shift, it seems reasonable to model $x_{2, {\rm Feb}}, \dots, x_{2, {\rm Jul}}$ as observations , and $x_{1, {\rm Jan}}, \dots, x_{1, {\rm Jun}}$ as independent variables to predict the August value for Coh2, $\hat x_{2, {\rm Aug}}= \beta_0 + \beta_1 x_{1, {\rm Jul}}$. However, how does one proceed with the next cohort? Intuitively, it seems like a better idea to predict the August values for Coh3 using the shifted values of Coh2 ($x_{2, {\rm Jan}}, \dots, x_{2, {\rm Jul}}$) than to go back to the less recent observations of Coh1 ($x_{1, {\rm Jan}}, \dots, x_{1, {\rm Jun}}$). While I am fine with estimating $\hat x_{3,{\rm Aug}}$ by plugging $x_{2, {\rm Jul}}$ into its corresponding equation, I am really reluctant to use $\hat x_{2,{\rm Aug}}$ to get $\hat x_{3,{\rm Sep}}$ - since I'd use a prediction (which has errors) to predict something else. Doesn't this create even more errors? TL; DR: What happens when I use a prediction instead of an actually measured value to get a prediction?


title: R limma voom function mean-variance trend

text: I am using the limma package in R to do some analysis on a count data matrix. I use the voom function and that normally creates a plot with the mean variance trend line in it. Now I created also a mean variance trend plot, but it doesn't look normal. Has anyone seen something like this before? ![enter image description here](http://i.stack.imgur.com/PuWZM.png)


title: Why is the Chow test with $R^2$ wrong?

text: The regular way to compute the F-value for a Chow forecast test is: $$F=\frac{(e_R'e_R-e_1'e_1)/g}{e_1'e_1/(n-k)}$$ My professor said something today about that a Chow forecast test using $R^2$ would be wrong, because the number of observations in the restricted model that gives $e_R$ is different from that in the model that gives $e_1$. Could anyone please explain me why this is so?


title: How to calculate whether a change in expenditure is statisically significant

text: Sorry, this is probably very basic... I have data on expenditure over 3 years. This expenditure is made up of 2 products L and H. L has been increasing and H has been decreasing, so overall L+H is about the same. Can I measure whether there has been a statistically significant decrease in H between the 3 time points? Thanks for your help, I'm fretting and can't see the wood for the trees!


title: step by step analysis for GDP forecast with PLS and SAS

text: I'm currently trying to do a forecast of GDP , although I'm new to the econometric field , with SAS and the Partial Least Square method. My question is the following: Does anyone have any articles that will give me a better understanding of what is needed in term of data preparation , the coding and any hints. I know I'm bit vague but my understanding is vague on how to use PLS in that occasion. Thank you in advance for any tips you will be able to give me.


title: How to choose nlme or lme4 R library for mixed effects models?

text: I have fit a few mixed effects models (particularly longitudinal models) using `lme4` in `R` but would like to really master the models and the code that goes with them. However, before diving in with both feet (and buying some books) I want to be sure that I am learning the right library. I have used `lme4` up to now because I just found it easier than `nlme`, but if `nlme` is better for my purposes then I feel I should use that. I'm sure neither is "better" in a simplistic way, but I would value some opinions or thoughts. My main criteria are:   1. easy to use (I'm a psychologist by training, and not particularly versed in statistics or coding, but I'm learning)   2. good features for fitting longitudinal data (if there is a difference here- but this is what I mainly use them for)   3. good (easy to interpret) graphical summaries, again not sure if there is a difference here but I often produce graphs for people even less technical than I, so nice clear plots are always good (I'm very fond of the xyplot function in lattice() for this reason). As usual, hope this question isn't too vague, and thanks in advance for any wisdom!


title: Brown-Forsythe and Welch f-ratios in two-way ANOVAs?

text: I understand that in One-Way ANOVA two alternative F-Ratios have been derived to be robust when homogeneity of variance has been violated. Tomarkin and Serlin (1986) review amongst other techniques the Brown-Forsythe and Welch F-Ratios and conclude that both control the type I error rate well. So far I have only come across W and B-Fs F Ratios in One-Way ANOVAS. Am I able to use them in ANOVAS with two factors? And if no, why not? Thanks


title: Find probability that a given number of events occur when each event has a different probability

text: Say I have three independent events with a probability of .25, .5, and .75 chance success for the individual events. How would I determine the probability that 0, 1, 2, or 3 of the events are successful?


title: How to use the LOGNORMALDIST function to generate a Cumulative Distribution Function?

text: In Excel or Google Docs you can readily construct a Cumulative Distribution Function for the Normal distribution using the NORMDIST function (X, Mean, Standard Deviation). I am trying to do the same for a Lognormal distribution using the LOGNORMALDIST function that has the same parameters (X, Mean, Standard Deviation). I have played around with using X or LN(X) on the mentioned parameters. But, I get absurd results. Do you know how to generate a Lognormal Cumulative Distribution Function using LOGNORMALDIST? Thanks in advance for any assistance. I will readily reward anyone with a good answer with a positive vote.


title: Where can I find good statistics quizzes?

text: Introductory, advanced, and even obscure, please. Mostly to test myself. I like to make sure I know what the heck I'm talking about :) Thanks


title: Combining Probability Density Functions

text: I am trying to predict the outcome of a random variable x, which is a real- valued number. In some cases I can observe another variable y1, which should approximate x. I model y1 as a Gaussian distribution with mean of 0 and an empirically estimated standard deviation, and use the probability density function of that Gaussian to predict x. In some cases I also have a second observation y2, which I can model similarly as a Gaussian. What is the appropriate way to combine y1 and y2 into an estimate on x? Should I add the distributions and model x as a mixture of Gaussians, or should I multiply them? Or, should I try both and pick the answer that maximizes the probability of sampling an independent set of data?


title: ARIMA forecasting life expectancy

text: I'm working on the forecasting of life expectancy actually. I have written code following the usual procedure. The results are not trustworthy because the life expectancy should have a positive slope (logically) but in my case it comes flat for the 50 years ahead. The ARIMA(1,1,2) is the best model with the lowest AICc. I'm in trouble with the point forecast where the straight line is flat. How do you think I can modify my code in order to get more powerful results ![enter image description here](http://i.stack.imgur.com/yLPE0.jpg) Here you find my code and the plot. The life expectancy is I(1). I have differenced it and the first difference is I(0). The first difference of life expectancy is stationary. Then I got the best model with the lowest AICc and forecast the the life expectancy level.               library(forecast)     AA<-Alberta$Male     AA1<-ts(AA,start=1921,end=2009,frequency=1)     A2<-diff(AA1,1)     fit1<-arima(A2, order=c(1,1,0))     fit2 <-arima(A2, order=c(1,1,1))     fit3<-arima(A2, order=c(1,1,2))     fit4<-arima(A2,order=c(2,1,0))     fit5<-arima(A2,order=c(2,1,1))     fit6<-arima(A2,order=c(2,1,2))     fit1     fit2     fit3     fit4     fit5     fit6     fit15<-arima(AA1,order=c(1,1,2))     ARIMA50ALBERTA<-forecast(fit15,50)     plot(ARIMA50ALBERTA)      These are the results:                Point     Forecast    Lo 80    Hi 80    Lo 95    Hi 95     2010       78.63251 77.74617 79.51885 77.27697 79.98806     2011       78.59585 77.37187 79.81984 76.72393 80.46778     2012       78.61099 77.16457 80.05740 76.39888 80.82309     2013       78.60474 76.95134 80.25814 76.07608 81.13340     2014       78.60732 76.77553 80.43911 75.80584 81.40880     2015       78.60625 76.60992 80.60259 75.55312 81.65938     2016       78.60669 76.45917 80.75422 75.32234 81.89105     2017       78.60651 76.31746 80.89557 75.10570 82.10732     2018       78.60659 76.18437 81.02880 74.90213 82.31104     2019       78.60656 76.05809 81.15502 74.70901 82.50410     2020       78.60657 75.93783 81.27531 74.52509 82.68805     2021       78.60656 75.82274 81.39039 74.34907 82.86405     2022       78.60657 75.71223 81.50090 74.18006 83.03307     2023       78.60656 75.60578 81.60734 74.01727 83.19586     2024       78.60656 75.50299 81.71014 73.86006 83.35307     2025       78.60656 75.40349 81.80964 73.70789 83.50524     2026       78.60656 75.30699 81.90614 73.56030 83.65283     2027       78.60656 75.21323 81.99990 73.41691 83.79622     2028       78.60656 75.12200 82.09113 73.27738 83.93575     2029       78.60656 75.03309 82.18004 73.14141 84.07172     2030       78.60656 74.94635 82.26678 73.00874 84.20439     2031       78.60656 74.86161 82.35152 72.87915 84.33398     2032       78.60656 74.77874 82.43439 72.75242 84.46071     2033       78.60656 74.69764 82.51549 72.62838 84.58475     2034       78.60656 74.61818 82.59495 72.50686 84.70627     2035       78.60656 74.54027 82.67286 72.38771 84.82542     2036       78.60656 74.46383 82.74930 72.27080 84.94233     2037       78.60656 74.38878 82.82435 72.15602 85.05711     2038       78.60656 74.31504 82.89809 72.04324 85.16989     2039       78.60656 74.24254 82.97059 71.93236 85.28077     2040       78.60656 74.17123 83.04190 71.82330 85.38983     2041       78.60656 74.10104 83.11209 71.71596 85.49717     2042       78.60656 74.03193 83.18119 71.61027 85.60286     2043       78.60656 73.96386 83.24927 71.50615 85.70698     2044       78.60656 73.89676 83.31637 71.40354 85.80959     2045       78.60656 73.83061 83.38252 71.30237 85.91076     2046       78.60656 73.76536 83.44777 71.20258 86.01055     2047       78.60656 73.70098 83.51215 71.10412 86.10901     2048       78.60656 73.63743 83.57570 71.00694 86.20619     2049       78.60656 73.57469 83.63844 70.91098 86.30215     2050       78.60656 73.51272 83.70041 70.81620 86.39693     2051       78.60656 73.45149 83.76164 70.72256 86.49057     2052       78.60656 73.39098 83.82214 70.63002 86.58311     2053       78.60656 73.33117 83.88196 70.53855 86.67458     2054       78.60656 73.27203 83.94110 70.44810 86.76503     2055       78.60656 73.21353 83.99960 70.35864 86.85449     2056       78.60656 73.15567 84.05746 70.27014 86.94299     2057       78.60656 73.09841 84.11472 70.18257 87.03056     2058       78.60656 73.04174 84.17139 70.09590 87.11723     2059       78.60656 72.98564 84.22749 70.01010 87.20303


title: Questions about KL divergence?

text: I am comparing two distributions with KL divergence which returns me a non- standardized number that, according to what I read about this measure, is the amount of information that is required to transform one hypothesis into the other. I have two questions: a) Is there a way to quantify a KL divergence so that it has a more meaningful interpretation, e.g. like an effect size or a R^2? Any form of standardization? b) In R, when using KLdiv (flexmix package) one can set the 'esp' value (standard esp=1e-4) that sets all points smaller than esp to some standard in order to provide numerical stability. I have been playing with different esp values and, for my data set, I am getting an increasingly larger KL divergence the smaller a number I pick. What is going on? I would expect that the smaller the esp, the more reliable the results should be since they let more 'real values' become part of the statistic. No? I have to change the esp since it otherwise does not calculate the statistic but simply shows up as NA in the result table... Thanks in advance, Ampleforth


title: Problem with ANOVA repeated measures: "Error() model is singular"

text: I need an help because I don´t know if the command for the ANOVA analysis I am performing in R is correct. Indeed using the function aov I get the following error: `In aov (......) Error() model is singular` The structure of my table is the following: subject, stimulus, condition, sex, response Example:               subject  stimulus condition sex    response     subject1    gravel  EXP1    M      59.8060     subject2    gravel  EXP1    M      49.9880     subject3    gravel  EXP1    M      73.7420     subject4    gravel  EXP1    M      45.5190     subject5    gravel  EXP1    M      51.6770     subject6    gravel  EXP1    M      42.1760     subject7    gravel  EXP1    M      56.1110     subject8    gravel  EXP1    M      54.9500     subject9    gravel  EXP1    M      62.6920     subject10   gravel  EXP1    M      50.7270     subject1    gravel  EXP2    M      70.9270     subject2    gravel  EXP2    M      61.3200     subject3    gravel  EXP2    M      70.2930     subject4    gravel  EXP2    M      49.9880     subject5    gravel  EXP2    M      69.1670     subject6    gravel  EXP2    M      62.2700     subject7    gravel  EXP2    M      70.9270     subject8    gravel  EXP2    M      63.6770     subject9    gravel  EXP2    M      72.4400     subject10   gravel  EXP2    M      58.8560     subject11   gravel  EXP1    F      46.5750     subject12   gravel  EXP1    F      58.1520     subject13   gravel  EXP1    F      57.4490     subject14   gravel  EXP1    F      59.8770     subject15   gravel  EXP1    F      55.5480     subject16   gravel  EXP1    F      46.2230     subject17   gravel  EXP1    F      63.3260     subject18   gravel  EXP1    F      60.6860     subject19   gravel  EXP1    F      59.4900     subject20   gravel  EXP1    F      52.6630     subject11   gravel  EXP2    F      55.7240     subject12   gravel  EXP2    F      66.4220     subject13   gravel  EXP2    F      65.9300     subject14   gravel  EXP2    F      61.8120     subject15   gravel  EXP2    F      62.5160     subject16   gravel  EXP2    F      65.5780     subject17   gravel  EXP2    F      59.5600     subject18   gravel  EXP2    F      63.8180     subject19   gravel  EXP2    F      61.4250     .....     .....     .....     .....      As you can notice each subject repeated the evaluation in 2 conditions (EXP1 and EXP2). What I am interested in is to know if there are significant differences between the evaluations of the males and the females. This is the command I used to perform the ANOVA with repeated measures:               aov1 = aov(response ~ stimulus*sex + Error(subject/(stimulus*sex)), data=scrd)     summary(aov1)      I get the following error:               > aov1 = aov(response ~ stimulus*sex + Error(subject/(stimulus*sex)), data=scrd)     Warning message:     In aov(response ~ stimulus * sex + Error(subject/(stimulus * sex)),  :     Error() model is singular     > summary(aov1)          Error: subject               Df  Sum Sq Mean Sq F value Pr(>F)     sex        1  166.71  166.72   1.273  0.274     Residuals 18 2357.29  130.96                         Error: subject:stimulus                   Df Sum Sq Mean Sq F value Pr(>F)         stimulus       6 7547.9 1257.98 35.9633 <2e-16 ***     stimulus:sex   6   94.2   15.70  0.4487 0.8445         Residuals    108 3777.8   34.98                        ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1           Error: Within                Df Sum Sq Mean Sq F value Pr(>F)     Residuals 420 9620.6  22.906                    >       The thing is that looking at the data it is evident for me that there is a difference between male and females, because for each stimulus I always get a mean higher for the males rather than the females. Therefore the ANOVA should indicate significant differences.... Is there anyone who can suggest me where I am wrong? Finally, I know that in R there are two libraries on linear mixed models called nlme and lme4, but I have never used it so far and I don´t know if I have to utilize it for my case. Is it the case to utilize it? If yes, could you please provide a quick R example of a command which could solve my problem? Thanks in advance! Best regards * * * Dear all, I am stuck now ;-( Indeed I understood everything you suggested me but still I don´t get significance in the ANOVA results, and definitively there is an error, because results cannot be non-significant. Indeed looking at the means for each stimulus, it is possible to notice that males gave always higher evaluations than females. To prove this I discarded for a moment the effect of the repeated measures, and I performed an ANOVA separately on both the two conditions (EXP1 and EXP2) during which the evaluations were given. What I get is significant differences between males and female, in both EXP1 and EXP2. Now, why when I perform the ANOVA with repeated measures I don´t get the same behavior? My design is the following: -sex is a between-subjects factor (with two levels) -stimulus is a within-subjects factor (with 3 assumed levels) -condition is a within-subjects factor (with 2 levels) -all factors are fully crossed I tried, both the ways suggested but without achieving significance:               mDf <- aggregate(response ~ subject + sex, data=scrd, FUN=mean)     summary(aov(response ~ sex, data=mDf))     # ANOVA with just the between-effect      and               aov1 = aov(response ~ sex*stimulus*condition + Error(subject/(stimulus*condition)), data=scrd)     summary(aov1)      Instead if I perform the ANOVA on the two subtables of EXP 1 and 2 I get significant differences.               table_EXP1 <- subset(scrd, condition == "EXP1")     table_EXP2 <- subset(scrd, condition == "EXP2")               fit_table_EXP1 <- lm(response ~ stimulus*sex, data=table_EXP1)      summary(fit_table_EXP1 )     anova(fit_table_EXP1 )               fit_table_EXP2 <- lm(response ~ stimulus*sex, data=table_EXP2)      summary(fit_table_EXP2)     anova(fit_table_EXP2)      ....how can this be possible?...it is a contraddiction.... HELP! Please enlighten me! Thanks in advance Cheers


title: Blinder-Oaxaca decomposition - Unexpained age unusually high

text: I'm doing a Blinder-decomposition (females as reference group). Is it normal that a variable has an explained or unexplained contribution above 100%? I know my calculation is correct. However, the unexplained age is 1.099 (log value); while the overall difference is 0.10. So unexplained age has an unexplained part of 1099%. The unusually high value comes from the big difference in slope. Males (0.06) vs females (0.02). So I calculate (0.06-0.02)0.04. Overall group men 7,809 Overall group women 7,706 Overall Difference 0,103 10,84% Overall Explained 0,043 Overall Unexplained 0,06 58% Explained Occupation 0,00 -1% Explained Age -0,02 -18% Explained Overtime -0,02 -15% Explained Faculty 0,08 75% Unexplained Occupation -0,01 -11% Unexplained Age 0,71 692% (here I combine age with age_squared) Unexplained Overtime -0,01 -11% Unexplained Faculty 0,00 0% Unexplained Constant -0,63 -611% -Thanks


title: Goodness of fit for a regression with multiple predictors

text: This is perhaps basic but I couldn't find a suitable reference. I have a regression model with a rather complicated link function. So $\vec{x}$ is a vector of continuous predictors, and $z$ is a binary variable such that according to the model: $Pr(z=1) = f(\vec{x})$ for some (known) function $f$. I observe data of the form $(\vec{x}^{(1)}, z^{(1)}), (\vec{x}^{(2)}, z^{(2)}), (\vec{x}^{(n)}, z^{(n)})$ and want to test the null hypothesis that the above model is the one generating the data - that is compute a statistic and reject the model if the statistic is too extreme. What would be a good goodness-of-fit test for this case? is there a 'standard' way to test for this? One possibility is binning the data points by the value of $f(\vec{x})$, (say to $10$ bins: $([0,0.1], ..[0.9,1])$ and performing a chi-square test for expected vs. observed proportion of $z$'s in each bin. Another is to bin the multidimensional space of the $\vec{x}$'s (say if $\vec{x}$ is two- dimensional, we can divide $R^2$ to $100$ squares and compute a chi-square for observed vs. expected for each square). Yet another one is not binning at all but just computing $\sum_i (z^{(i)} - f(\vec{x}^{(i)}))^2/f(\vec{x}^{(i)})$ but this seems to cause numerical issues since sometimes $f(\vec{x}^{(i)})$ is very small. Are there other known approaches? which test would be the most appropriate?


title: Modifying linear ballistic accumulator (LBA) simulation in  R

text: The "Linear Ballistic Accumulator" model (LBA) is a rather successful model for human behaviour in speeded simple decision tasks. Donkin et al (2009, PDF) provide code that permits estimating the parameters of the model given human behavioural data, and I've copied this code (with some minor formatting changes) to a gist here. However, I'd like to make a seemingly minor modification to the model but I'm not sure how to achieve this modification in the code. To start with the canonical model, LBA represents each response alternative as a competitor in a rather strange race such that the competitors can differ in the following characteristics:   * Starting position: this varies from race to race according to a uniform distribution bounded by U(0,X1).   * Speed: this is kept constant within a given race (no acceleration) but varies from race to race according to a Gaussian distribution defined by N(X2,X3)   * Finish line position (X4) Thus, each competitor has its own set of values for X1, X2, X3 and X4. The race is repeated many times, with the winner and their time recorded after each race. A constant of X5 is added to every winning time. Now, the modification I want to make is to swap the variability in the starting point to the finish line. That is, I want the start point to be zero for all competitors and all races, thereby eliminating X1, but I want to add a parameter, X6, that specifies the size of the range of a uniform distribution centered on X4 from which each competitor's finish line is sampled for each race. In this model, then, each competitor will have values for X2, X3, X4, and X6, and we still have the across-competitor value for X5. I'd be very grateful if anyone is willing to help with this. Oh, and to provide a mapping from the "X" named parameters described above to the variable names used by the LBA code I linked: X1 = x0max; X2 = driftrate; X3 = sddrift; X4 = chi; X5 = Ter.


title: Why does the adjusted r-squared of this model improve with addition of a statistically insignificant variable?

text: I stumbled on this while doing MLR, and was curious as to why this happens. The adjusted R-squared is (if I understand correctly) supposed to be a way of comparing the predictive quality of models with different numbers of explanatory variables. In the second model, I've added a statistically insignificant variable (weight), which has apparently improved the model. My only thought that this is because the point estimate of this value is not 0 - so there may be a 'significant' effect at a lower level. Is that right? Model 1:                   Model 1                                              Sum of           Mean          Source                   DF        Squares         Square    F Value    Pr > F               Model                     2     6017.30007     3008.65004      12.36    <.0001          Error                   424         103221      243.44647          Corrected Total         426         109239                                 Root MSE             15.60277    R-Square     0.0551                       Dependent Mean      120.03044    Adj R-Sq     0.0506                       Coeff Var            12.99901                                                  Parameter Estimates                       Parameter     Standard                        Variance     Variable   DF     Estimate        Error  t Value  Pr > |t|    Inflation    95% Confidence Limits          Intercept   1     75.85363      8.91793     8.51    <.0001            0     58.32478     93.38249     age         1      0.66112      0.15314     4.32    <.0001      1.00204      0.36011      0.96212     chol        1      1.86495      0.82213     2.27    0.0238      1.00204      0.24900      3.4809      Model 2 (with addition of insignificant variable):                                 Root MSE             15.58705    R-Square     0.0592                       Dependent Mean      120.03044    Adj R-Sq     0.0525                       Coeff Var            12.98591                                                  Parameter Estimates                       Parameter     Standard                        Variance     Variable   DF     Estimate        Error  t Value  Pr > |t|    Inflation    95% Confidence Limits          Intercept   1     57.69446     16.03325     3.60    0.0004            0     26.17970     89.20922     age         1      0.66180      0.15299     4.33    <.0001      1.00205      0.36110      0.96251     chol        1      2.02756      0.82993     2.44    0.0150      1.02320      0.39626      3.65885     weight      1      0.09687      0.07111     1.36    0.1738      1.02122     -0.04290      0.2366


title: Example to show that items in a strong association rule may actually be negatively correlated?

text: Give a short example to show that items in a strong association rule may actually be negatively correlated? Thanks in advance!


title: Wald test and Type III test don't agree?

text: ![enter image description here](http://i.stack.imgur.com/PFhg8.jpg) I am having trouble understanding the output here. In the green table, it seems like SAS is saying "Gender F" and "Gender M" are not significant. However, in the blue table, seems like it's saying that it would make a different if we just exclude the variable Gender. I am a bit confused because these two conclusions seem contradictory to each other. Any idea?


title: References for using networks to display correlations?

text: I just came by a post talking about networks for displaying correlations: ![enter image description here](http://i.stack.imgur.com/BhyaZ.png) Is this a known method? Can someone shed some insights into it? (I'm wondering about how useful it might be, and when.)


title: What is the curse of dimensionality?

text: Specifically, I'm looking for references (papers, books) which will rigorously show and explain the curse of dimensionality. This question arose after I began reading this white paper by Lafferty and Wasserman. In the third paragraph they mention a "well known" equation which implies that the best rate of convergence is $n^{-4/(4-d)}$; if anyone can expound upon that (and explain it), that would be very helpful. Also, can anyone point me to a reference which derives the "well known" equation?


title: What is the variance of the maximum of a sample?

text: I'm looking for bounds on the variance of the maximum of a set of random variables. In other words, I'm looking for closed-form formulas for $B$, such that $$ \mbox{Var}(\max_i X_i) \leq B \enspace, $$ where $X = \{ X_1, \ldots, X_M \}$ is a fixed set of $M$ random variables with finite means $\mu_1, \ldots, \mu_M$ and variances $\sigma_1^2, \ldots, \sigma_M^2$. I can deduce that $$ \mbox{Var}(\max_i X_i) \leq \sum_i \sigma_i^2 \enspace, $$ but this bound seems very loose. A numerical test seems to indicate that $B = \max_i \sigma_i^2$ might be a possibility, but I have not been able to prove this. Any help is appreciated.


title: model for continuous dependent variable bounded between 0 and 1

text: I'm attempting a multiple regression model where the predicted variable is runoff ratio - the ratio of watershed discharge to the precipitation input. This should generally be bounded [0,1], however, due to measurement error some values > 1 occur. Originally, I modeled this with the predicted variable un-transformed, but logistic regression has been suggested to me, I also have heard Beta regression suggested. I'm not sure how to proceed, and if these transformations are appropriate to my data: ![enter image description here](http://i.stack.imgur.com/hCOuF.jpg) My questions are: 1) Is a logistic regression appropriate for these data? and 2) If I were to proceed with logistic regression, would I need to convert the runoff ratios to proportions, or would I apply the logit to the values as they are? Sorry if these are obtuse questions - I'm new to logit and most of the information I have found is for binary response variables. **Edited for suggested additions:** As a simple version: I am modeling runoff ratio (rr) as an effect of precipitation (pcp) and antecedent water table position (ant): rr ~ pcp + ant rr is a continuous variable. I am not interested in the _probability_ of specific values, rather I'm interested in the values themselves - both to assess the significance of the predictors and as a predictive model. Conceptually, I was fine modeling it un-transformed. However, a simple linear regression allows predicted values outside of the physical range of [0,1]. As mentioned above, measurement error does lead to values >1, which I'll eventually have to deal with.


title: Variables occurring after target occurence, can they be used in regression model?

text: My study is to judge the appropriateness of target=1, not for predictive purpose. Can the variables occurring after target occurrence be used in regression or decision tree model? For example, a patient is admitted to ICU, he/she is likely to receive some treatments that often occur in ICU. Can I use the presence or not presence of those treatments(and other variables) to judge if the ICU admission is appropriate. I think the base of regression and decision tree is chi-sq test that describes the association between target and other variables, causal-effect relationship is not required if not for prediction. Am I right? Thanks


title: when to use panel time-series regression vs seeming unrelated regression

text: I am a bit confused on whether or not I have to use a fixed-effect panel time- series method or SUR (seemingly unrelated regression). To get a background of what I am trying to do, I have 10 panels of 25 weeks of data with four independent variables and one dependent variable and I am trying to find how these four independent variables effect the dependent variable. I am currently using R to do my analysis.


title: How to find confidence intervals for ratings?

text: Evan Miller's "How Not to Sort by Average Rating" proposes using the lower bound of a confidence interval to get a sensible aggregate "score" for rated items. However, it's working with a Bernoulli model: ratings are either thumbs up or thumbs down. > What's a reasonable confidence interval to use for a rating model which > assigns a discrete score of $1$ to $k$ stars, assuming that the number of > ratings for an item might be small? I think I can see how to adapt the centre of the Wilson and Agresti-Coull intervals as $$\tilde{p} = \frac{\sum_{i=1}^n{x_i} + z_{\alpha/2}^2\; p_0}{n + z_{\alpha/2}^2}$$ where either $p_0 = \frac{k+1}{2}$ or (probably better) it's the average rating over all items. However, I'm not sure how to adapt the width of the interval. My (revised) best guess would be $$\tilde{p} \pm \frac{z_{\alpha/2}}{\tilde{n}} \sqrt{\frac{\sum_{i=1}^n{(x_i - \tilde{p})^2} + z_{\alpha/2}(p_0-\tilde{p})^2}{\tilde{n}}}$$ with $\tilde{n} = n + z_{\alpha/2}^2$, but I can't justify with more than hand-waving it as an analogy of Agresti-Coull, taking that as $$\text{Estimate}(\bar{X}) \pm \frac{z_{\alpha/2}}{\tilde{n}} \sqrt{\text{Estimate}(\text{Var}(X))}$$ Are there standard confidence intervals which apply? (Note that I don't have subscriptions to any journals or easy access to a university library; by all means give proper references, but please supplement with the actual result!)


title: Variance of product of dependent variables

text: What is the formula for variance of product of dependent variables? In the case of independent variables the formula is simple: $$ {\rm var}(XY) = E(X^{2}Y^{2}) - E(XY)^{2} = {\rm var}(X){\rm var}(Y) + {\rm var}(X)E(Y)^2 + {\rm var}(Y)E(X)^2 $$ But what is the formula for correlated variables? By the way, how can I find the correlation based on the statistical data?


title: How do I solve $E\left[ E \left(X|Z \right) E\left( Y|Z \right)\right]$?

text: I am trying to solve $E\left[ E \left( \mathbf{X}|\mathbf{Z} \right) E \left( \mathbf{Y}|\mathbf{Z} \right) \right]$, (where $\mathbf{X}$, $\mathbf{Y}$, and $\mathbf{Z}$ are random variables) but I am confused - should the solution be $E \left( \mathbf{X}\mathbf{Y} \right)$ or $E\left(\mathbf{X} \right)E\left(\mathbf{Y}\right)$. Any explanation as to why this way or the other will be very helpful. Thanks.


title: What is the upside of treating a factor as random in a mixed model?

text: I have a problem embracing the benefits of labeling a model factor as random for a few reasons. To me it appears like in almost all cases the optimal solution is to treat all of the factors as fixed. First, the distinction of fixed vs random is quite arbitrary. The standard explanation is that, if one is interested in the particular experimental units per se, then one should use fixed effects, and, if one is interested in the population represented by the experimental units, one should use random effects. This is not of much help because it implies one can alternate between fixed and random views even if the data and experimental design remain the same. Also, this definition promotes an illusion that, if a factor is labeled as random, the inference drawn from the model is somehow more applicable to the population than in case when the factor is labeled as fixed. Finally, Gelman shows that the fixed-random distinction is confusing even at the definition level because there are four more definitions of what fixed and random effects are. Second, the estimation of mixed models is quite complicated. As opposed to a "purely fixed” model, there are more than a few ways to obtain the p-values. Prof. Bates who implemented REML estimation in the lme4 package in R went so far as to refuse to report the p-values altogether. Third, there is a murky issue of how many implicit parameters are introduced by a random factor. The following example is my adaptation of that in Burnham & Anderson, _Model Selection and Multi-Model Inference: A Practical Information-Theoretic Approach_. From the bias-variance tradeoff perspective, the role of random effects can be illustrated as follows. Consider a one-way ANOVA with $K$ treatments and $K$ main factor effects, of which $K - 1$ are estimable. The error term has $\mathcal N(0, \sigma^2)$ distribution. If the number of observations is fixed, the bias-variance tradeoff will deteriorate as $K$ goes up. Suppose we say that the $K$ main effects are drawn from $\mathcal N(0, \sigma_K)$ distribution. The corresponding model will have a complexity that is somewhere in between the fixed (overfitted) version and the underfitted model that contains the intercept only. The number of effective parameters in the fixed model is $$1 \:\:\mathrm{intercept} + (K - 1) \:\:\mathrm{main\: effects} + 1 \:\:\sigma = K + 1.$$ The number of effective parameters in the random model is at least three: $ \mathrm{intercept}, \sigma, \sigma_K$. In addition, the random model has a number of “hidden” parameters implied by the distributional (normal in this case) restriction imposed on the main effects. In particular, if there is a factor with two levels, it doesn’t make sense to call it random, even if we know for sure that its levels have been sampled at random from some population. That is because the fixed effect version has three parameters, and the random effect version has over three parameters. In this case, the random model turns out to have more complexity than the fixed version. Apparently, a switch from the fixed to random version is more grounded for larger $K$. However, the number of “hidden” parameters in the random model is unknown, so it is impossible to compare the fixed and random versions based on the information criteria such as AIC. Therefore, while this example illuminates the contribution of random effects (the possibility of a better bias-variance tradeoff), it also shows that it is hard to say when it is justifiable to relabel the factor from fixed to random. None of the above problems are present in a “purely fixed” model. Therefore, I am willing to ask:   1. Can anyone provide an example when something very bad happened when a random factor was used as if it were fixed? I believe there should be some simulation studies that address the issue explicitly.   2. Is there a proven quantitative method to decide when it makes sense to switch from fixed to random label?


title: Find interquartile range of these data

text: Take these data:               5, 8, 9, 14      R says the interquartile range is `3`:               IQR(c(5, 8, 9, 14))     # 3      ...but I make it `5`. What am I doing wrong? Here's the steps I've taken:   1. Find median, which is `8.5`   2. Split the data around the median into two groups, like this: `(5, 8), (9, 14)`   3. Find the median of these two groups, which is `6.5` and `11.5`   4. Subtract `6.5` from `11.5`, which is `5`


title: Why use the Mahalanobis distance

text: I understand in theory why the Mahalanobis distance is a good measure for mutlivariate outlier detection. However, everything I tend to read warns against calculating the inverse/pseudoinverse of a covariance matrix, which is needed to compute the mahalanobis distance. So, if nobody wants to compute the inverse, what distance measure should be used?


title: What is coskewness and how can it be calculated?

text: I would like to calculate coskewness of two random variables. However I couldn't find even basic information on this matter. Is there a standard definition? How to calculate it? If not what are my alternative options? Is it possible to have a normalized coskewness like correlation coefficient? Thanks


title: when should I need to be concerned about autocorrelation?

text: I have a series of wave heights observations (2observations per hour for each day in 7 months) and I'm trying to model a regression with wind data (same frequency of observation). Since waves observations, which are my response variable, have been taken continuously trought time, I started by thinking that this one of those case where your observations are correlated one another. however In the book (Experimental Design and Data Analysis for Biologists. Quinn, G.P. & Keough, M.J. (2002) Cambridge University Press, New York.) I read this passage at pag 93: "...A related situation is when we have a longer time series from one or a few units and we wish to fit a model where the predictor variable is related to a temporal sequence, i.e. a time series study (Diggle 1990). Error terms and Y-values that are non-independent through time are described as autocorrelated. A common occurrence in biology is positive first-order autocorrelation, where there is a positive relationship between error terms from adjacent observations through time, i.e. a positive..." here they specify that the predictor variable should be related to time. In my case my predictor variable are x and y wind components (basically summarizing wind speed and direction), maybe water depth, but I do not have any time variable in my formula. In this case do I have to take care of possible autoregressive structures? in more general terms: autocorrelation can be a problem ONLY when at least one of your predictors is some form of time (e.g. days, hours or month, years)? If the answer to this question is 'yes' (i.e. only when...), what happen if I introduce as a random effect a variable 'month'? do I have to be worried about Autocorr now or not since is not in the fixed part of the formula? thanks Filippo


title: How to choose negative training sample for Classification problem

text: Choosing positives sample is a relative straightforward task, but I'm having some problem on determine what should I use for the negative example. I'm working on a SVM binary classificator, trying to learn whether a pair of usernames belongs to the same individual. Is obvious what to use as a positive example, I've a set of usernames pairs I know belong to the same person. I'm shuffling randomly those same username to obtain my negative example, but i'm not sure whether is a good choice or not. Do you have any suggestion? Should I randomly create two string as fake usernames? Thank you


title: How many observations are taken for a pseudo-sample when bootstrapping?

text: I am trying to use function the `boot` function (in the R package, `boot`) and we want to change how many observations are resampled each iteration of the bootstrap. If it's not possible to change the number of observations in the pseudo- sample, how many observations are used when resampling?


title: Numerical Instability of calculating inverse covariance matrix

text: I have a 65 samples of 21-dimensional data (pasted here) and I am constructing the covariance matrix from it. When computed in C++ I get the covariance matrix pasted here. And when computed in matlab from the data (as shown below) I get the covariance matrix pasted here Matlab code for computing cov from data:               data = csvread('path/to/data');     matlab_cov = cov(data);      As you can see the difference in covariance matrices are minute (~e-07), which is probably due to numerical problems in the compiler using floating point arithmetic. However, when I compute the pseudo-inverse covariance matrix from the covariance matrix produced by matlab and the one produced by my C++ code, I get widely different results. I am computing them in the same way i.e.:               data = csvread('path/to/data');     matlab_cov = cov(data);     my_cov = csvread('path/to/cov_file');     matlab_inv = pinv(matlab_cov);     my_inv = pinv(my_cov);      The difference is so huge that when I am computing the mahalanobis distance from a sample (pasted here) to the distribution of the 65 samples by: $(65/64^2) \times ((sample-mean)\times {\sum}^{-1} \times (sample-mean)')$ using the different inverse covariance matrices (${\sum}^{-1}$) I get widely different results i.e.:                (65/(64^2))*((sample-sample_mean)*my_inv*(sample-sample_mean)')     ans =             1.0167e+05          (65/(64^2))*((sample-sample_mean)*matlab_inv*(sample-sample_mean)')     ans =            109.9612      Is it normal for the small (e-7) differences in covariance matrix to have such an effect on the computation of the pseudo-inverse matrix? And if so, what can I do to mitigate this effect? Failing this, are there any other distance metrics I can use that do not involve the inverse covariance? I use the Mahalanobis distance as we know for n samples it follows a beta distribution, which I use for hypothesis testing Many thanks in advance **EDIT: Adding C++ code for calculating covariance matrix below:** The `vector<vector<double> >` represents the collection of rows from the file pasted.               Mat covariance_matrix = Mat(21, 21, CV_32FC1, cv::Scalar(0));         for(int j = 0; j < 21; j++){             for(int k = 0; k < 21; k++){                 for(std::vector<vector<double> >::iterator it = data.begin(); it!= data.end(); it++){                     covariance_matrix.at<float>(j,k) += (it->at(j) - mean.at(j)) * (it->at(k) - mean[k]);                 }                 covariance_matrix.at<float>(j,k) /= 64;              }         }


title: Time series prediction - what is Autoregressive Tree model ? (Python)

text: Our problem: model evolution of values of a continuous variable over time. I came through a paper presenting an approach for predicting the next values for a time series. Whereas ARIMA model is more accurate for long term prediction, ARTXP model is preferred to infer the next values. Microsoft library for Data Mining algorithms implements ARTXP, a variation of Autoregressive Tree model. How does algorithm works? What is a Python implementation for this model ?


title: Is it possible to find out how unlikely a given $\mathrm{Beta}(\alpha,\beta) – \mathrm{Beta}(\alpha’,\beta’)$ is?

text: Let $X_i$ be the distance between two nodes in a graph. $X_i \sim \mathrm{NB}(r,p)$ where $p \sim \mathrm{Beta}(\alpha, \beta)$. After the posterior hyperparameters are obtained, $\mathrm{Beta}(\alpha,\beta) – \mathrm{Beta}(\alpha’,\beta’)$ could be found which would correspond to the change in expected distance. Is it possible to find out how unlikely a given $\mathrm{Beta}(\alpha,\beta) – \mathrm{Beta}(\alpha’,\beta’)$ (change in expected distance) is?


title: Is there a conventional definition of $\bumpeq$ in statistics?

text: I was reading a paper on Bayesian curve fitting (Dimatteo et. al. Bayesian curve-fitting with free-knot splines, 2001) and came across the symbol $\bumpeq$. It is used a few times throughout the paper but never explicitly defined. After a few google and stackexchange searches, it appears as though the symbol is neither widely used nor conventionally defined. Below I give an example with context from the cited paper. I apologize in advance for not defining any of the other symbols, but doing so would amount to copying over large portions of text from the paper I have linked to and would be of little use to the question. From p1059 (equation 8): > Incidentally, we can also see this in the likelihood ratio approximation for > the normal model in equation (6) by > > > $$\frac{p(y|k^c,\xi^c)}{p(y|k,\xi)}\bumpeq\frac{1}{\sqrt{n}}\left(\frac{(y-B_{k,\xi}\hat{\beta})^T(y-B_{k,\xi}\hat{\beta})}{(y-B_{k,\xi^c}\hat{\beta^c})^T(y-B_{k,\xi^c}\hat{\beta^c})}\right)^{n/2}=exp(-\text{BIC}/2)$$ From context it seems that $\bumpeq$ represents an approximation. If this is case, then is it synonymous with more conventional symbols for an approximation like $\approx$ or $\sim$? or is it being used to represent a particular kind of approximation for which $\approx$ or $\sim$ would be insufficient or misleading?


title: In multiple regression, a coefficient gets opposite sign in subsample than in the full sample

text: In a multiple regression, the regression coefficient of the two subsamples is opposite to that of the full sample. Both coefficients are significant. The sum of the number of obs in the two subsamples is equal to that of the full sample. Is this normal? What are the possible causes? Thanks.


title: Fidelity of measurement using conditional probabilities

text: Let me begin by saying that I'm not entirely sure if this is the correct forum, or if Mathematics would be more suitable. The problem I'm about to describe is statistical in nature, so I suppose it fits here. I have a system $m_I$ that is either in state 1, 0 or -1. In order to determine which state it is in, I measure a different (related) system $N$ times, and I count how many 'positive occurances' $k$ I get. The chance of a positive occurrence is $p_i$, so this is binomially distributed. Depending on the state $m_I$, the there is a different $p_i$, so in principle I should measure three binomial distributions. These three binomial distributions have two intersections, $k1$ and $k2$. In order to decide what state of $m_I$ I'm dealing with, I simply see what distribution it was most likely to belong to by seeing on what side of the values $k_1$ and $k_2$ my measurement is. I want to optimize the fidelity of being correct about $m_I$ which I call $F_{avg} = \frac{1}{3}(F_1+F_0+F_{-1})$, which is defined using the cumulative distribution function of the binomial distribution $ P(X \leq k) = \textrm{Bincdf}(k,n,p) = \sum_{i=1}^{k} \left (\begin{matrix} n \\ i \end{matrix}\right )p^i(1-p)^{n-i} $ and $ F_{1} = 1 - P(X \leq k_2) $ $ F_{0} = P(X \leq k_2) - P(X \leq k_1) $ $ F_{-1} = P(X \leq k_1) $ The fidelities are thus simply the chance that a certain $k$ belongs to the interval of one of the three binomial distributions. Here the probabilities $p$ are fixed, but $k1,k2$ vary when I vary $n$. Now, $k_1$ and $k_2$ I can numerically calculate for a specific $N$, and I can thus also just evaluate $F_{avg}$ for various $N$. This shows, as is to be expected, that it converges to 1 very quickly as the distributions become more narrow for higher $N$. But here's where the problem comes in. What I have not included in this model is that every time I make a measurement (so I make $N$ in total), there is a chance $P$ that something goes wrong and that my approach is no longer valid. Basically, by measuring I have a chance of changing $m_I = 1$ into $m_I = 0$, and the other permutations. Now, for simplicity lets say that 1 can only go to 0, 0 to -1 and -1 to 1, all with the same chance P. My problem is, how do I include this in my fidelities. What I first thought of was just simply multiplying each $F$ by a factor of $(1-P)^N$, which is the chance that the system has not changed in N measurements. But my feeling is that this is not correct, and that the approach is a little more intricate. Intuitively the notion of a convolution comes to mind, but I'm not sure how applicable that is? I apologize if my story is vague. I've tried rewriting it about 4 times now, and this is as legible as I can currently make it, but that's also because I know the situation very well. So if there's any part that is particularly unclear, please let me know and I'll try to rephrase it!


title: Including confounders in a model

text: Suppose that you are performing a linear regression examining the main effect $x_1$ and want to adjust for possible confounders $x_2, x_3, x_4$. Is it better to have an unadjusted model and a model adjusted for **all** potential confounders? Or should you also consider models adjusted for only some of the confounders (e.g. $x_2$, $x_2$ and $x_3$, etc.)?


title: How would you visualize a segmented funnel? (and could you do it with Python?)

text: I saw this post on Moz which presented a segmented marketing funnel: ![enter image description here](http://i.stack.imgur.com/t9Y4g.gif) This kind of thing would have quite a bit of value in my job. What I have no idea is how to visualize raw data to show a segmented funnel like this one. The idea is that sales leads come from different sources (which we use to segment the data by) and go through several stages by the time they convert to a deal. From each stage to another some drop off. The width of each slice is determined by the absolute number of leads in each. [ **EDIT** : Notice the image used here for reference is misleading when it comes to the numbers specified on the right of each slice. There appears to be no relationship between the width of the slice and the number. The image should only be taken as a reference to the design of the segmented funnel]. Anyway, any idea how to visualize it? If possible, I'd love to have a way to do so in Python. Here's a Google Doc with some dummy data if anybody needs some... Looking forward to your insights. Thanks!


title: Are clusters in social sciences globular?

text: Many clustering algorithms require globular (gaussian) clusters. Main example is k-means. If clusters are not globular, these methods can bring wrong results. In social sciences, data is often clustered by these algorithms (market segmentation is very common example). My question then is whether there are some reasons why we should regard clusters from social surveys as globular (and therefore appropriate for methods like k-means). An example of social data that could be clustered is mentioned in my previous question. I asked my students to rate selected interest (1-5 scale) and the task then could be to cluster students according to similar interests.


title: Instrumental variable Tobit in R

text: I have a data generating process of the form:               res1 <- rnorm(N);     res2 <- res1*0.5 + rnorm(N)     x    <-  z[,1]*2 + res1;     ys   <- x*b + res2;     d    <- (ys>0); #dummy variable     y    <- d*ys;      (blatantly stolen from an old RHelp thread, where the question remained unanswered back in 2006) That is, my independent variable `x` is endogenous but there's an exogenous instrument `z`, and my dependent variable `y` is censored below at zero. Are there any ways I can consistently estimate this model in R? What I'm looking for is essentially an instrumental variable Tobit (as implemented in Stata's ivtobit command).


title: Predicting ecommerce variables based on visitation data - R or Python

text: I have 1.5 years worth of e-commerce data (orders and revenue) for 2 countries (US and UK) with daily metrics for visits, unique visitors, and page views. I would like to model these transactions and apply the model to 1.5 years worth of visitation data (visits, unique visitors, and page views) to predict e-commerce for 10 additional countries. I tried to use the `glm()` function in R, but from what I gather, the outcome variable has to be a binomial. What would be the best way to solve this problem in R or even Python?


title: Understanding covariance of errors in regression

text: I am having a hard time understanding the elements of an error covariance matrix for a class. Can someone clarify? First, the diagonal. The variance is $E(e_i^2) - E(e_i)^2$. $E(e_i) = 0$, so it's just $E(e_i^2)$. Ok . . . so that is the square of the actual value of $e_i$? $y_i - ŷ_i$? A non-zero number, i.e., how much that error value "varies" from the expected zero? So, how does the E work in this? The off-diagonal terms = $E(e_ie_j)$? The expected value of their product? What is that telling you? My notes say: "The graph suggests that the variance around the regression line is constant: $\sigma^2 = \sigma_i^2$ " ![enter image description here](http://i.stack.imgur.com/Eoe0K.png) I can't see it. Where is the constant? The ei values are all over the place, though their mean/expected value does look like zero. And why would the off- diagonal terms be zero? I would expect constant error variance to look something like: ![enter image description here](http://i.stack.imgur.com/81Yu6.png)


title: Reference needed - Who first introduced linear models, fixed and random effects models?

text: I am writing an essay which briefly discusses linear models as well as models with fixed and random effects. I am googling since ages to find a reference which says who first introduced linear models and especially later mixed effects models. I found on wikipedia that in 1918 Ronald Fisher introduced random effects models to study the correlations of trait values between relatives.[1] Also, as far as I can see Laird, Nan M.; Ware, James H. (1982). "Random-Effects Models for Longitudinal Data". was the first one who applied those models to longitudinal data. My data is also of type longitudinal and I therefore can use that reference later on. However, it would be great to found evidence of who first introduced linear models as wells the fixed effets models.


title: Prediction of a binary variable

text: I am establishing a model for prediction of a binary variable (Yes/No) depending on three continuous variables ($A$,$B$,$C$). I applied logistic regression analysis for a learning dataset vith the Tanagra software, and the results were good with high prediction accuracy. My question is: is it possible to get the probabilities of the prediction via logistic regression? Something like (0.7 Yes)? If not, what test do I have to use to get such a result?


title: Better feature selection for disease prediction

text: I am doing disease classification (i.e. a person is classified as normal or abnormal) using naive Bayes and SMO classifiers. I have around 30 attributes. Out of these I need to select the most relevant attributes. I am trying to choose the best attributes using Information Gain. Is this a better method? Please help me out. Thanks in advance.


title: LASSO vs forward selection

text: I have two questions:   1. I use cross validation to select a LASSO model, does the step in which a particular variable enter, indicate its relative importance?   Let's age enter in step 1 and sex enter in step 2, country enter in step 3. Can I make the claim that age is more important than sex , which in turn is more important than country. (Sorry about my bad language, I wonder what would be the precise statement I can make here.)   2. How does the procedure in question 1 compare to forward selection, and what is its theoretical advantage of weakness?


title: Must I normalize inputs into a perceptron that uses a sigmoid activation function?

text: I am building a neural network. Each perceptron in the network uses a sigmoid activation function. Must I normalize my inputs (which currently range form 0 to 1200)? I ask this because the sigmoid function approaches 1 as the input to it approaches infinity. Hence, if my input is not between 0 and 1, I'm afraid that the sigmoid will always return 1.


title: What method should I use for this optimization / feature selection project

text: I'm going to describe a problem and I'm not sure how to best solve it. I will describe the situation. When answering please recommend a method and maybe a software library. I'm using Python for my data analysis. I'm trying to figure out what makes people on time. So the input data would be a bunch of variables that describes the person. Some of these variables are continuous like age, some are discrete categories like sex and method of transportation. The output result is the amount of time this person is early or late expressed as an error term vs. expected time of arrival. `0` is exactly on time. `+5` is 5 minutes late. All the data is cross sectional and independent (I think) from each other. I would like to test interaction between variables, say being `male` and `college educated`, etc. How do I go about figuring out what makes people late? I've been told to use `scikit-learn`. If that is appropriate, could you please show me a little code to get me started? The example in `scikit-learn` is a univariate example. I was hoping someone could post a simple multivariate example. The 2nd part of this question is on how to predict who's going to be late, but less important.


title: Multiplying a matrix by a scalar which has a prior distribution in OpenBUGS

text: So I am having a problem specifying my model in OpenBUGS. A set of vectors in a linear regression model is given a multivariate normal prior with a constant mean vector and a constant precision matrix which is multiplied by a parameter with a hyperprior of a gamma distribution. The code is as follows:               model{     #likelihood     for(i in 1:93){         BegSal[i]~dnorm(mu[i],tau)         mu[i]<-b[1]+b[2]*Sex[i]+b[3]*Educ[i]+b[4]*Exper[i]+b[5]*Time[i]     }     #priors     tau~dgamma(.00001,.00001)     covNew[1:5,1:5]<-tau*covPri[1:5,1:5]     b[1:5]~dmnorm(meanPri[1:5],covNew[1:5,1:5])     }          #data     list(meanPri= c(13625,4500,-14388.3333,55,308.3333),          covPri = structure(.Data = c(362.5,362.5,452.5,452.5,362.5,                                       362.5,365,452.5,452.5,362.5,                                       452.5,452.5,565,565,452.5,                                       452.5,452.5,565,25565,452.5,                                       362.5,362.5,452.5,452.5,1802.5),                                       .Dim = c(5, 5)) )      I'm receiving an error "Expected Multivariate Node," when I try and compile. If anyone has any suggestions I would greatly appreciate it. I left the rest of the data out as it is fairly long. Thanks for any help.


title: Standard error for a statistic obtained via simulation

text: I am simulating the effect of certain conditions on estimates obtained using OLS regression. In running 100 replications I get 100 sets of standard errors for the regression coefficients. How do I calculate the standard error of the estimates? Is it simply the standard deviation of the estimated standard errors for a given coefficient divided by the square root of the number of replications? Thanks in advance.


title: Nominal response - censored predictor

text: What's the best way to study the association between nominal response (3 non- ordered categories) and one ordinal predictor variable (most observations are interval scaled, but some are censored (only known to be above x)) in smallish sample (50). I think that Mantel-Haenszel row means score test (SAS) might be appropriate, but it's rarely used and equivalent to Kruskal-Wallis test and I would be more inclined to use Kruskal-Wallis test if my response was ordinal and predictor nominal, not the other way around. Do I have better options than Kruskal- Wallis test?


title: Classification of cluster-correlated data

text: I'm dealing with a supervized binary classification issue. My dataset is composed of 1500 individuals, living in 600 households. I have approximately 4000 variables to classify my subjects as "infected/uninfected". I was wondering how would it be possible to account for the hierarchical nature of my data in a data mining classification method, such as CART, MARS or other methods, as it is done for instance in mixed-effects models ? I suppose that the hierarchical structure of the data cannot be ignored, because the risk of a individual to be infected is higher is there is already an infected individual in his household. Thank you Yohann


title: Plot explicit cdf instead of ecdf in R

text: I have adjusted the parameters (lambda, mu, sigma) for a mixture of two normals fitted to my data. Now I would like to plot the cdf of this model using the explicit function instead of the ecdf. Is there any way to do this or I do I have to simulate data so then I can use again ecdf? The explicit function is something like:               ipc_values_EM\$lambda[1] * dnorm(x, ipc_values_EM\$mu[1], ipc_values_EM\$sigma[1])      +     ipc_values_EM\$lambda[2] * dnorm(x, ipc_values_EM\$mu[2], ipc_values_EM\$sigma[2])       (as you can note, is the mixture of two normals different mus and different sigmas)


title: Performance metric for algorithm predicting probability of low probability events

text: I am writing a program to predict the click through rates of online ads. Two important notes about this problem:   * click through rates are very small (like 0.1%)   * click through rates depend on several parameters (like size of the ad, country in which ad is shown, whether the user has seen this ad earlier etc) As I employ statistical/machine learning techniques, I would like to measure the performance of my techniques on historical data. I cannot employ metrics like precision or accuracy since all I predict is probability rather than predicting occurrence or non-occurrence of event. In such a case, how should I measure the performance of my algorithm?


title: Is there a way to maximize/minimize a custom function in R?

text: I'm trying to minimize a custom function. It should accept five parameters and the data set and do all sorts of calculations, producing a single number as an output. I want to find a combination of five input parameters which yields smallest output of my function.


title: How to calculate confidence interval of the mean age of a sample?

text: How do I establish the +/- in terms of a set of mean ages. My mean age is 27.2. The ages are 20 23 24 43 22 26 18 32 18 41 22 20 26 46 21 27 19 19 39 40 19 39 18 38 24 24 23 30. Thank you.


title: Best imputation method for stochastic noisy data?

text: What is the best imputation method for a dataset consisting of stochastic data? For example, let's say you have a table of security returns. In some cases the missings are random, in other cases are not. For example, a new IPO would have a relatively short time-series. Suggestions:   * Expectations-Maximization algorithm   * Nearest neighbor algorithm   * Linear interpolation (I would want to avoid this since it is deterministic) Thanks!


title: Comparing backtesting returns with real trading returns

text: I have 10 years of backtested simulated performance of some trading strategy (using historical prices), and N months of actual trading performance. What statistical test can I do find out if I'm on target with the backtesting numbers? (both, in terms of expected annual returns and expected annual Sharpe ratio)


title: MCMC when the density involves integration over a simplex

text: I have the following setup.   * Parameters $W$ with density $\pi(w)$.   * Observed data $X_1,...,X_n$ iid.   * Density of $X_i|W=w$ is $f(x_i|w) = \int_{\Delta(x_i)} f(\mathbf c|w) \,d\mathbf c$.   * The simplex $\Delta(x_i) = \{\mathbf c \geq 0 : c_1 + \cdots + c_m = x_i\}$. I want to find the distribution of $W|X_1=x_n,\dots,X_n=x_n$. Although $f(\mathbf c|w)$ is a fairly cheap to evaluate, I don't think it is possible (after a lot of effort!) to get a closed or simple form for $f(x_i|w)$. Therefore it's also difficult to find out much about the posterior distribution of $W$, analytically, either. Instead I have been considering MCMC for this purpose. It is straightforward (assuming I can find a good enough $q$) to write down a basic Metropolis- Hasting algorithm:   * Come up with a proposal distribution $q(w|w')$ and pick $w_0$.   * For $t=1,2,3,\dots$:      * Sample $w$ from a $q(\cdot|w_{t-1})$.     * Let $r = \frac{ q(w|w_{t-1}) \pi(w) \prod_{i=1}^n f(x_i|w) } { q(w_{t-1}|w) \pi(w_{t-1}) \prod_{i=1}^n f(x_i|w_{t-1}) }$.     * Let $w_t = w$ with probability $\min(r,1)$, and let $w_t = w_{t-1}$ otherwise. The trouble is that, as was already mentioned above, the only way I know to evaluate $f(x_i|w)$ is via a (multi-dimensional) numerical integration. And in the above algorithm, this has to be done $2n$ times per iteration! So my $\textbf{questions}$ are, in order of specificity:   * Is there a way to combine the numerical integration required for $f(x_i|w)$ and the overall MCMC algorithm? (Intuitively, suppose we use an MCMC algorithm to evaluate each integral $f(x_i|w)$, then there would be two levels of nested Markov chains, and perhaps there is a way to combine these chains to get better overall performance.)   * Is there a (obvious, standard, or otherwise) better way to solve my problem than the basic MH algorithm above? Thanks in advance for any replies! [P.S. In case it matters, $f(\mathbf c|w)$ is the density of a probability distribution on $\mathbb R^m$.]


title: How to visualize the true dimensionality of the data?

text: I have a dataset that's nominally 16-dimensional. I have about 100 samples in one case and about 20,000 in another. Based on various exploratory analyses I've conducted using PCA and heat maps, I'm convinced that the true dimensionality (i.e. the number of dimensions needed to capture most of the "signal") is around 4. I want to create a slide to that effect for a presentation. The "conventional wisdom" about this data, which I'm looking to disprove, is that the true dimensionality is one or two. What's a good, simple visualization for showing the true dimensionality of a dataset? Preferably it should be understandable to people who have some background in statistics but are not "real" statisticians.


title: Difference between Hausdorff and earth mover (EMD) distance

text: I have two dataset that i want to compare. each dataset contain the weight of 10 different person measured for 3 different day. I am interested in measuring the probabily that the two sample originate from the same population. People seem to suggest doing a Kolmogorov-Smirnov test but i need a measurement. I was thinking doing the EMD to compare the distribution for each day EMD(dataset1.day1,dataset2.day1) + EMD(dataset1.day2,dataset2.day2) + EMD(dataset1.day3,dataset2.day3) where dataset1.day1 is the histogram of the value for day1 in dataset 1... But i could probably take each person as a 3d datapoint and do the EMD in 3d. One other possibility was to do the Hausdorff distance but doing the average of the distance for each point instead of taking the maximum distance. The two dataset have very different skewness so i was also considering using the Mann-Whitney-Wilcoxon_test. What are the main difference between the two technique.


title: Distribution of a random segment on a string

text: I have a linear string of unit length, and I randomly sample two locations a and b from Uniform(0, 1). Then I cut the string at these two locations to get a sub-string. What is the distribution for the size of the sub-string (i.e. |a-b|)? Otherwise, I want to at least know the mean and variance.


title: Nest-Varying Parameters in Nested Logit Model Estimation (mlogit package)

text: I am trying to replicate a Stata manual example (-mlogit-) of nested logit estimation,               webuse restaurant     nlogitgen type = restaurant(fast: Freebirds | MamasPizza, family: ///      CafeEccell | LosNortenos | WingsNmore, fancy: Christophers | MadCows)          nlogittree restaurant type, choice(chosen)     nlogit chosen cost distance rating || type: income kids, base(family) ///      || restaurant:, noconst case(family_id)      The output from this is as below:               . nlogit chosen cost distance rating || type: income kids, base(family) ///     >  || restaurant:, noconst case(family_id)          tree structure specified for the nested logit model           type    N      restaurant    N   k      -------------------------------------      fast   600 --- Freebirds    300  12                  +- MamasPizza   300  15      family 900 --- CafeEccell   300  78                  |- LosNortenos  300  75                  +- WingsNmore   300  69      fancy  600 --- Christophers 300  27                  +- MadCows      300  24     -------------------------------------                          total  2100 300          k = number of times alternative is chosen     N = number of observations at each level          Iteration 0:   log likelihood = -541.93581       Iteration 1:   log likelihood = -517.95909  (backed up)     Iteration 2:   log likelihood = -511.99261  (backed up)     Iteration 3:   log likelihood =  -506.5559  (backed up)     Iteration 4:   log likelihood = -501.33127  (backed up)     Iteration 5:   log likelihood = -492.93121  (backed up)     Iteration 6:   log likelihood = -491.75884       Iteration 7:   log likelihood = -490.81911       Iteration 8:   log likelihood = -489.56195       Iteration 9:   log likelihood = -488.58891       Iteration 10:  log likelihood = -485.68122       Iteration 11:  log likelihood = -485.59216       Iteration 12:  log likelihood = -485.56202       Iteration 13:  log likelihood = -485.51992       Iteration 14:  log likelihood = -485.48062       Iteration 15:  log likelihood = -485.47384       Iteration 16:  log likelihood = -485.47333       Iteration 17:  log likelihood = -485.47331            RUM-consistent nested logit regression         Number of obs      =       2100     Case variable: family_id                       Number of cases    =        300          Alternative variable: restaurant               Alts per case: min =          7                                                                   avg =        7.0                                                                   max =          7                                                            Wald chi2(7)    =      46.71     Log likelihood = -485.47331                       Prob > chi2     =     0.0000          ------------------------------------------------------------------------------           chosen |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]     -------------+----------------------------------------------------------------     restaurant   |             cost |  -.1843847   .0933975    -1.97   0.048    -.3674404   -.0013289         distance |  -.3797474   .1003828    -3.78   0.000    -.5764941   -.1830007           rating |    .463694   .3264935     1.42   0.156    -.1762215     1.10361     ------------------------------------------------------------------------------     type equations     ------------------------------------------------------------------------------     fast         |           income |  -.0266038   .0117306    -2.27   0.023    -.0495952   -.0036123             kids |  -.0872584   .1385026    -0.63   0.529    -.3587184    .1842016     -------------+----------------------------------------------------------------     family       |           income |          0  (base)             kids |          0  (base)     -------------+----------------------------------------------------------------     fancy        |           income |   .0461827   .0090936     5.08   0.000     .0283595    .0640059             kids |  -.3959413   .1220356    -3.24   0.001    -.6351267   -.1567559     ------------------------------------------------------------------------------     dissimilarity parameters     ------------------------------------------------------------------------------     type         |        /fast_tau |   1.712878    1.48685                     -1.201295    4.627051      /family_tau |   2.505113   .9646351                       .614463    4.395763       /fancy_tau |   4.099844   2.810123                     -1.407896    9.607583     ------------------------------------------------------------------------------     LR test for IIA (tau = 1):           chi2(3) =     6.87   Prob > chi2 = 0.0762     ------------------------------------------------------------------------------      Here is the R code that I wrote to replicate that:               library(mlogit)     library(foreign)     dfNL = read.dta('http://www.stata-press.com/data/r12/restaurant.dta')          nlRest = mlogit(chosen ~ - 1 + cost + distance + rating | 0 + income + kids,            data = dfNL, un.nest.el = TRUE,            nests = list(fast = c('Freebirds', 'MamasPizza'),                          family = c('CafeEccell', 'LosNortenos', 'WingsNmore'),                         fancy = c('Christophers', 'MadCows')),            shape = 'long', choice = 'chosen', chid.var = 'family_id',            alt.var = 'restaurant')     summary(nlRest)      How can I specify that the parameters on "income" and "kids" should vary at the nest level rather than at the choice level? Should I impose constraints or is there a more natural way to do this?


title: When to use F-Test versus Chi-Square Test in Analysis of Deviance

text: I'm choosing a GAM, and when I try to do               anova(model.1, model.2, test="F")      It gives me a warning message saying that the F test is inappropriate for binary outcomes. Why is this? I get virtually the same p-value either way.


title: Effect of high number of support vectors

text: I am doing Rainfall-Runoff modeling. I have 4014 inputs and 4014 outputs. I am confused about the of support vectors. Suppose I have a model having 2000 support vectors and I have another model having 3900 support vectors. Both perform well, and the second model has slightly better performance. Which is better? What is the relation between number of support vectors and over fitting?


title: Panel regressions with an interaction term between a time dummy variable and a time invariant variable

text: I have estimated the coefficients of the following equation, using the fixed- effect model: $Y_{it}=\alpha _i+ \rho _t + \beta _1 X_{it}+\beta _2 C_i*D_t+\epsilon_{it}$ I have observations from 1980 to 2010 for $Y_{it}$ and $X_{it}$.I am interested in the interaction term. The dummy variable $D_t$ is equal to $1$ for years 2000 to 2010 (and $0$ otherwise) and is interacted with the continuous variable $C_i$ which is time-invariant (but varies over the units $i$). My question is: can $\beta _2$ be interpreted as the differentiated effect of $C$ on $Y$ during the period 2000-2010 compared to 1980-2000? Thanks!


title: Fitting a straight line to components of complex numbers

text: I have a strange problem that I'm not sure how to solve: I have complex data points in a time series. The amplitude of these complex numbers in the time series forms a straight line, which I have fit for. But now what I want to do is _subtract_ this fit from the data points. Except I have no idea how to do that since I have fit in amplitude space, and I have to subtract it from (real, imaginary) space. Is there any way to do this, or am I completely misguided?


title: Model selection in mixed-model context using lmer

text: Assume I have two factors A and B potentially predicting my outcome Y. Now I would like to test for fixed-effects using likelihood ratio test to find the best model.               fm1 <- lmer(Y~1+A*B+(1|subject))     fm2 <- lmer(Y~1+A+B+(1|subject))     fm3 <- lmer(Y~1+A+(1|subject))     fm4 <- lmer(Y~1+B+(1|subject))     fm5 <- lmer(Y~1+(1|subject))          anova(fm1, fm2, fm3, fm4, fm5)      However I would also like to test for different random-effect specifications.               # all possible random-effects specifications for fm1     fm6 <- lmer(Y~1+A*B+(1+A*B|subject))      fm7 <- lmer(Y~1+A*B+(1+A|subject)+(1+B|subject))     fm8 <- lmer(Y~1+A*B+(1+A|subject))     fm8 <- lmer(Y~1+A*B+(1+B|subject))      fm8 <- lmer(Y~1+A*B+(1|subject))       **Is it a valid approach to first detect which fixed-effects specification is most predictive [`anova(fm1, fm2, fm3, fm4, fm5)`] and then to try different random-effects specification with this model? Or would only the most comprehensive approach be valid in which I specify all possible combinations of fixed- and random-effects to compare all these model to find the best fit?**


title: Second survey after knowing how different strata respond?

text: We did customer survey in last fall from which we know how different stratas answer to survey questions. We did notice that amongst many strata having younger people who move often the non-response rate was much higher than amongst more established of our customers. Original survey was based on simple random sampling without replacement amongst strate and sample sizes strata-wise were based on proportional allocation. Should we use this information next time when we do customer survey? Should we allocate more for those strata which have more non-response. We could perhaps use some re-weighting to correct marginal distributions afterwards. Could this be reasonable approach? We have not done customer surveys earlier ourselves and previous survey was financed together with our industry association and we got only summary statistics. And survey was done by outside consultant firm. So we do not have any technical reports what they had actually done in that survey.


title: Generalized Linear Mixed Model

text: I have three independent variables (all binary) and two interaction variables that I would like to test for their effect on the dependent variable "Awareness_level". This dependent variable is also binary. I assigned participants randomly to one of 8 different conditions. In the assigned condition they were exposed to 1 of 8 advertisements (advertisements were varied based on different combinations of the 3 independent variables). All participants were asked the same questions to assess awareness: one question to assess recall and one question to assess recognition. I restructured the data in SPSS so that for each subject there are two rows of data. The first row (index variable awareness_type is 1) shows recall and the second row (index variable awareness_type is 2) shows recognition. I created the variable awareness_level as a target variable. See picture (might be more clear). ![enter image description here](http://i.stack.imgur.com/OXxOT.png) To add an intercept to the model that is different for every participant, I choose to fit a binomial logistic mixed model by running a Generalized linear mixed model with a binomial regression in SPSS. I think I should add all the independent variables (used to create the different conditions) + the interaction terms as fixed factors, while adding a random effect for subject ID. **Should I also add the index variable "awareness_type" as a independent variable?** I used this manual from SPSS as a basis (although it is for multinomial logistic mixed model) http://pic.dhe.ibm.com/infocenter/spssstat/v20r0m0/index.jsp?topic=%2Fcom.ibm.spss.statistics.cs%2Fglmm_cablesurvey_howto.htm They interpret the fit of the model by looking at the accuracy / classification table. But if I'm correct you should compare the accuracy rate to the proportional by chance accuracy rate right? **How do I calculate this one?** **And wouldn't it be better to assess the model fit based on the information criterion instead of accuracy rate? And which one is better than, Bayesian or Akaike corrected?** Would be great if you can help me!!