title: Help identifying time series model

text: I am trying to fit an arima model to data below. The best model seems to be ARMA(2,2) with the first AR and MA components zero to me. However, it does not quite satisfy me. I am looking for any ideas to lead me to a better model. Data plot, acf, pacf and fitted model results are below                 [1]  3.800731572  1.964861846 -0.330186755 -2.074280974  0.045544766  1.414932643       [7] -1.772733971 -0.269238758  2.201851709 -0.923099841 -2.022964462  1.244593286      [13]  0.138510229  0.207393784 -1.179208337 -2.600458137  1.597338790  2.672883672      [19]  0.021512541 -0.239735960 -0.241002298 -2.652076605 -1.006429172  4.069421128      [25]  1.567593649 -3.809127815 -1.098172055  0.461986823 -0.371368587  2.377791538      [31] -0.343117083 -1.627939880 -1.495075368  1.089181437  0.559939059 -1.628253076      [37]  2.110155105  1.016467745 -2.427617328  1.244370318  0.964439525 -1.102202024      [43]  0.165851643 -0.229573936 -1.624649814  0.295206971  0.068025881 -1.134450650      [49] -0.028201597  2.503538284  1.114222419 -2.833496083 -0.522543695  0.050113179      [55]  0.125369098  2.105797843  2.830923699 -4.856194214 -3.273881724  3.279800792      [61]  2.239616685 -0.908085837  0.139333522 -1.450735818 -0.803266180  2.190218597      [67]  1.173394051 -1.246471581  1.242838802  0.072678970 -2.275862413 -0.107232220      [73] -0.536815960 -0.194192692  4.280008694 -0.434172239 -4.423742767 -0.086640065      [79]  1.953905466  0.980978972  0.718926929 -2.118235181 -3.086168990  3.488865741      [85]  3.599936662 -3.739336414 -0.543348582  0.301020550 -1.055976766  2.578110734      [91]  1.595147532 -1.380082733 -1.554146884  1.611852730 -0.024437323 -0.695824440      [97]  1.212412282 -0.905360495 -1.731530168  0.175207750  2.228134579 -0.249534514     [103]  0.317061339  0.301498249 -2.915450651  0.700227471  4.178040621  0.300634992     [109] -2.775388196 -0.813531442  1.682700373  0.170627572 -0.304814504  0.593898664     [115]  0.705283783  0.247763384 -0.973869829  0.797997103 -0.210376025 -0.664768398     [121]  0.007792289 -0.123965987 -1.053204162  1.317299536  3.605380859 -0.207812137     [127] -2.753767618 -0.953270567  1.302774452  0.475310161  0.439760253 -0.519127147     [133] -2.713290201  0.528059101  1.590144928  0.290186322 -1.140784855 -0.290639801     [139]  0.451966832 -1.829276852  0.616356479  3.382008422 -0.372334534  0.586538986     [145]  1.065393107 -1.045478951 -1.297666091 -0.804821865 -0.894058890  0.945981384     [151]  1.865434286 -2.347928952 -1.145754072  4.593493727 -0.509447663 -1.593344586     [157]  0.161678089 -1.102974728  1.219269299  2.142194154 -2.169974430 -0.847790444     [163]  2.112060982 -1.994966429 -1.549014700  2.163255646 -0.794707401 -0.255578617     [169]  1.691674982 -0.901877878 -0.336901862  2.345808096 -1.869011631  0.415512418     [175]  1.417162468 -0.189171118  0.604161930 -2.354043433 -3.087017949  1.391578651     [181]  1.607888511 -2.601455766  1.653942000  1.040655043 -2.547348811 -0.811884895     [187]  1.291304418  1.346237124  0.646396579 -0.888786433 -1.238057892 -0.801669576     [193]  0.360090751  0.751050406  0.273321234 -0.919838506 -0.791943921  1.771953858     [199] -0.306807634 -1.033265249       ![enter image description here](http://i.stack.imgur.com/cY8ys.png) ![acf](http://i.stack.imgur.com/dqn6c.png) ![enter image description here](http://i.stack.imgur.com/Sgj5v.png)               ARIMA(2,0,2) with non-zero mean           Coefficients:           ar1      ar2  ma1      ma2  intercept             0  -0.3677    0  -0.6407    -0.0044     s.e.    0   0.0767    0   0.0536     0.0222          sigma^2 estimated as 1.366:  log likelihood=-316.07     AIC=638.14   AICc=638.57   BIC=657.93


title: Practical meaning of expected value, standard deviation & correlation

text: _We've got given annual results of two stock companies described with following values:_ **_Company X:** expected value $\mu_X=0.05$, standard deviation $\sigma_X=0.02$_ **_Company Y:** expected value $\mu_Y=0.07$, standard deviation $\sigma_Y=0.03$_ _We also know that a cross correlation between results of company X and Y is $\rho=-0.5, \rho=\frac{cov(X,Y)}{\sqrt{var(X)var(Y)}}$. How do we should divide **z** dollars between these companies to obtain:_ _a) maximum revenue (expected return on a portfolio)?_ _b) minimum risk (variance of a portfolio)?_ I don't even expect a full solution, just a brief explanation which would help me to understand a practical meaning of above values.


title: What are some guidelines for sample size when constructing a Bayesian Network?

text: Both general advice and advice specific to my case below is welcome. I'm trying to contruct a Bayesian network with 20-30 continuous variables. I have 15 locations sampled, with 14 time points each. The dataset will grow over time, but I have a feeling that there will be severe underfitting for quite a while. For my purposes, I have grouped the variables into 6 groups and limited which groups may be connected to other groups, including direction. I don't know how this affects the power of the algorithm finding the BN. A user said here that 30-100 minimum is preferred, although I'd feel better if there were a good source for this. I've seen "30 minimum" referenced in terms of the minimum amount of data you should have for a normal distribution (page 18, 2nd paragraph), but I'm not sure how to use this claim with respect to my data. In the future, when the sample size is larger, I also want to use time-lagged data to build a Bayesian network similar to a VAR model, but I would like to know if I will have the same issues in the variable selection process, or if there's a way to appropriately truncate the possibilities to make up for a smaller sample size.


title: Upper Bound on Mutual Information

text: I was advised to post this question here rather than in the math stack exchange. So here it is: I am interested in an upper bound on mutual information that I have been encountering frequently in the statistics and probability literature. I have yet to see the "purest" form of the inequality, so I will attempt to state it first, and then I will provide some literature that use leverage it. Ideally someone will be able to help me understand where this bound comes from and how it is derived. Here is the idea: Let $M = \left\{\theta_1,\theta_2,\ldots,\theta_m\right\}$ be a set of "models". Additionally, let $\phi : \mathbb{R}^{n\times m} \to \left\{1,2,\ldots,m\right\}$. Indeed, $\phi\left(X^n\right)\in \left\{1,2,\ldots,m\right\}$ where $X^n$ may be thought of as a data (or design) matrix of $n$ observations. In other words, we are trying to recover, through the estimator $\phi$ one of the $m$ models which best describes our data. I am interested in the following quantity: $I\left(\phi; X^n\right)$, where $I$ represents the mutual information. Then the following bounds $I$ from above: $$ I\left(\phi,X^n\right) \leq \frac{n}{m^2} \sum_{i=1}^m \sum_{j = 1}^m S\left(\mathbb{P}_{\theta_i} \left|\right| \mathbb{P}_{\theta_j}\right) $$ Here, $S\left(\cdot\left|\right|\cdot\right)$ denotes the symmetric KL- divergence. Further, $\mathbb{P}_{\theta_i}$ represents the underlying probability distribution for the model $\theta_i$. This bound appears in Exploiting the Limits of Structure Learning via Inherent Symmetry and Information-theoretic bounds on model selection for Gaussian Markov random ﬁelds. In the former, see equation (2). In the latter, see equation (18). For the latter paper though, I believe there is a typo as there seems to be a factor of $2n$ missing. (Please correct me if I am wrong.) Therefore, it may be best to consult the first paper.


title: Conceptual question on image pattern representation

text: I have a basic question regarding pattern learning, or pattern representation. Assume I have a complex pattern of this form, could you please provide me with some research directions or concepts that I can follow to learn how to represent (mathematically describe) these forms of patterns? in general the pattern does not have a closed contour nor it can be represented with analytical objects like boxes, circles etc. ![binary discrete pattern](http://i.stack.imgur.com/nuSqr.jpg)


title: Anonymizing survey data, plus an ethics question

text: This is a fairly simple question about raw data in Excel. I have a survey asking respondents to rate 17 questions on a rating scale from 1 to 5. The goal is for individuals who answer the survey to be anonymous. There are also demographic questions. For the sake of anonymity, we don't want someone to be able to match a single response to all of its demographic answers. Dissecting by all the demographic indicators could probably get you an individual person, and we don't want to report groups smaller than 10 people (arbitrary limit, management picked it). If a survey response is paired with a single demographic indicator, it usually creates groups of responses larger than 10. There are about 40,000 responses to this survey. So my first question is: Say I have 3 demographic questions. I copy all results of the 17 rating questions on to 3 tabs, and on each tab I include one of the 3 demographic indicators next to the survey data. Then I randomize the order of responses on all 3 tabs. I assume it's possible to de-anonymize this, so **What would be a better way to anonymize the data?** My concern is that there are about a trillion different possible permutations of survey answers, so someone could use the order of survey responses as a unique identifier. BUT: there are going to be repetitions, so they could never de-anonymize it completely. About 17% of the responses have at least 1 clone. Second question: **If we guarantee respondents anonymity, would we be maintaining that promise if we reported the data as proposed?** I recognize that there is a risk of the data being de-anonymized, but the data will be going mostly to people who are not very data savvy, and who honestly probably don't care enough to try to de-anonymize it. Thoughts?


title: Interpolating time series

text: what are best ways to interpolate time series? I have three data points(1980, 1990 and 2001) and I need to interpolate them. Using R na.approx doesn't seem to be what I need since the data I need to interpolate is population and it is very unlikely for it to move linearly.


title: How to plot correlation coefficients from multiple regressions while leaving out some of the variables from the plot?

text: Lets say that I have "y" that I want to model with linear regression. "x" and "z" are the things I'm interrested in showing folks, but I also have things that I want to adjust for, but not really show in my plot. Now, I would like to show my coefficients in a plot, but I would like to keep the things I adjusted for out of it. So perhaps this could be coronary artery calcification modeled as "CAC ~ SomeBloodStuff + Bloodpressure + BMI + Smoking_status". The BMI and Smoking_status would be something that I would want to take into account, but just note that I have adjusted for them. I gather this is how to do the model:               MyModel <- lm( CAC ~ SomeBloodstuff + BP + BMI + Smoking_status, data=MyData)     summary(MyModel)      coefficients(MyModel)     confint(MyModel, level=0.95) # Seems like a legit model.      The coefplot function in arm library gives just the sort of presentation that I want, but shows all the dimensions.               library(arm)     coefplot(MyModel)      How could I leave those few variables out of the plot while keeping them in the regression and get a plot that looks like what the coefplot produces?


title: Is there a version of Buhlmann-Straub credibility that uses an non-fixed $\theta_i$?

text: Everything I've read about Buhlmann-Straub credibility assumes a fixed $\theta_i$ (the unknown parameter that $X_{ij}$, the variable of interest, depends on). Does anyone know of a version where theta can vary, say with time?


title: How to handle edge cases in instant runoff voting?

text: I want to make forum software with instant runoff voting instead of polls, but I don't know what to do in edge cases. Should I allow voters to not weigh certain choices at all? If there are 5 choices and they only care about 3, and don't want weigh the remaining choices, should I allow that? If I allow people to not weigh choices, and all of their selected choices lose, does that mean that I would eliminate the voter completely? I wouldn't know what to do with them otherwise. Should I allow people to weigh multiple choices equally? If so I don't know how to handle that and need help. I think I have my head wrapped around instant runoff, but I've been completely wrong about statistics before.


title: Quasi maximum likelihood estimation versus pseudo MLE

text: If I'm not wrong both "quasi" and "pseudo" denote the same thing, namely the optimization under wrong distributional assumptions. Moreover I think that the terms are not restricted to the assumption of normality. Is there an experienced reader who can confirm this? Cheers!


title: How do I reduce the number of data points in a series?

text: I haven't studied statistics for over 10 years (and then just a basic course), so maybe my question is a bit hard to understand. Anyway, what I want to do is reduce the number of data points in a series. The x-axis is number of milliseconds since start of measurement and the y-axis is the reading for that point. Often there is thousands of data points, but I might only need a few hundreds. So my question is: How do I accurately reduce the number of data points? What is the process called? (So I can google it) Are there any prefered algorithms (I will implement it in C#) Hope you got some clues. Sorry for my lack of proper terminology. * * * Edit: More details comes here: The raw data I got is heart rate data, and in the form of number of milliseconds since last beat. Before plotting the data I calculate number of milliseconds from first sample, and the bpm (beats per minute) at each data point (60000/timesincelastbeat). I want to visualize the data, i.e. plot it in a line graph. I want to reduce the number of points in the graph from thousands to some hundreds. One option would be to calculate the average bpm for every second in the series, or maybe every 5 seconds or so. That would have been quite easy if I knew I would have at least one sample for each of those periods (seconds of 5-seconds-intervals).


title: which test should be used to detect the effect of one variable on another?

text: I have conducted a research on four independent groups of students; one experimental and one control group at two levels. The design of the study was like a pre-test and post-test. I want to detect the effects of acquiring how to use one variable (metacognitive strategies) on another variable (self- regulation). To obtain the index of the students' use of the first variable, a self-reported 10-item likert scale was distributed among the students of the groups. To assess self-regulation of the students a self-reported 32-item likert scale was also distributed. Then an intervention program related to the first questionnaire was employed for the experimental groups. then the same questionnaires were filled by the students as a post-test. Comparing means, the results showed a significant difference in the means of the post-tests. Now, I want to use a test to detect the effects of the first questionnaire (or program) on the second, both of which are likert type with pre- and post treatment data. how can I show the effects of this intervention program? is correlation a correct solution?


title: Factor analysis results

text: I carried out factor analysis on 31 5 likert-scale questions which represent my 6 constructs. The results showed 6 factors, consistent with my hypothesis but the analysis grouped my questions differently from what I expected. For example when I wrote the survey I dedicated the first 5 questions as measuring one variable, the next five to measure another, etc. The results of the analysis grouped questions 1.2.7.8.15.16 in one variable and other questions in other group, so what happened is that the factor analysis mixed my classification. What should I do in this case? Should I continue my analysis and test the relationship between my variables based on my own classification or what ?


title: Uniform minimum variance unbiased estimator

text: Let $X_1, X_2, ..., X_n$ be an iid random sample from a Poisson$(\lambda)$ distribution: a) find the UMVUE of $\theta$ = $\lambda^k$ for $k > 0$ a known integer b) find the UMVUE of $\tau = e^{-\lambda} = \mathbb{P}(X = 0)$. **Edit:** I have found that $T = \sum X_i$ is a complete sufficient statistic. It has a Poisson$(n\lambda)$ distribution. Now to find the UMVUE, we need a function $g(t)$ such that $E[g(T)] = \theta$, then $g(t)$ will be UMVUE. For, that $E[g(T)] = \theta$ implies $$e^{-n\lambda}\sum_{t=0}^\infty g(t) \frac{(n\lambda)^t}{t!} = \theta =\lambda^k.$$ Now I don't know how to proceed here. If it were the continuous case, I could use the differential regularity condition. The same thing in b) too.


title: Kernel selection intuition

text: Which problems are best solved using which kernels and why? Can you give a simple toy problem that isn't linearly separable in input space but is linearly separable in feature space, using an RBF kernel? A polynomial kernel?


title: Specifying the correlation structure of an unevely spaced time series in GEE with geepack

text: I have counts of plants from different sites over a number of years. In each census year, all sites were surveyed, but the gaps between census years vary (between 1 and 4 years between consecutive surveys). I would like to use a GEE to account for the correlation structure of the data (but really, my purpose is just in estimating the long-term population trend, and some estimate of its uncertainty). I think I can specify a model like this:               library(geepack)     model1<-geeglm(Count~Year,id=Site,family=poisson,data=mydata,corstr="ar1",waves=Year)      But I am not sure if this is quite right--I think the correlation parameter estimate it gives is the correlation between any two consecutive counts, i.e. without accounting for the time between counts. My question relates to a previous one below, but I could not comment on it. The waves argument makes no differences to my results (the data is already ordered by site, and time within site). Can anyone point me in the right direction? Or comment on how much I should worry about this? How to analyze GEE with unevenly spaced observations?


title: Likelihood of censored data

text: Let $X_1,X_2,\ldots, X_{n_1}$ be IID with PDF $f(x-\theta) $, for $-\infty<x<\infty$ and $-\infty<\theta<\infty$. Denote the CDF of $X_i$ by $F(x-\theta)$. Let $Z_1,Z_2, \ldots, Z_{n_2}$ denote the censored observations. For these observations we only know that $Z_j >a $ for some $a$ that is known and that the $Z_j$s are independent of the $X_i$s. Then the "observed" and "complete" , that is for both the observed and the censored data, likelihoods are given by: $$L(\theta| \mathbf{x} )= [1-F(a-\theta)]^{n_2} \prod_{i=1}^{n_1} f(x_i -\theta) $$ $$ L( \theta |\mathbf{x,z})=\prod^{n_1}_{i=1}f(x_i -\theta) \prod_{i=1}^{n_2}f(z_i-\theta) $$ My question is why do the likelihoods take that form? I believe for the $Z_j$s we have a left censoring as they only take greater than $a$ values, as opposed to the $X_is$ that are not restricted in that way. We also know that they are mutually independent. Could it be that the "observed" likelihood is obtained after using that independence? We have $n_2$ observations greater than $a$, times the joint density of the $X_i$s, right?. A little more explanation on these two equations would go a long way. Thank you in advance.


title: Understanding statistical significance in variance and correlation?

text: I have little background in statistics. For a dozen items, I have vectors of ratings for these items. I have calculated the variance, and Pearson correlation for these items (using numpy). Now I would like to compare them. My adviser asked me to see if the differences in variance and correlation were "statistically significant" using the 95% rule. How can I do this? How does this statistical significance relate to scipy's p-value which is returned with Pearson correlation (it makes no mention of "95%")? > The p-value roughly indicates the probability of an uncorrelated system > producing >datasets that have a Pearson correlation at least as extreme as > the one computed from >these datasets. The p-values are not entirely > reliable but are probably reasonable for >a datasets larger than 500 or so.


title: Optimal to sub-optimal sampling patterns

text: In the compressed sensing (CS) framework an incoherence property $\mu$ is used to quantify the correlation between a sensing matrix $\Phi$ and a representation matrix $\Psi$. This incoherence relationship is given by ![enter image description here](http://i.stack.imgur.com/gCcHJ.jpg) where $n$ is the number of elements in the matrix (signal). As far as I am able to tell, CS assumes that $\Phi$ is a matrix designed to be as orthogonal as possible to the representation basis, hence both are usually known prior to the sampling step and $\mu$ can be calculated; $\Psi$ may be chosen as a basis which best represents the data in a sparse/compressed sense (e.g. Fourier basis, Wavelet basis, some fixed basis). I have a set of data samples collected from an image. The samples are lined up along parallel diagonal straight lines of various lengths that go across the image. Is it possible to deduce a sensing matrix $\Phi$ from this pattern so as to define an incoherence? My goal is to try to measure a "distance" by which the sampling pattern I have been given is sub-optimal compared to an optimal sensing matrix. In this way I want to be able to show how one sampling pattern might be better than another. If that's not possible, is the next best thing to use a Restricted Isometry Property (RIP) $\delta_k$ to define whether $\Phi$ meets the RIP? ![enter image description here](http://i.stack.imgur.com/8mMtU.jpg) where $y$ = $\Phi v$ = samples, and $v$ is the exhaustive dataset (i.e. all of the data, which I have)


title: Inconsistent Z-scores from multiple comparisons of factors in glmmadmb models

text: First time poster, long time desperate reader. This place has been so helpful that I can usually find that someone has had the same problem as myself and had it answered here. However, I haven't found an answer to my problem below in a way that I understand. I am having trouble finding consistent z-scores with `glmmadmb` when used in conjunction with `glht`. The z-scores attached to the reported p-values of a 6-level categorical variable match up with the p-values pinpointed on a z-table using the z-scores from the output of a glmmadmb model. But after a multiple comparisons, they do not. When testing a simpler 2-level factor the z-scores seem to be accurate for the original model and a multiple comparisons. Examples below. My query is three-fold:   1. I'm confused. Why is this happening?   2. Which test-statistic should I be reporting when using these comparisons?   3. How do I find the values of those test statistics? I have been using glmmadmb (glmmADMB) to model bee visit duration (minutes) of bees to flowers with a categorical fixed factor (hour of morning: 6AM,7AM,8AM,9AM,10AM,11AM) and random factor (site). I chose to use a gamma distribution because the response variable is right-skewed, non-zero, and continuous. Code and output below:               sb<- glmmadmb(min~hour + (1|site), family = "gamma", data=sbsub)     summary(sb)           Call:     glmmadmb(formula = min ~ hour + (1 | site), data = sbsub, family = "gamma")          AIC: 134           Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)    0.117      0.291    0.40  0.68708         hour7         -0.794      0.322   -2.47  0.01363 *       hour8         -1.112      0.315   -3.53  0.00042 ***     hour9         -0.642      0.323   -1.99  0.04660 *       hour10        -1.178      0.311   -3.79  0.00015 ***     hour11         0.350      0.329    1.07  0.28646         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Number of observations: total=93, site=11      Random effect variance(s):     Group=site                 Variance StdDev     (Intercept)   0.3774 0.6143          Gamma shape parameter: 1.4184 (std. err.: 0.19914)          Log-likelihood: -58.9969       The p-values in the output can be found in a z-table after using the listed z-score from the outputs and multiplying the z-table value by 2...But, I did not just want to see how bee visit duration at each hour compares to the base level of 6AM. I wanted to know how visit duration compares across all combinations of hours. So I used `glht` and `cld` from the `multcomp` package and got the output below. Looks sexy, and informative. I could plop those letters right on a bar chart.               sb2<-glht(sb, linfct = mcp(hour = "Tukey"))     summary(sb2)               Simultaneous Tests for General Linear Hypotheses          Multiple Comparisons of Means: Tukey Contrasts               Fit: glmmadmb(formula = min ~ hour + (1 | site), data = sbsub, family = "gamma")          Linear Hypotheses:                  Estimate Std. Error z value Pr(>|z|)        7 - 6 == 0   -0.79358    0.32168  -2.467   0.1206        8 - 6 == 0   -1.11241    0.31528  -3.528    <0.01 **     9 - 6 == 0   -0.64218    0.32272  -1.990   0.3215        10 - 6 == 0  -1.17755    0.31065  -3.791    <0.01 **     11 - 6 == 0   0.35035    0.32868   1.066   0.8793        8 - 7 == 0   -0.31883    0.44769  -0.712   0.9767        9 - 7 == 0    0.15140    0.45244   0.335   0.9993        10 - 7 == 0  -0.38398    0.44428  -0.864   0.9469        11 - 7 == 0   1.14393    0.46486   2.461   0.1232        9 - 8 == 0    0.47023    0.45019   1.044   0.8881        10 - 8 == 0  -0.06515    0.44124  -0.148   1.0000        11 - 8 == 0   1.46276    0.45547   3.212   0.0151 *      10 - 9 == 0  -0.53537    0.43020  -1.244   0.7919        11 - 9 == 0   0.99253    0.46182   2.149   0.2399        11 - 10 == 0  1.52790    0.46197   3.307   0.0110 *      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1     (Adjusted p values reported -- single-step method)          cld(sb2)        6    7    8    9   10   11        "b" "ab"  "a" "ab"  "a"  "b"       But, these z-scores don't seem to mean anything. Looking them up in the z-table doesn't reveal the same p-value/2 in the output. So, reporting them seems moot. However, I feel I should report some kind of test-statistic to accompany the p-values? This also seems to happen when using `glmer` with a binomial distribution modeling predation and parasitism. In that case, the z-scores from a multiple comparisons between a 3-level categorical variable and a 2-level categorical variable do not have accurate z-scores. I tried it with another smaller pollination model testing visit duration of bees in male and female flowers.               sb<- glmmadmb(min~sex + (1|site), family = "gamma", data=sbsub)     summary(sb)          Call:     glmmadmb(formula = min ~ sex + (1 | site), data = sbsub, family = "gamma")          AIC: 154.4           Coefficients:                 Estimate Std. Error z value Pr(>|z|)       (Intercept)   -0.154      0.257   -0.60    0.550       sex1          -0.367      0.215   -1.71    0.087 .     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Number of observations: total=93, site=11      Random effect variance(s):     Group=site                 Variance StdDev     (Intercept)   0.4517 0.6721          Gamma shape parameter: 1.102 (std. err.: 0.15293)          Log-likelihood: -73.1965                sb2<-glht(sb, linfct = mcp(sex = "Tukey"))     summary(sb2)          Fit: glmmadmb(formula = min ~ sex + (1 | site), data = sbsub, family = "gamma")          Linear Hypotheses:                Estimate Std. Error z value Pr(>|z|)       1 - 0 == 0  -0.3669     0.2147  -1.709   0.0875 .     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1     (Adjusted p values reported -- single-step method)          cld(sb2)        0    1           "a"  "a"        The p-values/2 in the output of both the original model and the multiple comparison can be found after using the listed z-score from the outputs to find a value in the z-table. Not sure what's happening here. Any thoughts? Much appreciated, Nava


title: Modeling non-exclusive choices

text: I am trying to model votes by individuals on a number of bills over time using bill-specific characteristics as explanatory variables. Imagine that the same bill is offered to them albeit with minor changes from year to year. **I want to determine what sort of changes in text of the bill will make individuals vote yes, if they had voted no before. (i.e. what is it about the bill features that makes individual switch from no(0) to yes (1) for a certain bill).** My dependent variable is 0-1 whether a person voted yes at a given time. My main explanatory variables are characteristics of each bills (bills are scored on a number of dimensions). Ideally, I would want to model both bill-specific and individual-specific independent variables (perhaps via interaction?). **I am really confused about what model i should use. A simple individual- fixed effects logit is not a perfect model. I was also thinking of conditional logit, but the choices are non-exclusive and at any given time an individual is choosing whether to vote yes on only one bill. I would really appreciate any model suggestions.** My dataframe looks like this (but can be reworked):               I year   bill1 bill2 bill3 bill4 billfeature1 billfeature2  billfeature3     1   1990  0     1    1       0         5                3              1       1   1991  0     0    1       1         6                3              3       1   1992  0     0    0       0         7                4              6       1   1993  0     0    1       1         7                5              7       1   1994  1     1    1       1         8                6              8       2   1990  0     1    1       1         5                3              1           2   1991  0     0    1       0         6                3              3             2   1992  0     0    0       0         7                4              6       2   1993  0     0    1       1         7                5              7     2   1994  0     0    1       1         8                6              8


title: Inference with Gaussian Random Variable

text: Let $X = N(0,\frac{1}{\alpha})$, $Y = 2X + 8 + N_{y}$, and $N_{y}$ be a noise $N_{y} = N(0,1)$. Then, $P(y|x) = \frac{1}{\sqrt{2\pi}}exp\{ -\frac{1}{2}(y - 2x - 8)^{2} \}$ and $P(x) = \sqrt{\frac{\alpha}{2\pi}}exp\{-\frac{\alpha x^{2}}{2}\} $. The mean vector is: $$\mathbf{\mu} = \left( \begin{array}{c} \mu_{x}\\ \mu_{y}\end{array} \right)= \left( \begin{array}{c} 0\\ 8\end{array} \right).$$ The question is how to calculate the variance of Y. I know that the correct answer is $$\frac{4}{\alpha} + 1, $$ but don't know how to get from $$var(Y) = E[(Y-\mu_{y})^{2}] = E[(2X+N_{y})^{2}] $$ to $$\frac{4}{\alpha} + 1. $$ Can anybody help? UPDATE: Thank you All for answers


title: Software for easy-yet-robust data exploration

text: In my attempts to fight spreadsheet mayhem, I am often evangelical in pushing for more robust tools such as true statistics software (R, Stata, and the like). Recently, I was challenged on this view by someone who stated flat out that they simply will not learn to program. I would like to provide them with data analysis tools that require no programming (but ideally which would extend to programming if they decide to dip a toe into the water later). What packages are out there for data exploration that I can recommend with a straight face?


title: How to do a linear model in R?

text: The objective of this research was to investigate the long-term effects of irrigation with treated waste water on some chemical soil properties. The investigation was carried out by comparison of soil properties in two different fields: one irrigated with the effluent from Parkan Waste water Treatment Plant over a period of six years, and the other one irrigated with water over the same period of time. Soil samples were taken from different depths of 0-15, 15-30, 30-60, 60-100 and 100-150 cm in both fields, and analyzed for various chemical properties. For visual summaries, I am going to plot depths of soil and element (with line and bar plots). However, I also need to fit a linear model with one categorical and one continuous predictor. How can this be done in R? Thanks for your answer!


title: How to statistically compare two algorithms across three datasets in feature selection and classification?

text: **Problem background:** As part of my research, I have written two algorithms that can select a set of features from a data set (gene expression data from cancer patients). These features are then tested to see how well they can classify an unseen sample as either cancer or non-cancer. For each run of the algorithm, a solution (a set of features) is generated and tested on Z unseen samples. Percentage accuracy of the solution is calculated like this: `(correct classifications / Z) * 100`. I have two algorithms: algorithm X & algorithm Y I have three separate (different cancers) data sets: data set A, data set B & data set C. These data sets are very different from each other. They don't have the same number of samples or same number of measurements (features) per sample. I have run each algorithm 10 times on each data set. So, algorithm X has 10 results from data set A, 10 from data set B and 10 from data set C. Overall, algorithm X has 30 results. **My problem:** I would like to see if algorithm X's combined performance across all three data sets is statistically significantly different from algorithm Y's combined performance. Is it possible for me to combine results for algorithm X from each data set into a single set of results? This way, I would have 30 standardized results for algorithm X and 30 standardized results for algorithm Y. I can then use the t-test to see if there is a significant difference between the two methods. **Edit -** These are Evolutionary Algorithms, so they return a slightly different solution each time they are run. However, if there's a feature in a sample that when present can strongly classify the sample as either being cancer or non-cancer, then that feature will be selected almost every time the algorithm is run. I get **_slightly_** different results for each of the 10 runs due to the following reasons:   * These algorithms are randomly seeded.   * I use repeated random sub-sampling validation (10 repeats).   * The datasets that I use (DNA microarray and Proteomics) are very difficult to work with in the sense that there are many local optima the algorithm can get stuck in.   * There are lots of inter-feature and inter-subset interactions that I would like to detect.   * I train 50 chromosomes and they are not restricted to any particular length. They are free to grow and shrink (although selection pressure guides them towards shorter lengths). This also introduces some variation to the final result. Having said, the algorithm almost always selects a particular subset of features! Here's a sample of my results (only 4 runs out of 10 for each algorithm is shown here):               Dataset/run     Algorithm X     Algorithm Y     A 1             90.91           90.91     A 2             90.91           95.45     A 3             90.91           90.91     A 4             90.91           90.91          B 1             100             100     B 2             100             100     B 3             95.65           100     B 4             95.65           86.96          C 1             90.32           87.10     C 2             70.97           80.65     C 3             96.77           83.87     C 4             80.65           83.87      As you can see, I've put together results for two algorithms from three datasets. I can do Kruskal-Wallis test on this data but will it be valid? I ask this because:   * I'm not sure accuracies in different data sets are commensurable. If they are not, then putting them together like I've done would be meaningless and any statistical test done on them would also be meaningless.   * When you put accuracies together like this, the overall result is susceptible to outliers. One algorithm's excellent performance on one dataset may mask it's average performance on another dataset. I cannot use t-test in this case either, this is because:   * Commensurability - the t-test only makes sense when the differences over the data sets are commensurate.   * t-test requires that the differences between the two algorithms compared are distributed normally, there's no way (at least that I'm aware of) to guarantee this condition in my case.   * t-test is affected by outliers which skew the test statistics and decrease the test’s power by increasing the estimated standard error. What do you think? How can I do a comparison between algorithm X & Y in this case?


title: How to assign ranks given the distribution above mean?

text: Say I have the following data,   * % students below average in math: S1   * % of students at proficient level in math: S2   * % of students at above proficient level: S3 S1 + S2 + S3 = 100% If I have this distribution across multiple classes and multiple years, how can I do a comparitive assesment ? A naive score would be a weighted mean like S1 + 2S2 + 3S3 (so a bias towards strong performers). Anybody has pointers on better estimation ?


title: R Language what is diffrence between rnorm and runif

text: R Language what is diffrence between functions rnorm and runif


title: What is the difference between "two-ways" and "individual" effects in GMM estimation for panel analysis?

text: I am starting to use GMM estimation of dynamic panels with fixed-effects in R, but couldn't find proper information on the differences between "two-ways" and "individual" effects. Both can be easily estimated in R, but I couldn't get the right conceptual difference. Any help or indication here will be most welcome.


title: Unexpected error from running predict variance after mgarch dcc

text: I am getting an unexpected error from running predict variance after mgarch dcc command. I am using Stata 13 for Windows. I tried to find the dynamic correlation between two time series by using mgarch dcc step by step.               arch var1, arch(1) garch(1)     predict res1, residual     predict h1, var     gen std1=sqrt(h1)     gen stdres1=res1/std1       Do the same for var2,get stdres2. Then I typed:               mgarch dcc ( stdres1 stdres2)     predict H*, variance      then Stata sent error message:               MGARCH_predict_DCC(): 3301 subscript invalid     <istmt>: - function returned error*       Normally when I run predict H*, variance after running the command: mgarch dcc (var1, var2), arch(1) garch(1), there is no error. Do you know where I were wrong? How can I fix it? I already searched everywhere but still can not find the answer or even someone has similar problem.


title: Bayesian mean square error

text: Given a i.i.d sample $X_{1},..,X_{n}$ of bernoulli random variables test 2 hypotheses $H_{0}:p=2/3$ and $H_{1}:p=1/3$. Bayesian prior is $\pi(2/3)=1/3$ and $\pi(1/3)=2/3$. Find the bayesian criterion for acceptng $H_{0}$, find the bayesian mean square error for the test and for $n=8$ compute this mean square error using normal approximation I have found the bayesian criterion for acceptance as $\sum_{i=1}^{n}x_{i}{\geq}\frac{n+1-log_{2}(\alpha^{-1}-1)}{2}$. where $\alpha$ is a value is chosen prior to the test. How do you do the other two parts? Thanks


title: The union, intersection and complement of events

text: On the Probability chapter of a 1995 mathematical statistical book I am reviewing I have found the following exercise: > Let A and B be arbitrary events. Let C be the event that either A occurs or > B occurs, but not both. Express C in terms of A and B using any of the basic > operations of union, intersection and complement. Now, the book suggested answer is to describe the entire sample space as: $$\Omega=(A\cap B)^{C}\cap(A\cup B)$$ I think the correct answer is: $$\Omega=(A\cap B)^{C}\cup(A\cap B)$$ and $$C=(A\cap B)^{C}$$ Where is the error that I have made?


title: Calculating Log Prob. of Dirichlet distribution in High Dimensions

text: I'm interested in calculating the log probability of data drawn from a Dirichlet distribution. In particular, I'm interested in calculations that are stable in high dimensions, perhaps 1000 dimensions or more. I've tried to use the "dirichlet_logProb.m" function from Tom Minka's fastfit toolbox: http://research.microsoft.com/en- us/um/people/minka/software/fastfit/. However, in high dimensions, this seems to give unexpected large _positive_ density values to _every_ data vector I provide. Certainly, when a distribution is sharply peaked, we can have arbitrarily large (greater than one) density values at some points in the support. For example, think of a zero-mean Normal distribution with very small variance... it will give very large density (>>1) to x close to zero, but vanishingly small density (<<1) to x more than a few standard deviations away. The same should be true of _any_ distribution: we should always be able to find points that have large density and small density, since the density needs to integrate to one over the support. However, in high dimensions I cannot get my Dirichlet log probability calculations to work out this way... every single point I try has density greater than one. Here's a sketch of my experiment: I consider a symmetric Dirichlet with concentration parameter <1. This should prefer data that are sparse, where only some dimensions have significant probability mass, and all others are near zero. I then look at three possible draws, which I illustrate here using K=4 dimensions.   1. Extremely sparse: x = [1 0 0 0]   2. Halfway uniform: x = [0.5 0.5 0 0]   3. Uniform: x = [0.25 0.25 0.25 0.25] Given concentration parameter lambda < 1, we should expect 1 to have much larger density than 2, and 2 to have much larger density than 3. Since the uniform is basically the opposite of sparsity, we should also expect this outcome to have negligible (less than one) density... which corresponds to negative log probability. Here is a result of the computed log probabilities for these three outcomes, across many values of the dimension K, with concentration=0.01. Remember, these are _log_ probabilities, so positive values imply density > 1, and negative values density < 1.                       [1 0...0]  [unif 0..0] [ unif ]     K=   4   9.18e+01   5.59e+01  -9.71e+00     K=  10   2.77e+02   1.43e+02  -2.09e+01     K=  50   1.52e+03   7.85e+02  -3.58e+01     K= 100   3.07e+03   1.64e+03  -4.04e+00     K= 500   1.55e+04   8.98e+03   7.80e+02     K=1000   3.11e+04   1.87e+04   2.25e+03      You can reproduce this table with this code: http://gist.github.com/3699842 The problem, of course, is that with K=500 and K=1000 _all_ of the outcomes 1,2,3 have very positive log probability (density >> 1). This seems troublesome to me, since I'd expect the most unlikely outcome (the uniform) to have density less than one... am I mistaken somehow? Can anybody suggest an answer? Is this a numerical issue? Or do I misunderstand something about the Dirichlet? My prime suspect is the fact that the support region is a bit funny (a simplex, rather than the reals).


title: Forecasting discrete non normal data

text: Today I was thinking about how one would treat what I'm guessing must come up constantly in statistics so apologies if it's trivial. Some statistics such are easy to forecast, such as length, thanks to their normalish properties. But what if I'm interested in forecasting something like the number of children someone will birth? This seems to be a rather weird case to me. It's not only skewed, but it also has a lower bound of zero. Also, even if we're disregarding the time of pregnancy, the probability of birthing a child in a given year isn't independent, i.e. $P(X>0) \neq P(X>1|X>0)$, so it shouldn't be Poisson either. How would one treat a phenomenon like this?


title: Rigorous Bell-CHSH?

text: The usual derivations by physicists of the Bell-CHSH inequalities as a limit imposed on any possible local hidden variables description of quantum systems seems so sloppy, that I'm tempted to ask you statisticians if there is a rigorous mathematical derivation of these results. Also, what would be a proper statistical analysis of the corresponding experiments (such as the one by Aspect et al)?


title: Probability of class membership given univariate normal distribution

text: Assuming a class is well described with a normal distribution of _u_ and _s_ , is it reasonable to calculate the probability of membership as: $Pr(x)=Pr(|x-u|)=2.0*\text{cdf}(|x-u|,s)$? I've briefly researched unary classification mechanisms, but I believe them to be out of scope since the probabilities computed above are expected to be part of a Bayesian inference system. Suggestions? If the above is deemed reasonable, does this method apply to bivariate and trivariate probability of membership if independence is assumed (i.e. covariance matrix is diagonal)?


title: Comparative Sentence Classification

text: I've been having a difficult time looking for a solution to identifying the comparative sentence. I like to extract the first noun (NP1) or second noun (NP2) from the given sentence fragment with comparaed word. First, I want to extract the comparative sentence, then I want to extract the entity(NP1).


title: Heteroskedasticity in a regression model

text: $\newcommand{\var}{\mathrm{Var}}$In my textbook I am studying about heteroskedasticity in a linear regression model, and for doing the hypothesis test, it says our hypothesis is: $H_0: \var(u|x_1,x_2,...,x_k)=\sigma ^2$, however than it says: > Because we are assuming that $u$ has a zero conditional expectation, > $\var(u|\mathbf{x})=E(u^2|\mathbf{x})$ (where $\mathbf{x}$ is all the > regressors $x_1,x_2,x_3,...,x_n$), and so the null hypothessis of > homoskedasticity is equivalent to: $$H_0: > E(u^2|x_1,x_2,...,x_k)=E(u^2)=\sigma^2$$ Why is this true? I dont understand why we use $u^2$ and then regress the independent variables on it.


title: Expressing observations in regression

text: Suppose we are using logistic regression on a binary outcome variable $Y$ with two continuous predictors $X_1$ and $X_2$. Suppose the observation $Y = 1, X_1 = 10, X_2 = 15$ occurs 1000 times. In SAS, how would you account for this frequency in the regression without having to manually type the observation 1000 times? **Edit.** I figured it out.


title: Fuzzy Record Linkage of Spatial Datasets

text: I have two datasets describing real-estate properties Dataset 1 describes building characteristics; it includes the location of the entrance to the building along with building descriptions and demographics;               Longitude     Latitude          Quality Index       Age of building (years) (.....hundreds)     13.3166960844 52.4720710276     1                    5     13.3724798666 52.5319328086     4                    11     13.3810401552 52.53146403       5                    12     (.... millions)      Dataset 2 describes real estate transactions; it includes locations, transaction price and date               Longitude     Latitude          Transaction Date     Transaction Price     13.3943043389 52.5342395928     01.05.12             100000     13.3946001820 52.53146403       10.03.10             324500        13.3927096776 52.532384714      04.12.05             400000     (..... millions)      I would like to join the two datasets, so each record in dataset 2 has the corresponding data from dataset 1, but there are several issues -   1. The location in dataset 1 is very accurate, but the one in dataset 2 is less accurate and may be 10-1000 meters off the actual location   2. The properties in dataset 1, along with those in dataset 2, should be able to explain at least 90% of the variance in the transaction price, if you could only match the records in the two datasets I heard there are good record linkage algorithms out there - but have no experience with them. Otherwise, I am quite experienced with (geo-) regression techniques. But frankly, I have no idea where I should even start.


title: Absolute convergence requirement for Expectation

text: In calculating the expectation of a discrete random variable $X$, we not only require that $\Sigma x_iP(X=x_i)$ converges, but also converges absolutely. I understand this requirement as probably stemming from the fact that a rearrangement of non-absolutely-convergent countably-infinite sum can have a different sum. I was wondering if a similar requirement exists for purely continuous random variable $X$, when we compute the expectation using the Riemann integral $\int_{-\infty}^\infty x f(x) dx$. [I have a similar question for computing expectation using Stieltjes integrals - is there some sort of "absolute convergence" requirement?] I understand that the most general definition of expectation involves Lebesgue integrals, but I am not very familiar with Lebesgue theory, so to be concrete (if you intend to reply via Lebesgue theory): In the special case of a purely continuous random variable, does the Lebesgue integral, when reduced to Riemann integral, have any form of "absolute convergence" requirement, or it is automatically satisfied in some sense? What about the case of Stieltjes integral? How does the "absolute convergence" requirement manifest itself for the discrete case?


title: Correlation between scale and categorical variable

text: I am trying to find the correlation or association between a categorical variable - classes from 1 -4 from a GIS ( _geographic information system_ ) classification - and a continuous variable with values between 0 and 9. What would be the best test to use for this?


title: Odds at least 1 person is born in January?

text: There are 100 people in a room. What are the odds at least 1 person is born in January? What is the best way to calculate this? I used Binomial Distribution. $E(x) = np = 100 \frac{1}{12} = 8.3$ $std.dev. = \sqrt{npq} = 2.76$ $Z(\text{1 person}) = (1 - 8.3) / 2.76 = -2.64$ $p = P(Z<-2.64) = .004$ or $.4%$ odds that less than 1 person is born in Jan Can a confidence interval calc apply here? If so, how? (With 95% certainty, the number of people born in Jan is b/w x & y)


title: Kernel Bandwidth: Scott's vs. Silverman's rules

text: Could anyone explain in plain English what the difference is between Scott's and Silverman's rules of thumb for bandwidth selection? Specifically, **when** is one better than the other? Is it related to the underlying distribution? Number of samples? P.S. I am referring to the code in SciPy.


title: How to calculate average time in service

text: We want to work out a value for "Average time in service" for our users. However, as we have many subscribers still active, how can we do that? If we just look at average time of deactivated subscribers, that won't show the picture. Is there some "correct" way to do this? We know our churn rate by month, so maybe we need to work out from those still on the service, when they're likely to leave? and use that?


I'm fitting a logistic regression to a simple two-feature dataset. The feature values range from about 1 to 100. When I scale the features (using scikit- learn StandardScaler), the in-sample performance improves significantly, but the out of sample performance drops. Why would feature scaling cause overfitting?


title: Is a random sample of a Poisson distribution also Poisson distributed?

text: Car analogy: Assume the traffic (number of cars per hour) on a road has a Poisson distribution, and the time between cars has the matching exponential distribution. If the chance of each car being red is independent from both time and the color of other cars, will the number of red cars per hour also have a Poisson distribution? I strongly suspect so. Furthermore I expect E(red cars)=p(red) * E(cars), and because of the Poisson distribution σ(red cars) = p(red) * σ(cars). But how would I (dis)prove this?


title: Why is AUC higher for a classifier that is less accurate than for one that is more accurate?

text: I have two classifiers   * A: naive Bayesian network   * B: tree (singly-connected) Bayesian network In terms of accuracy and other measures, A performs comparatively worse than B. However, when I use the R packages ROCR and AUC to perform ROC analysis, it turns out that the AUC for A is higher than the AUC for B. Why is this happening? The true positive (tp), false positive (fp), false negative (fn), true negative (tn), sensitivity (sen), specificity (spec), positive predictive value (ppv), negative predictive value (npv), and accuracy (acc) for A and B are as follows.               +------+---------+---------+     |      |    A    |    B    |     +------+---------+---------+     | tp   | 3601    | 769     |     | fp   | 0       | 0       |     | fn   | 6569    | 5918    |     | tn   | 15655   | 19138   |     | sens | 0.35408 | 0.11500 |     | spec | 1.00000 | 1.00000 |     | ppv  | 1.00000 | 1.00000 |     | npv  | 0.70442 | 0.76381 |     | acc  | 0.74563 | 0.77084 |     +------+---------+---------+      With the exception of sens and ties (spec and ppv) on the marginals (excluding tp, fn, fn, and tn), B seems to perform better than A. When I compute the AUC for sens (y-axis) vs 1-spec (x-axis)               aucroc <- auc(roc(data$prediction,data$labels));      here is the AUC comparison.               +----------------+---------+---------+     |                |    A    |    B    |     +----------------+---------+---------+     | sens vs 1-spec | 0.77540 | 0.64590 |     | sens vs spec   | 0.70770 | 0.61000 |     +----------------+---------+---------+      So here are my questions:   * Why is the AUC for A better than B, when B "seems" to outperform A with respect to accuracy?    * So, how do I really judge/compare the classification performances of A and B? i mean, do i use the AUC value? do i use the acc value? and why?    * Furthermore, when I apply proper scoring rules to A and B, B outperforms A in terms of log loss, quadratic loss, and spherical loss (p < 0.001). How do these weigh in on judging classification performance with respect to AUC?   * The ROC graph for A looks very smooth (it is a curved arc), but the ROC graph for B looks like a set of connected lines. Why is this? As requested, here are the plots for model A. ![model A naive bayes net](http://i.stack.imgur.com/vuPGZ.jpg) Here are the plots for model B. ![model B regular bayes net](http://i.stack.imgur.com/IriH8.jpg) Here are the histogram plots of the distribution of the probabilities for A and B. (breaks are set to 20). ![histogram plot](http://i.stack.imgur.com/ruejz.jpg) Here is the scatter plot of the probabilities of B vs A. ![scatter plot](http://i.stack.imgur.com/zYfCS.jpg)


title: Kohonen self organizing maps: determining the number of neurons and grid size

text: I have a large dataset I am trying to do cluster analysis using SOM. The dataset is huge (~ billions of records) and I am not sure what the number of neurons should be or the SOM grid size to start with. Any pointers to some material that talks about estimating the number of neurons and grid size would be greatly appreciated.


title: Fitting a Generalized Linear Model (GLM) in R

text: I am learning about Generalized Linear Models and the use of the R statistical package, but, unfortunately, I am unable to understand some fundamental concepts. I am trying to develop a GLM - Poisson model but using a specific log link function. The function is of the form $$\ln(E(y_i)) = \ln(\beta_1) + \beta_2 \ln(\text{exp}_1) + \beta_3 \ln(\text{exp}_2).$$ In this equation, $\text{exp}_1$ and $\text{exp}_2$ are measures of exposure in the model. From my understanding, in R, I would first load all the data and ensure it was properly set-up. I then believe I should be running:               model = glm(formula = Y~exp1+exp2, family=poisson(link="log"),data=CSV_table)      As I am new to GLMs and R, I am not exactly sure what specifying poisson(link="log") does. I hope this question isn't too trivial. I have been trying to google clear concise explanations online for hours; however many answers/links assume a level of knowledge higher than mine.


title: Early papers using statistical models for word completion?

text: I am trying to find references about (early, as in pre-2000, but i would really love something pre-1990) statistical language models in word completion (like T-9, or Google's search autocomplete), but find it hard to find such papers. The best match i found was using statistical models to support semi- automated translations, but are there also works that went into that area "just to make (word processing) life a little easier"?


title: Calculating Effect Size in R

text: I've been trying to learn about effect size in relation to linear regression and am wondering how to implement it in R. Sure, I have p-values and they indicate how "predictive" the explanatory variable is. However, for each variable in a linear model, I was wondering how to compute a standardized score for how much it impacts the response variable. Some sample data if you can present a R solution.               x1 = rnorm(10)     x2 = rnorm(10)     y1 = rnorm(10)          mod = lm( y1 ! x1 + x2)     summary(mod


title: Help with determining model for cross-sectional data where all variables are dummies

text: I am currently working on my dissertation project where my data are essentially all dummies. From my dependent to my independent variables, everything is a dummy variable (0,1), at least for the preliminary model. The data are cross-sectional. I was wondering what the best model to use would be. I have decided to do a simple regression followed by logit and probit. Does that sounds like it covers everything? I have some interaction terms too. Also, what tests I can use?


title: Finding area of k-means

text: My problem is that I have different points (4) in a 3D space each point allows to a class. Other 4 points in a 3D space each point allows to a class. I want to know which of these spaces is more sparse. My question are:   1. How I can to get the area of svm with 2-dimensional features?   2. How I can to get the volume of svm with 3-dimensional features?


title: Deriving the bivariate Poisson distribution

text: I've recently encountered the bivariate Poisson distribution, but I'm a little confused as to how it can be derived. The distribution is given by: $P(X = x, Y = y) = e^{-(\theta_{1}+\theta_{2}+\theta_{0})} \displaystyle\frac{\theta_{1}^{x}}{x!}\frac{\theta_{2}^{y}}{y!} \sum_{i=0}^{min(x,y)}\binom{x}{i}\binom{y}{i}i!\left(\frac{\theta_{0}}{\theta_{1}\theta_{2}}\right)^{i}$ From what I can gather, the $\theta_{0}$ term is a measure of correlation between $X$ and $Y$; hence, when $X$ and $Y$ are independent, $\theta_{0} = 0$ and the distribution simply becomes the product of two univariate Poisson distributions. Bearing this in mind, my confusion is predicated on the summation term - I'm assuming this term explains the the correlation between $X$ and $Y$. It seems to me that the summand constitutes some sort of product of binomial cumulative distribution functions where the probability of "success" is given by $\left(\frac{\theta_{0}}{\theta_{1}\theta_{2}}\right)$ and the probability of "failure" is given by $i!^{\frac{1}{min(x,y)-i}}$, because $\left(i!^{\frac{1}{min(x,y)-i!}}\right)^{(min(x,y)-i)} = i!$, but I could be way off with this. Could somebody provide some assistance on how this distribution can be derived? Also, if it could be included in any answer how this model might be extended to a multivariate scenario (say three or more random variables), that would be great! (Finally, I have noted that there was a similar question posted before (Understanding the bivariate Poisson distribution), but the derivation wasn't actually explored.)


title: What does Class Complexity mean in Weka?

text: When running Weka on my dataset, in the results printout I get the following rows:               Class complexity | order 0           44138.7117 bits      8.9205 bits/instance     Class complexity | scheme            46697.7417 bits      9.4377 bits/instance     Complexity improvement     (Sf)      -2559.03   bits     -0.5172 bits/instance      What does Class complexity mean? And what is the difference between "order 0" and "scheme"? Does it have to do with entropy and information gain?


title: Great model prediction, but no significant variables

text: I am performing a binary logistic regression. I have developed a simple model which I am testing using the SPSS application. This first determines the predictive ability of a baseline model without using the independent variables and then adds in the variables to eatablish if the "new" model increases the % prediction. The baseline model prediction % is already at 94.7%. The model with variables added pushes this up to 98.6%. This is great, but none of the preditive variables is significant which confuses me. Can you offer an explanation for this lack of significance?


title: Foundamental limit theorem of Markov chains with higher order chains?

text: I have a second order Markov chain with 4 states {A,T,C,G} (the 4 DNA nucleotides). the transition matrix looks like this:                   A    T    C    G     AA[0.1, 0.6, 0.2, 0.1]     AT[0.3, 0.1, 0.5, 0.1]     AC[0.5, 0.3,  0,  0.2]     AG[..., ..., ..., ...]     TA[..., ..., ..., ...]     TT[..., ..., ..., ...]     TC[..., ..., ..., ...]     TG[..., ..., ..., ...]     CA[..., ..., ..., ...]     CT[..., ..., ..., ...]     CG[..., ..., ..., ...]     GA[..., ..., ..., ...]     GT[..., ..., ..., ...]     GC[..., ..., ..., ...]     GG[..., ..., ..., ...]      I wanted to calculate the stationary probability vector for the 4 states to which this matrix converges. The Markov chain is regular. In case of first order Markov chains this is easily done by calculating the limit of $P^n$ with $n\rightarrow \infty$. I do not know how to approach the problem in case of second order Markov chains. Also, having a limited dataset from which to determine the transition matrix, can I consider the stationary distribution of the 4 nucleotides as being the theoretical distribution I would have if I had a much larger pool from which to draw (with the same transition matrix)? In other words, can I consider the stationary distribution like an estimation of the theoretical nucleotide frequency given the transition matrix obtained from limited data?


title: Test whether time of maximum differs across two groups

text: Suppose I have 10 repeated measurements of a quantitative variable $X$ of the individuals within two groups (let's call them A and B). For each group, I have the estimated means of variable $X$ at the 10 time points, that is, $$\hat{\mu}_1^A, \hat{\mu}_2^A, \ldots, \hat{\mu}_{10}^A$$ and $$\hat{\mu}_1^B, \hat{\mu}_2^B, \ldots, \hat{\mu}_{10}^B,$$ and the corresponding $10 \times 10$ variance-covariance matrices of these estimates. Tests of $$H_0: \mu_t^A = \mu_t^B$$ (for $t = 1, \ldots, 10$) can be easily conducted using Wald-type tests (with an appropriate correction for multiple testing) with: $$z = \frac{\hat{\mu}_t^A - \hat{\mu}_t^B}{\sqrt{Var(\hat{\mu}_t^A) + Var(\hat{\mu}_t^B)}}.$$ Now I want to test whether the time point of the maximum within each group differs across the two groups. I suppose the null hypothesis could be written as $$H_0: t = t' \; \text{for} \max({\mu}_t^A) \; \text{and} \max({\mu}_{t'}^B).$$ What would be a way of testing that? **EDIT** : Some R code to simulate data to play with.               library(MASS)          ### number of individuals in the two groups     n1 <- 50     n2 <- 50          ### assume a flat trajectory with a single peak at times 4 and 6 in the respective groups     mu1 <- rep(0,10)     mu2 <- rep(0,10)     mu1[4] <- 1     mu2[6] <- 1          ### make the raw data autocorrelated with an AR1 structure     V1 <- 5 * toeplitz(ARMAacf(ar=.7, lag.max=9))     V2 <- 5 * toeplitz(ARMAacf(ar=.7, lag.max=9))          ### simulate raw data     X1 <- mvrnorm(n1, mu1, V1)     X2 <- mvrnorm(n2, mu2, V2)          ### observed means over time in each group     m1 <- apply(X1, 2, mean)     m2 <- apply(X2, 2, mean)          ### find time of the maximum mean in each group     tmax1 <- which(m1 == max(m1))     tmax2 <- which(m2 == max(m2))          ### plot means over time in each group     plot(1:10, m1,  type="o", pch=19, ylim=range(c(m1,m2)), col="red", xlab="Time", ylab="Mean")     lines(1:10, m2, type="o", pch=19, ylim=range(c(m1,m2)), col="blue")     points(tmax1, max(m1), pch=19, cex=2, col="red")     points(tmax2, max(m2), pch=19, cex=2, col="blue")      And an example of what the resulting plot looks like. ![figure of means over time in two groups](http://i.stack.imgur.com/58eiB.png)


title: Sample size calculation for paired t-test single arm using r

text: I am planning a study on alteration in AUC of a drug due to a surgery of gastrointestinal tract. It's a pre-post study. A group of patients will undergo blood sampling for AUC (area under curve) both before and after surgery. How can I calculate sample size using R?


title: Bootstrapping importance sampling estimates

text: I've got two algorithms, $A$ and $B$, I want to evaluate. Both algorithms differ by the distribution of points they produce. Let us suppose that $A$ generates points from $P_S$ distribution while B generates points according to $P_T$ distribution. Further, I've got a dataset of 21,000 points generated by algorithm $A$, i.e. from $P_S$ distribution. Each point has it's value $v_i$ in {-1, 0, +1}: $v_i \in \{-1, 0, 1\}$. Next, I am generally interested in studying two quantities: (1) mean value over points the algorithms produce, (2) once deployed, what is probability that algorithms $A$ and $B$ will result in positive mean after observing $k = \{100, 200, ..., 20000\}$ points? In order to answer (1) I run importance sampling. Since the points were generated by A, $mean_A$ is calculated straightforwardly $$ mean_A = \frac{1}{N} \sum_i v_i \frac{P_S}{P_S} $$ For $B$ it's a little bit more tricky: $$ mean_B = \frac{1}{N} \sum_i v_i \frac{P_T}{P_S} $$ It turns out $B$ outperforms $A$, i.e. it puts more weight on positive points. The problem arises then I approach the question (2). Naive idea was to bootstrap $k$ points from the dataset uniformly ( = from $P_S$) and run the same importance sampling re-weighting as before. And the results are discouraging: for any $k$, no matter how many bootstraps I run, the importance sampling estimates for $B$ are positive less frequently, than for $A$, i.e. $A$ outperforms $B$, despite that $B$ puts more weight on positive points than $A$ (as proved by experiment (1)). A possible problem is that for a considerable part of the points $P_T = 0$, i.e. they are never produced by $B$. In the literature I examined, it is claimed that this is not a problem for IS, but obviously this is a problem for naive bootstrapping I do: adding zero-weighted points will force $B$ to perform worse, since they are more and more likely to fill in the bootstrap sample. Moreover, putting infinite weight in a single positive point will again result into $B$ under-performing $A$. The only solution I see so far is to filter out all the "impossible for $B$" points before running bootstrapping and calculate IS estimate with respect to source distribution "$P_S$ with the impossible points removed". Is this the only possible solution? Is the problem considered somewhere in the literature? Is there anything I can reference to? Or, should I replace this uniform bootstrap with something else (non-uniform)? Any suggestions are appreciated. PS. It seems to me that I can replace the naive bootstrap I'm doing with the sampling/importance resampling algorithm (SIR): given $N$ points sampled from $A$ (the dataset) I can calculate importance weights as in IS and use them (after normalisation) as probabilities to sample again. The resulting sample is going to be iid sample from $B$ and thus can be (?) used as a bootstrap sample from $P_T$. Am I right?


title: In R, coefficients of MA function are wrong?

text: I'm currently sifting through my copy of Analysis of Financial Time Series 2nd Edition by Ruey Tsay, and one of the sections involves fitting a MA model to certain data (data set is here). Here's the fit with exact maximum likelihood according to the text, with certain insignificant parameters removed: rt = 0.013 + a(t) + 0.181a(t−1) − 0.121a(t−3) + 0.122a(t−9) σ(a) = 0.0724 However, when I try to fit it with R...               > mew = read.table("m-ew.dat")     > arima(mew,order = c(0,0,9),fixed = c(NA,0,NA,rep(0,5),NA,NA),method = "ML")     Call:     arima(x = mew, order = c(0, 0, 9), fixed = c(NA, 0, NA, rep(0, 5), NA, NA),      method = "ML")          Coefficients:             ma1  ma2      ma3  ma4  ma5  ma6  ma7  ma8     ma9  intercept           0.180    0  -0.1318    0    0    0    0    0  0.1373     0.0132     s.e.  0.031    0   0.0362    0    0    0    0    0  0.0327     0.0029          sigma^2 estimated as 0.005282:  log likelihood = 1039.1,  aic = -2068.21      As you can see, the ma1 coefficients are the same, but ma3 and ma9 are different, even with method = "ML", i.e. maximum likelihood. Why is this? Also, from a practical standpoint, while ma2 and ma4-ma8 may be 0 (their 95% confidence intervals overlap with 0), removing them from the model raises the AIC, lowers the p-value with regards to the Ljung-Box test on the residuals, and also lowers the log-likelihood value. Is it even worth removing these parameters if such things happen?


title: comparing same negative binomial models built with subset data

text: I have a dataset made our of several stacked datasets (one for each state). I want to check whether a zero inflated negative binomial model with data from an individual state is different from the ZINB made from the whole data (national). Are there any suggestions on how to do a test like this? One idea I had was to calculate the likelihood of the national model using only the single state data and then compare them, but the models would not be nested so I would be left with zero degrees of freedom. Thanks for any help,


title: Example of a subset of $\mathbb{R}$ which is not measurable?

text: Can anyone please give me an example of non-measurable real set?


title: PROCESS Hayes | Interaction Insignificant | 2/3 Conditions Significant

text: When using PROCESS by Hayes I get a insignificant interaction effect but 2/3 of the conditional effects are significant. Is there still something I can say about this in my results? So for instance, there is no overall interaction significance but there is on the low levels? ![Output](http://i.stack.imgur.com/sR5u8.png) Also, when I plot the output I see the following (I'm looking if regulatory institutions moderate the effect of labour cost on FDI): ![Plot](http://i.stack.imgur.com/HM4cP.png) I see that for high labour cost the regulatory institutions do matter. And not for low labour cost. So can I say that for in case of high labour cost the regulatory institutions do matter? Or because the overall interaction effect is significant I can't say anything about it?


title: How could a t-test be used for comparing two ML algorithms?

text: I am new to ML, and I am reading a paper about ML comparison. It underlines the comparison between two algorithms based on the t-test but it does not describe the process further. I know what a t-test is, but I don't understand how to implement it for ML algorithm comparison.


title: Poisson is to exponential as Gamma-Poisson is to what?

text: A Poisson distribution can measure events per unit time, and the parameter is $\lambda$. The exponential distribution measures the time until next event, with the parameter $\frac{1}{\lambda}$. One can convert one distribution into the other, depending on whether it is easier to model events or times. Now, a gamma-poisson is a "stretched" poisson with a larger variance. A Weibull distribution is a "stretched" exponential with a larger variance. But can these two be easily converted into each other, in the same way Poisson can be converted into exponential? Or is there some other distribution that is more appropriate to use in combination with the gamma-poisson distribution? The gamma-poisson is also known as the negative binomial distribution, or NBD.


title: Is it right to use a 2-way ANOVA for 3 experimental and 2 control groups? And should the controls be included in the ANOVA?

text: I'm currently working on our thesis and I got stuck with the analysis part for I am not sure which statistical test I should use. My work had 5 groups. 1 positive and 1 negative control and 3 experimental groups. Is it right to use a 2-way ANOVA? And should the controls be included in the ANOVA? I'd like to compare if the experimental groups significantly varied from the control groups. What would be the best fitting test? Thank you so much.


title: Iterative maximum likelihood estimation of history

text: The model is the following:   * You receive a hidden integer h from 1 to N (uniform distribution)   * You take an action a numbered between 1 to M. Describing the entire model is a MxN probability matrix with each entry saying "The probability of taking action a given hidden number h" The rows of this model all sum to 1 (given some hidden number h the probability of taking some action is 1). An example of what I had in mind would be the following. You receive a poker card from a deck containing (Ace, King, Queen, Jack, Ten, 9, 8, 7, 6, 5, 4, 3, 2) which we number from 1 to 13 (N=13). The available actions are Fold, Call, Raise (so M=3). In the model generating the history we receive, there might be properties like P(Raise | Ace)=1 or P(Call | 2)=0. Now we get a large amount of history of seen actions. Sometimes we get to see the hidden number generating the action, most often we do not. How do we use an iterative maximum likelihood estimation tool to solve this problem? The presented is a sub-problem. In the real problem the action seen from this model (and the probability distribution over h it generates) is the input to another model of the same kind as this. However, if I figure this out I'm certain I will be able to solve the rest. I've attempted to use gradient ascent with lagrange multipliers to handle the constraints, however I'm screwing up somewhere and cannot figure it out. I'm open for any advice and suggestions. Thanks in advance.


title: How to deal with unbalanced group sizes in mixed design analysis?

text: I have 2 x 2 x 2 mixed design with two between subject factors (sex, organizer) and one within subjects factor (task). The group sizes of the 'sex'-factor is unequal. When I perform a factorial repeated measures ANOVA, I get the following warning message: `Warning: Data is unbalanced (unequal N per group). Make sure you specified a well-considered value for the type argument to ezANOVA().` I used the following model in R: `model <\- ezANOVA(data=df, dv=.(top_start), wid=.(id), between=.(sex, org), within=.(task), type = 3, detailed = TRUE)` I used `type = 3` for the Anova because, as I understood, it is suited for unbalanced group sizes. I have the following questions:   * Do I need to code contrasts when all my factors have only two levels?   * Did I use the right type of Anova?   * Are there other ways to do this analysis?


title: Most effective way to learn time series with poor quant background

text: End goal will be practical application (model building) by using time series analysis to analyze/forecast macroeconomic/finance data. Background: I have taken stats, introductory econometrics, intermediate macro, micro in the past couple years, but is has been 10 years since I took my only calculus course. I'm confused whether I should (a) go do all the required math, (b) just find a math econ refresher course, then to read a book focused on theory or (c) just try to build models. I realize I'm starting from a low base, but I'm hoping to become reasonably proficient in a year (unless this is unrealistic?). Thanks in advance for your help.


title: Katz backoff for n-gram language models

text: I have defined unigram, bigram and trigram (language) models with my training data. Now I am checking their fitness towards a test data by using Katz backoff with a fixed discount of 0.5 (mo Good-turing). If I start with a trigram, when there is a zero count, the trigram falls back to bigrams and if the bigram also has a zero count, the bigram falls back to an unigram (all unigrams defined in the test data). Each time I fall back, I multiply the (n-1) data with a parameter alpha. If I start with a bigram, I may fall back to an unigram. That is, I get results for trigrams and bigrams that are comparable. My question is : How should I process the unigrams, if I am only analyzing the unigrams? I have defined a closed vocabulary in a such a way that all tokens in the test text are defined in the training text.


title: A good way to show lots of data graphically

text: I'm working on a project that involves 14 variables and 345,000 observations for housing data (things like year built, square footage, price sold, county of residence, etc). I'm concerned with trying to find good graphical techniques and R libraries that contain nice plotting techniques. I'm already seeing what in ggplot and lattice will work nicely, and I'm thinking of doing violin plots for some of my numerical variables. What other packages would people recommend for displaying a large amount of either numerical or factor-typed variables in a clear, polished, and most importantly, succinct manner?


title: Minimal I-maps induced by sets of scopes. Clarification needed

text: I have a question about Prop. 9.1 on page 307 in "Probabilistic Graphical Models" (link to google books) (Koller / Friedman). I don't see why $\mathcal{H}_{\Phi}$ is a minimal I-map: If I have only one factor in $\Phi$ over two Variables $X$ and $Y$ and define it in a way $X$ and $Y$ are independent; then $\mathcal{H}_{\Phi}$ will be a 2-node clique (?). But that's not a minimal Map for the defined distribution (which would be 2 isolated nodes). Would be pleased if anyone could clarify at what point my thinking is wrong.


title: Why are Gaussian "discriminant" analysis models called so?

text: Gaussian discriminant analysis models learn $P(x|y)$ and then apply Bayes rule to evaluate $$P(y|x) = \frac{P(x|y)P_{prior}(y)}{\Sigma_{g \in Y} P(x|g) P_{prior}(g) }.$$ Hence, they are generative models. Why then is it called discriminant analysis? If it is because we finally derive a discriminant curve between the classes, then that happens for all generative models.


title: Multiplying Gaussian distributions of different dimensions

text: The multiplication of multivariate Gaussian distributions defined over some parameter vector of a given dimension can be achieved by the following. Assuming that the Gaussian is parametrized by the precision matrix ($\Lambda$) and mean ($\mu$). We can compute the new precision as: $$ \Lambda = \Lambda_1 + \Lambda_2 $$ Also, the new mean can be computed using the identity: $$ \Lambda \mu = \Lambda_1 \mu_1 + \Lambda_2 \mu_2 $$ So, this is very convenient. However, I have a situation where I need to perform such a multiplication where one of the Gaussian distribution depends on $n$ number of variables. However, the second distribution depends on a subset of these $n$ variables. So, to perform this multiplication, I need to somehow make the second distribution also depend on all the $n$ variables, though in reality only a subset of these $n$ variables are affecting the underlying distribution. The exact situation is that I have a posterior distribution $q(\theta)$ where $\theta$ is n-dimensional and I have the likelihood of the following form: $$ L(\theta) = \prod_{i=1}^{k}L(\theta_i) $$ Here, $L(\theta_i)$ is a subset of $n$ variables and $L(\theta_i)$ is independent of $L(\theta_j)$ for $i \neq j$. I need to multiply $q(\theta)$ with a $L(\theta_i)$. So I somehow need to express this $L(\theta_i)$ as a function of $\theta$. Each $L(\theta_i)$ is a 3-dimensional normal distribution i.e. $\theta_i$ is 3-dimensional. Does anyone know how this can be done?


title: Does a high likelihood of duplicate samples invalidate my data or certain operations on it?

text: We have data from random polling within a group of people, sampled over a long period of time. As the pool is relatively small and anonymity is crucial there is a high likelihood that we have duplicate samples (though we can't tell for certain because we'd expect responses to vary over time). Does this somehow invalidate the data, or mean that certain operations on it will not be meaningful? Or is it OK to proceed as normal so long as I state this likelihood? * * * I'm happy to provide more detail if needed - just say what would help in a comment. Also I'm making this CW. Please feel free to edit the question if there are other relevant implications of duplicate data that would be worth specifying.


title: Mixed effects log-linear models

text: There are occasions where I would like to fit log-linear models where the independence assumption between observations is violated. It is the normal case that I have multiple observations from each subject. Is there a mixed- effects extension to log-linear models like there are for linear and generalized linear models? I'm interested in the existence of such a thing in principle, but also its concrete existence in an R library.


title: How can I display empirical pdf of my 100x1 vector data in Matlab?

text: I have a data which is 100x1 vector. How can I display its empirical pdf in Matlab? Also, if I want to compare the pdf of three vectors on the same graph, then how to do that? Right now I am using pdfplot.m file to plot my empirical pdf, however when I want to compare the 3 distributions by using 'hold on', then firstly its not working and secondly all the distributions are in same color. Thanks! Also, I don't want to plot cdf or histogram.


title: Hierarchical classification

text: I'm currently working on the classification with massive amount of data. Similar to the kaggle one. Data input consist of features and multiple labels that can be hierarchically aligned. At first I flattened the data and tried to learn multi-label classifier. That involved both, methods that are specially designed for it (like multi-label kNN) and One-vs-All methods. However these approaches didn't yield much result and were really complex. Afterwards I found structured SVM library. Description says, that this library actually minds the structure in the data. However due to its complexity and lack of examples, I didn't have enough time to try it myself. I was thinking maybe someone could shed some light to other hierarchical classification methods out there (in R, Python, C or Java), especially the guys that were tackling that problem on kaggle. What's your approach to this?


title: Cross-section and fixed effects models

text: Imagine 2 models: 1 for production, 1 for revenue, have cross-section and time dummies. Is it permissable to use CS in one and time in the other? How do I justify this? F-test for fixed effects is indicating this and also dummies I want to omit in each case are not really significant.


title: Understanding replication and lack-of-fit in regression modeling

text: Regarding regression analysis, I once heard that "Replicate at various predictor values will allow for lack-of-fit analysis". I do not understand how to get this kind of statement?


title: What is the difference between a nested and non-nested model in CFA?

text: I thought it had to do with varying numbers of factors but based on the literature I have read it seems to do more with fixing/freeing of parameters. Any info would be greatly appreciated!


title: How SAS calculates PROC GLM

text: I'm trying to do a pairwise test to calculate p-values between various products that are being evaluated by several assessors. I'm trying to reproduce what SAS is doing in R, and I'm running into an issue in one particular situation: when the items being compared (in my case SAMPLE) are complete for one item in the pairwise comparison but not the other. Here's an example:               data raw;     INPUT JUDGE REP SAMPLE X1;     datalines;     1 1 1 3.0     1 2 1 4.0     2 1 1 7.0     2 2 1 2.0     1 1 2 1.0     1 2 2 6.0     2 1 2 3.0     ;          run;               proc glm data = raw;     class judge sample rep;     model      X1 = judge sample rep judge * sample / ss3;          lsmeans sample / stderr tdiff alpha=0.10 e=judge*sample;          run;      In this case we have a JUDGE (2) who is missing REP 2 for SAMPLE 2. SAMPLE 1 has a complete set of data. When I calculate the pairwise comparison between SAMPLE 1 and 2, I get a p-value of 0.48973, whereas in SAS it's an even 0.5. To get my value I'm doing the following: mse = interaction_sum_of_squares / degrees of freedom = 0.66666666 / 1   constant = sum of (1/n1j + 1/n2j) = (1/2 + 1/2 + 1/2 + 1) = 2.5   sd = sqrt( (mse / num_judges^2) * constant)) = 0.6454972 t-value = lsmean1 - lsmean2 / sd = (4.0 - 3.3333333) / (0.6454972) = 1.032796 Meanwhile, the SAS t-value is 1.0. I'm matching the lsmeans that SAS is getting, and the interaction_sum_of_squares and the degrees of freedom match what SAS is putting out as well. That really only leaves the constant, or something else SAS is doing which I can't figure out. I match SAS perfectly when the data is balanced, it's only when it's unbalanced that it's not working.


title: Time series and cross correlation

text: I have two time series:   1. numbers of subscriptions per hours for 10 days, for a website.   2. the time where TV ads, for a specific TV channel, have been displayed e.g. 10:25AM on Monday 12 of August, TV AD1 displayed on CNBC as an example. The issue with that time series is there is only the time and the date but nothing else. I would like to do a cross correlation to see the impact, or the lack of impact, of the TV ads on the numbers of subscriptions. But with only the time of the adverts, is this cross correlation doable? If yes, how? I'm not asking for a specific solution designed by any of you, rather website links where I can find some ideas or any hints as I'm stuck at the moment.


title: In maximum likelihood, pattern matrix, one case $>1.0$

text: How bad is this for my model? I've used ML (maximum likelihood) with promax rotation and Kaiser on a dataset ($n>300$) to extract underlying factors to build a model. Eigenvalue selection was reset from $.1$ to $.5$, so as to focus on cases (questions) that explain the most covariance according to underlying factors. I hope my explanation makes sense, as I`m new to this. Anyways, 9 clean factors emerged, with only this one troubling number, in the pattern matrix. I tested reliability of the factors with ANOVA and reliability was excellent (most factors above $.900$). I`ve already used the information to refine and collect data with a questionnaire, but this one number is still troubling me. Does it invalidate the model? Is it a serious flaw in my model? What are the implications of a case number in the pattern matrix being $>1.0$? None of the communalities were $>1.0$, so I don`t think this is a Heywood case. This is potentially a serious problem, as it might necessitate me redoing the model and adding a year to my studies - the questionnaire is a longitudinal one, taking one academic year. Also, I am interested in publishing the resulting questionnaire to aid in others` related research, but if this one number is a serious error, then publication is out of the question. I would greatly appreciate some advice regarding this problem.


title: Testing equality of relative frequency ratio in binomial sample

text: I am sorry if this question is too basic or was somewhere else here on Cr- Validated, but i didn't find anything suitable. I study patients with cancer who were made a byopsy taking to test presence of cancer cellcs 12 prostate core samples. When core samples are explored, they are divided into two groups whether there was found 1-2 prostatic intraepithelial neoplasia samples or 3-4 prostatic intraepithelial neoplasia samples (or PIN samples). I want to test the hypothesis that cancer detectability ( 0 states for no cancer cells found and 1 states for presence of cancer cells) is higher among patients with 3-4 PIN samples rather than 1-2. To visualise - this is how sample could look like               PIN Group Cancer     1   I     0     2   I     0     1   I     1     ...     4   II    1     3   II    1     4   II    0     ....      Mathematically speaking I have to test hypothesis of relative frequency ratio equality. How do I do that? Simultaneously I opened two sources - one asserting that I have to calculate two CI for ratios and if they intersect then it means that ratios can be equal, and other source claiming it is wrong method. What do you think? I will be eternally grateful for any clarifications!


title: AIC and likelihood function applied to a dataset

text: **Update: First part answered. Directly jump to the second part for those wanting to answer!** **First Part** I am trying to have a better understanding of what is a likelihood function and what is AIC. Here is my current understanding: If we toss a coin 11 times and get 5 heads and 6 tails, the likelihood function is $$L(p) = P(data|p) = p^5(1-p)^6$$ p being the hyopthesis of the function which is the probability to get head. the Maximum Likelihood is the value of p that satisfies this equation: $$\frac{dL(p)}{dp} = 5p^4(1-p)^6-6p^5(1-p)^5=0$$ Therefore, $\hat p=\frac{5}{11}$ The AIC is given by: $$AIC=2k-2ln(L)$$ k being the number of parameters (one in our example). Therefore: $$AIC=2*1-2*ln(0^5*(1-0)^6) = 2-2*ln(L)=2$$ I guess I missunderstand something! Is there a mistake? Where is it? **Second Part** I am also wondering how does this AIC calculations apply when performing regressions. For example we have this simple set of data (writing in R format):               respond_variable = c(1,2.2,3.1,4)     predictor_variable = c(1,2,6,10)      We perform a linear and a quadratic model on them. How is the AIC calculated? Thanks a lot for your help!


title: libSVM - probability estimates in multi-class problems

text: I am using libSVM for classification on a 3-class dataset. Using the option "-b 1" - for getting probability estimates for prediction - gives me surprising results. I am reproducing the issue with a much simpler dataset here. The option "-b" can take only two values: "-b 1" denotes we need probability estimates of the predicted label in the output, and "-b 0" indicates we don't need the estimates. "-b 0" is the default i.e. not specifying an option is equivalent to saying "-b 0" This is my sample training data (let's say the file is called _simple_ ):               1 1:0.1 2:0.1 3:0.1     1 1:0.15 2:0.15 3:0.15     2 1:0.5 2:0.5 3:0.5     2 1:0.55 2:0.53 3:0.49     3 1:0.9 2:0.92 3:0.93     3 1:0.88 2:0.91 3:0.97      It's easy to see what I am doing:   * for vectors with each dimension ~0.1, the label is 1.    * for vectors with each dimension ~0.5, the label is 2.    * for vectors with each dimension ~0.9, the label is 3. Here's my test data (the file is called _simple.t_ ):               1 1:0.1 2:0.13 3:0.11     2 1:0.49 2:0.55 3:0.56     3 1:0.9 2:0.95 3:0.99      Commands run with probability enabled:               ./svm-train  -b 1 simple     ./svm-predict -b 1 simple.t simple.model output     Accuracy = 0% (0/3) (classification)      Output file:               labels 1 2 3     3 0.0447161 0.226854 0.728429     1 0.49332 0.248142 0.258538     1 0.713506 0.24226 0.0442344      Commands run with probability disabled:               ./svm-train  simple     ./svm-predict simple.t simple.model output     Accuracy = 100% (3/3) (classification)      Output file:               1     2     3      I find this very surprising, if not absurd! Why is enabling probability estimates changing the way the classifier works? _Drastically_ : the accuracy drops from **100%** to **0%**. What am I missing here?


title: Comparing two data sets for similarity with variable data quality

text: I have to compare two data sets for similarity. I'm designing a piece of software which will help improve a drummers accuracy. The software displays a rhythm which the drummer then attempts to play. The first set of data is the reference set of the timings of beats which have been played. The second set represents the beats detected by the user's drum playing. If the user played perfectly the two sets of data would be identical. I need to find a way to score the user's playing as a percentage. i.e. how accurate were their timings. Ideally I'd just compare the two sets of data one by one however this doesn't work because the beat detection algorithm I'm using isn't 100% reliable. Sometimes it misses a beat and other times it records spurious beats. This means that the two sets of data can be different lengths. Example:               1) Refr: 0.000000     2) Refr: 0.187500     3) Refr: 0.375000     4) Refr: 0.562500     5) Refr: 0.750000     6) Refr: 0.937500     7) Refr: 1.125000     8) Refr: 1.312500          1) Beat: 0.000000     2) Beat: 0.208984     3) Beat: 0.374745     4) Beat: 0.394745     5) Beat: 0.580505     6) Beat: 0.766266     7) Beat: 0.940414     8) Beat: 1.149391     9) Beat: 1.335152      Is there a robust statistical technique for this kind of problem?


title: Probability of reaching node A from node B in exactly X steps

text: I have a three-node matrix with two edges (A-B and A-C). I would like to determine what the probability is starting from B and ending at C in exactly 100 steps. I have only written out probabilities:               P(A|B) = 1     P(B|A) = 0.5     P(A|C) = 1     P(C|A) = 0.5      But there are so many combinations of ways to get from B to C in exactly 100 steps using these probabilities. Any suggestions on how to continue this problem?


title: Finding critical region of statistic when the distribution of the transform of the statistic is known

text: Let $X_1, ..., X_n$ be a random sample from $\mathrm{Exp}(\lambda_1)$ and $Y_1, ..., Y_m$ be a random sample from $\mathrm{Exp}(\lambda_2)$. I have found the generalized likelihood ratio test (GLRT) statistic $\Lambda$ for testing $H_0:\lambda_1 =\lambda_2$ vs $H_1:\lambda_1 \neq \lambda_2$ is $\Lambda=\frac{n^n m^m}{(m+n)^{m+n}} T^{-n}(1-T)^{-m}$ where $T = \frac{\sum_i X_i}{\sum_i X_i + \sum_j Y_j}$. Now I need to find the critical region of the test by the asymptotic distribution of the test $\Lambda$. I know the relation $W = 2 \log(\Lambda) \sim \chi_d^2$ where d is the number of parameters to be tested. How can I use this relation to find the critical region in the above problem?


title: Guess level of measurement of unknown data

text: i'm writing a few convenience functions for exploring datasets. For that, I have to guess the type of data from the values. What I'm looking for are some good rules of thumb. Some easy cases:   * `['spam', 'eggs']` is nominal, because the dataset contains strings.   * `[true, false, false]` is nominal, too.   * `[0, 1, 0, 0, 0, 1, 1]` is either nominal or ordinal.   * `[0, 1, 4, 3, 1, 2, 2]` too, but it's probably ordinal. Many values, and they have numbers, but not too many.   * `[0.2, 2, 1.6, 3]` is cardinal, as the level of precision suggests there could possibly be many other values. As you see the judgment is absolutely vague, but that's not the point: I just need heuristics that work in 80% of the cases. Or better, 95%. More ideas?


title: What post-hoc test can be used on data with missing values?

text: I have run 13 data-stream classification approaches on 8 different datasets, giving a 8x13 result table. I was hoping to use the Friedman test to test for a significant difference between the approaches. However one of the approaches can not be used on two of the datasets so there are two missing values in my result table. I have read that the Skillings-Mack statistic can be used when there a missing values so I used the Matlab implementation I found at: http://www.geog.ucsb.edu/~pingel/code/mackskill/skillmack.html which gave a T value of 48.58 and a p-value < 0.0001. My question is in regards to post-hoc tests. Are there any recommended post- hoc tests for data with missing values or can any of the usual post-hoc test be used?


title: How to 'sum' a standard deviation?

text: I have a monthly average for a value and a standard deviation corresponding to that average. I am now computing the annual average as the sum of monthly averages, how can I represent the standard deviation for the summed average ? For example considering output from a wind farm:               Month        MWh     StdDev     January      927     333      February     1234    250     March        1032    301     April        876     204     May          865     165     June         750     263     July         780     280     August       690     98     September    730     76     October      821     240     November     803     178     December     850     250      We can say that in the average year the wind farm produces 10,358 MWh, but what is the standard deviation corresponding to this figure ?


title: Visualizing PCA time series

text: I have the following experiment. My predictors are: Plant area, 6 types of plants, 5 time points and 12 replicates for class. Time points are not equally spaced. A normal plot of Area vs. time will show $6\times 12\times 5$ points in the plot which doesn't reflect the fact that 5 of those points correspond to the same plant at different time point but whose Area has changed with time. Is there a plot that shows some sort of derivative approach?


title: Problem of "clustering" into most similar groups

text: I do not think that the following problem can be solved with k-means clustering. I am not sure though. Okay, let me describe the problem. I need to find a way or an algorithm that groups members of a given data set (of positive integers) so that the difference between group means is **minimized** (not maximized, as usual). There are two constraints:   1. The number of groups should not exceed log(N) base 2. N is the input array size. Let us assume that N = 16, always.   2. The size of the group should be at least log(N) base 2. Here it is 4. Please see below example. Guidance toward a solution to this problem is appreciated.                  for input array = (12, 14, 16, 16, 18, 19, 20, 21, 24, 26, 27, 29, 29, 30, 31, 32)      One of the solutions could be the following:               set    group           mean     1   (12, 14, 31, 32)    22.25     2   (16, 16, 29, 30)    22.75     3   (18, 19, 27, 29)    23.25     4   (20, 21, 24, 26)    22.75